#' print_jams_report(opt=opt)
#'
#' JAMSalpha function
#' @export
print_jams_report <- function(opt = opt, outputdir = NULL, elements = c("readplots", "taxonomy", "SixteenSid", "resfinder", "abricate", "plasmidfinder", "vfdb", "func"), plotucobias = FALSE, report_style = NULL){
if (is.null(report_style)){
report_style <- opt$analysis
}
#See what is available
elementsinopt <- elements[elements %in% names(opt)]
sampletype <- as.character(opt$projinfo[which(opt$projinfo$Run_info == "Run_type"), ]$Run_value)
samplename <- as.character(opt$projinfo[which(opt$projinfo$Run_info == "Sample_name"), ]$Run_value)
#Authorship message
authors <- as.character(as.person(packageDescription("JAMS")$Author))
authorshipmessage <- c(packageDescription("JAMS")$Title, paste("JAMS version", packageVersion("JAMS")), authors, paste("Contact:", "john.mcculloch@nih.gov"), "National Cancer Institute", "National Institutes of Health", "Bethesda, MD, USA")
if (missing(outputdir)){
outputdir <- opt$outdir
}
pdffn <- file.path(outputdir, paste(paste(opt$prefix, "JAMSalpha", "report", sep="_"), "pdf", sep = "."))
pdf(pdffn, paper = "a4r")
#Header with run info
plot.new()
grid.table(opt$projinfo, rows = NULL, cols = NULL, theme = ttheme_default(base_size = 12))
plot.new()
grid.table(authorshipmessage, rows = NULL, cols = NULL, theme = ttheme_default(base_size = 15))
#Print read and assembly plots, if available
if ("readplots" %in% elementsinopt){
print(opt$readplots)
}
#Print taxonomic analysis, if applicable
if (("taxonomy" %in% elements)){
taxlvlspresent <- colnames(opt$LKTdose)[colnames(opt$LKTdose) %in% c("Domain", "Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species", "LKT")]
if (report_style %in% c("metagenome", "metatranscriptome")){
for (taxl in taxlvlspresent){
print(pareto_abundance(LKTdose = opt$LKTdose, samplename = samplename, taxlevel = taxl, assemblystats = opt$assemblystats, showcompletenesstext = TRUE, showcompletenessline = FALSE, showcontigline = FALSE, showparetocurve = TRUE, showpercentagetext = TRUE))
}
#Print an extra pareto chart for LKT with genome completeness
print(pareto_abundance(LKTdose = opt$LKTdose, assemblystats = opt$assemblystats, samplename = samplename, taxlevel = "LKT", showcompletenesstext = TRUE, showcompletenessline = TRUE, showcontigline = FALSE, showparetocurve = TRUE, showpercentagetext = TRUE))
} else if (report_style == "isolate"){
donut_list <- lapply(c("LKT", "Species", "IS1"), function (x) { plot_genome_taxonomy_donut(opt = opt, taxlevel = x, percentage_threshold = 99, ntop = 10, brewer_palette = 3) } )
names(donut_list) <- c("Species", "IS1", "LKT")
Isolateassemblystatstable <- ggtexttable(subset(opt$assemblystats, TaxLevel == "Isolate")[ , 2:ncol(opt$assemblystats)], rows = NULL, theme = ttheme("mBlue", base_size = 8))
donuts_top <- ggarrange(plotlist = donut_list, ncol = 3, nrow = 1, align = "h")
donuts <- ggarrange(donuts_top, Isolateassemblystatstable, ncol = 1, nrow = 2, align = "hv")
print(donuts)
}
if ("ucoplots" %in% names(opt)){
print(opt$ucoplots)
}
}
#If transcriptome, report proportion of CDS to rRNA
#if(opt$analysis %in% c("metatranscriptome", "isolaternaseq")){
RNArelabund <- subset(opt$featuredose, Analysis == "FeatType")
RNArelabund <- RNArelabund[, c("Accession", "Description", "NumBases")]
totbases <- sum(RNArelabund$NumBases)
RNArelabund$PPM <- round((RNArelabund$NumBases/totbases)*1000000, 0)
dfRNA <- RNArelabund[,c("Accession", "PPM")]
colnames(dfRNA) <- c("FeatType", "PPM")
dfRNA$Pct <- round(dfRNA$PPM/10000, 1)
dfRNA$lbls <- paste0(dfRNA$FeatType, " = ", dfRNA$Pct, "%")
dfRNA$lbls <- as.factor(dfRNA$lbls)
dfRNA <- subset(dfRNA, Pct > 0)
Featcolours <- c("green", "red", "blue", "yellow", "brown", "black")[1:nrow(dfRNA)]
rnaP <- ggplot(dfRNA, aes(x = "", y = Pct, fill = lbls)) + geom_bar(width = 1, stat = "identity") + scale_fill_manual(values = Featcolours) + coord_polar("y") + labs(title = paste(opt$prefix, "% Feature type found in Contigs", sep=" - "))
rnaP <- rnaP + labs(fill = "Feature Type")
rnaP <- rnaP + theme( axis.title.x = element_blank(), axis.title.y = element_blank(), panel.border = element_blank(), panel.grid=element_blank(), axis.ticks <- element_blank(), plot.title=element_text(size=14, face="bold"))
print(rnaP)
#}
#Print consolidated 16S table, if available
if (("SixteenSid" %in% colnames(opt$featuredata)) && ("SixteenSid" %in% elements)){
taxa_16S_cons <- subset(opt$featuredata, SixteenSid != "none")
if (nrow(taxa_16S_cons) > 0){
taxa_16S_cons <- taxa_16S_cons[, c("Feature", "LengthDNA", "Contig", "SixteenSid")]
taxa_16S_cons <- left_join(taxa_16S_cons, opt$contigsdata)
taxa_16S_cons <- taxa_16S_cons[, c("Feature", "LengthDNA", "Contig", "Length", "SixteenSid", "LKT")]
colnames(taxa_16S_cons) <- c("Gene", "Length_16S", "Contig", "Contig_Length", "LKT_dada2", "LKT_kraken")
ti <- "Consolidated 16S rRNA table"
print_table(tb=taxa_16S_cons, tabletitle=ti, fontsize=5, numrows=20)
} else {
print("No 16S rRNA features found in contigs")
}
}
availblastanalyses <- elementsinopt[!(elementsinopt %in% c("readplots", "SixteenSid"))]
batp <- NULL
for (ba in availblastanalyses){
batp <- opt[[ba]]
if (ba == "resfinder"){
contig2length <- opt$contigsdata[, c("Contig","Length")]
batp <- batp[ , c("Accession","Class","Pident","Contig","LKT")]
batp <- left_join(batp, contig2length, by = "Contig")
colnames(batp) <- c("DB_Hit","Function","Percent_ID","Contig","LKT", "Length")
batp <- batp[,c("DB_Hit","Function","Percent_ID","Contig","Length","LKT")]
} else if (ba == "abricate"){
contig2length <- opt$contigsdata[, c("Contig","Length")]
batp <- batp[ , c("Database", "Accession", "Phenotype", "Pident", "Contig", "LKT")]
batp <- left_join(batp, contig2length, by = "Contig")
colnames(batp) <- c("Database", "DB_Hit", "Phenotype", "Percent_ID", "Contig", "LKT", "Length")
batp <- batp[ , c("Database", "DB_Hit", "Phenotype", "Percent_ID", "Contig", "Length", "LKT")]
} else {
batp <- batp[ , c("Accession","Product","Pident","LKT")]
colnames(batp) <- c("DB_Hit","Function","Percent_ID","LKT")
}
#Shorten hit because it is too verbose from resfinder and vfdb
batp$DB_Hit <- gsub(":.*", "", batp$DB_Hit)
ti <- switch(ba, "resfinder" = "Known antibiotic resistance genes", "vfdb" = "Virulence factors present in VFDB \n www.ncbi.nlm.nih.gov/pubmed/15608208", "plasmidfinder" = "Plasmid Replicon-associated genes present in PlasmidFinder \n https://www.ncbi.nlm.nih.gov/pubmed/24777092")
print_table(tb = batp, tabletitle = ti, fontsize = 6, numrows = 20)
}
if("func" %in% elements){
funcanalyses <- unique(opt$featuredose$Analysis)[!(unique(opt$featuredose$Analysis) %in% c(availblastanalyses, "napdos", "serofinderO", "serofinderH"))]
for(fa in funcanalyses){
plot_wordcloud_sample(opt=opt, analysis=fa, removeunclassifieds=TRUE)
}
}
plot.new()
grid.table(authorshipmessage, rows = NULL, cols = NULL, theme = ttheme_default(base_size = 15))
dev.off()
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.