utils.R

#!/usr/bin/Rscript

## IMPORT
suppressMessages(library(ggsci)) # colors

##############################
# INPUT
read_nanostats <- function(fname){
    df <- read.csv(
        file=fname,
        sep="\t",
        header=TRUE,
        check.names=FALSE,
        stringsAsFactors=FALSE
    )
    df_cols <- c("stat"="Statistic", "value"="Value")
    colnames(df) <- df_cols[colnames(df)]
    return(df)
}

read_quast <- function(fname){
    df <- read.csv(
        file=fname,
        sep="\t",
        header=TRUE,
        row.names=1,
        check.names=FALSE,
        stringsAsFactors=FALSE
    )
    testit::assert(all(colnames(df) %in% names(ASM_TOOL_NAMES)))
    colnames(df) <- ASM_TOOL_NAMES[colnames(df)]
    df <- df[QUAST_VARS, ASM_TOOL_NAMES]
    return(df)
}

read_plasflow <- function(fname){
    df <- read.csv(
        file=fname,
        sep="\t",
        header=TRUE,
        stringsAsFactors=FALSE,
        check.names=FALSE
    )
    testit::assert(all(df$tool %in% names(ASM_TOOL_NAMES)))
    df$tool <- ASM_TOOL_NAMES[df$tool]
    df <- df[df$label %in% names(PLASFLOW_NAMES$labels),]
    df$label <- PLASFLOW_NAMES$labels[df$label]
    return(df)
}

read_rgi <- function(fname){
    df <- read.csv(
        file=fname,
        sep="\t",
        header=TRUE,
        stringsAsFactors=FALSE,
        check.names=FALSE
    )
    colnames(df) <- sapply(colnames(df), function(x){ ifelse(x %in% names(ASM_TOOL_NAMES), ASM_TOOL_NAMES[x], x) })
    testit::assert(all(df$col %in% names(RGI_NAMES$col)))
    df$col <- RGI_NAMES$col[df$col]
    return(df)
}

##############################
# PLOTS

plot_rgi_overlap <- function(df, ctype, col){
    df <- df[df$type == ctype & df$col == col,]
    df_list <- lapply(ASM_TOOL_NAMES, function(x){ df[df[,x] > 0,"label"] })
    names(df_list) <- ASM_TOOL_NAMES[names(df_list)]
    UpSetR::upset(
        data=UpSetR::fromList(df_list),
        # overlap order
        order.by="degree",
        decreasing=FALSE,
        # y-label title
        mainbar.y.label=sprintf("Intersection size (%s hits, %s)", ctype, col),
        # text size
        text.scale = c(1.5, 1.5, 1.5, 1.5, 1.5, 1.5),
        # colors
        set.metadata=list(
            data=data.frame(
                sets=names(df_list),
                Tool=names(df_list),
                stringsAsFactors=FALSE
            ), 
            plots=list(list(type="matrix_rows", column="Tool", colors=ASM_TOOL_COLORS, alpha=0.7))
        )
    )
}

##############################
# CONST

###############
# Assemblers

# names
ASM_TOOL_NAMES <- c(
    "flye"="Flye",
    "megahit"="MEGAHIT",
    "metaspadeshybrid"="metaSPAdes (H)",
    "metaspades"="metaSPAdes",
    "operams"="OPERA-MS"
)
# colors
ASM_TOOL_COLORS <- ggsci::pal_nejm("default", alpha=1)(length(ASM_TOOL_NAMES))
names(ASM_TOOL_COLORS) <- ASM_TOOL_NAMES

###############
# Gene tools
# GENE_TOOL_NAMES <- c(
#     "prodigal_partial"="Prodigal (partial)",
#     "prodigal_total"="Prodigal (total)",
#     "cdhit_unique"="CD-HIT (unique)",
#     "cdhit_total"="CD-HIT (total)"
# )

###############
# CRISPR tools
CRISPR_TOOL_NAMES <- c(
    "minced"="MinCED",
    "casc"="CasC"
)

###############
# PlasFlow

# names
PLASFLOW_NAMES <- list(
    statstype=c(
        count="Sequence count",
        sum="Cumulative sequence length [bp]",
        count_pct="Sequence count [%]",
        sum_pct="Cumulative sequence length [%]"
    ),
    labels=c(
        chromosome="Chromosome",
        plasmid="Plasmid",
        unclassified="Unclassified"
    )
)
# colors
PLASFLOW_COLORS <- list(
    labels=ggsci::pal_nejm("default", alpha=1)(4)[c(2,3,4)]
)
names(PLASFLOW_COLORS$labels) <- PLASFLOW_NAMES$labels

###############
# RGI
RGI_NAMES <- list(
    col=c(
        "Best_Hit_ARO"="ARO term",
        "ARO"="ARO",
        "Drug Class"="Drug class",       
        "Resistance Mechanism"="Resistance mechanism",
        "AMR Gene Family"="Gene family"
    )
)

###############
# QUAST
QUAST_VARS <- c(
    "# contigs",
    "Largest contig",            
    "Total length",
    "N50",
    "L50",                       
    "# N's per 100 kbp"
)