In uzh/ezRun: An R meta-package for the analysis of Next Generation Sequencing Data

knitr::opts_chunk$set(echo = TRUE,warning=FALSE, message=FALSE, out.width = "49%")
## This report requires the phyloseq  R object generated by the 16S App
debug <- FALSE

QCChimeraObject <- readRDS("QCChimeraObject.rds")
physeqFullObject <- readRDS("physeqFullObject.rds")
physeqObjectNoTreeUnfilt <- readRDS("physeqObjectNoTreeUnfilt.rds")
param = readRDS("param.rds")
numTopRanks <- param$numTopRanks 
rawCount <- param$rawCount
rank <- param$taxonomicRank
sampleFraction <-param$sampleFraction
isGroupThere = param$grouping != ""
if (isGroupThere) {
        group <- param$grouping
        if (param$sampleGroup == "" | param$refGroup == ""){
            stop("Both sample and reference groups must be specified")
        } else{
            sampleGroup <- param$sampleGroup
            refGroup <- param$refGroup
        }
}
if (isGroupThere){
    areThereMultVar <- TRUE
    } else {
    areThereMultVar <- FALSE
}

filtStepPlot <- ggplot(QCChimeraObject$filtStep,aes(x=fStep,y=Freq, group=sample)) + geom_line(aes(color=sample)) + geom_point(aes(color=sample)) 
filtStepPlot <- filtStepPlot + xlab("Filtering Step") + ylab("Read Count") +
    theme(axis.text.x = element_text(angle = 270))

chimPlot <- chimeraSummaryPlot(QCChimeraObject$chimera)

### create plots: 1. abundance
abundPlotSampleList <- list()
abundPlotGroupList <- list()
abundPlotSampleListUnfilt <- list()
abundPlotGroupListUnfilt <- list()
isAbundPlotSample <- list()
isAbundPlotGroup <-list()
isAbundPlotSampleUnfilt <- list()
isAbundPlotGroupUnfilt <- list()
taxRanksVec <- attr(physeqFullObject@tax_table@.Data,"dimnames")[[2]]
taxRanks <- setNames(as.list(taxRanksVec),taxRanksVec)
for (rankNames  in taxRanks){
    tempAbundPlotFilt <- abundPlot(rankNames,physeqFullObject,xAesLogic="S",numTopRanks,group)
    if (tempAbundPlotFilt$stop==TRUE) {
        isAbundPlotSample[[rankNames]] <- FALSE
        next
    } else {
        isAbundPlotSample[[rankNames]] <- TRUE
        percentFrac <- sampleFraction*100
        plotTitle <- paste("Community composition in each sample\nat the rank", rankNames, "for the filtered dataset. Only taxa observed at least",rawCount, "\ntimes in at least", percentFrac,"% of the samples are kept.")
        abundPlotSampleList[[rankNames]] <- tempAbundPlotFilt$abPlot +labs(title=plotTitle) + 
            theme(plot.title=element_text(size=11,hjust=0.5))
    }
}
for (rankNames  in taxRanks){
    tempAbundPlotUnfilt <-  abundPlot(rankNames,physeqObjectNoTreeUnfilt,xAesLogic="S",numTopRanks,group)
    if (tempAbundPlotUnfilt$stop==TRUE) {
        isAbundPlotSampleUnfilt[[rankNames]] <- FALSE
        next
    } else {
        isAbundPlotSampleUnfilt[[rankNames]] <- TRUE
        plotTitle <- paste("Community composition at the rank", rankNames, "\nfor the unfiltered dataset.")
        abundPlotSampleListUnfilt[[rankNames]] <- tempAbundPlotUnfilt$abPlot+labs(title=plotTitle) + 
            theme(plot.title=element_text(size=11,hjust=0.5))
    }
}
if (isGroupThere) {
    for (rankNames  in taxRanks){
        groupAes <- colnames(sample_data(physeqObjectNoTreeUnfilt))[1]
        tempAbundPlotFilt <- abundPlot(rankNames,physeqFullObject,xAesLogic="G",numTopRanks,group)

        if (tempAbundPlotFilt$stop==TRUE) {
            isAbundPlotGroup[[rankNames]] <- FALSE
            next
        } else {
            isAbundPlotGroup[[rankNames]] <- TRUE
            plotTitle <- paste("Community composition in each group\nat the rank", rankNames, "for the filtered dataset. Only taxa observed at least",rawCount, "\ntimes in at least", percentFrac,"% of the samples are kept.")
            abundPlotGroupList[[rankNames]] <- tempAbundPlotFilt$abPlot +labs(title=plotTitle) + 
                theme(plot.title=element_text(size=11,hjust=0.5))
        }
    }

    for (rankNames  in taxRanks){
        tempAbundPlotUnfilt <-  abundPlot(rankNames,physeqObjectNoTreeUnfilt,xAesLogic="G",numTopRanks,group)
        if (tempAbundPlotUnfilt$stop==TRUE) {
            isAbundPlotGroupUnfilt[[rankNames]] <- FALSE 
            next
        } else {
            isAbundPlotGroupUnfilt[[rankNames]] <- TRUE
            plotTitle <- paste("Community composition in each group at the rank", rankNames, "\nfor the unfiltered dataset.")
            abundPlotGroupListUnfilt[[rankNames]] <- tempAbundPlotUnfilt$abPlot +labs(title=plotTitle) + 
                theme(plot.title=element_text(size=11,hjust=0.5))
        }
    }
}

top20taxaWithPerc <- communityPercSummTable(physeqFullObject,rank)

### create plots: 2. ordination by taxa
plotTitle <- "Clustering of the taxa using using Bray-Curtis."
ordPlotTaxa <- ordPlot(rank,physeqFullObject,"taxa",areThereMultVar,numTopRanks,isGroupThere)+
    labs(title=plotTitle) + 
    theme(plot.title=element_text(size=11,hjust=0.5))

### create plots: 2. ordination by sample
plotTitle <- "Clustering of the samples."
ordPlotSample<- ordPlot(rank,physeqFullObject,"samples",areThereMultVar,numTopRanks,isGroupThere)+
    labs(title=plotTitle) + 
    theme(plot.title=element_text(size=11,hjust=0.5), legend.key.size = unit(0.4, "cm"))

### create plots: 3. richness. Here we use only the unfiltered dataset.
plotTitle <- "Estimate of the community's richness the samples (unfiltered dataset)."
plotRichnessBySample <- plot_richness(physeqObjectNoTreeUnfilt, measures=c("Shannon")) +labs(title=plotTitle) + 
    theme(plot.title=element_text(size=11,hjust=0.5))
if (isGroupThere) {
    plotTitle <- "Boxplots of the community's richness in the groups."
    plotRichnessByGroup <- groupModRichPlot(physeqObjectNoTreeUnfilt, x=group,measures=c("Shannon"))+labs(title=plotTitle) + 
        theme(plot.title=element_text(size=11,hjust=0.5))
}

### rarefaction: 3. rarefaction plots.
abundTableRar <- data.frame(t(otu_table(physeqFullObject)))
iNEXT_data <- iNEXT(abundTableRar, q = 0, datatype = "abundance")
#rarefactionPlot2 <- ggiNEXT(iNEXT_data, type=1)
#https://cran.r-project.org/web/packages/iNEXT/vignettes/Introduction.html
# plotTitle_1 <- paste("Rarefaction plots based on the\nselection of the most abundant", numTopRanks, "ranks")
# plotTitle_2 <- paste("Saturation plots based on the\nselection of the most abundant", numTopRanks, "ranks")
# rarefPlot_1 <- rarefactionPlot(abundTableRar,1) +labs(title=plotTitle_1) + 
#     theme(plot.title=element_text(size=11,hjust=0.5))
# rarefPlot_2 <- rarefactionPlot(abundTableRar,2) +labs(title=plotTitle_2) + 
#     theme(plot.title=element_text(size=11,hjust=0.5))

### create plots:5. pheatmap
show_pHeatmap <- heatmapForPhylotseqPlotPheatmap(physeqFullObject,areThereMultVar,
                                                 isGroupThere,rank)

### run deseq part as we need the output
if (isGroupThere){
    deseqResults <- phyloSeqToDeseq2_tableAndPlots(physeqFullObject,rank,group=group,
                                                   sampleGroup,refGroup)
    diffAbundTableName <- "diffAbundanceTable.txt"
}

Started on r format(Sys.time(), "%Y-%m-%d %H:%M:%S")

Phyloseq report {.tabset}

QC and chimera report

par(mfrow=c(1,2), las=1)
plot(filtStepPlot)
plot(chimPlot)

Abundance barplots

par(mfrow=c(2,2), las=1)
if (isAbundPlotSample[[rank]]) {
    plot(abundPlotSampleList[[rank]])
}else{
    message <- paste("After filtering, no taxa remaining wih a", rank,"annotaton. Nothing to show.")
    format(message)
}
if (isAbundPlotSampleUnfilt[[rank]]) {
    plot(abundPlotSampleListUnfilt[[rank]])
}else{
    message <- paste("No taxa annotated at the", rank,"level. Nothing to show.")
    format(message)
}
if (isGroupThere) {
    if (isAbundPlotGroup[[rank]]) {
        plot(abundPlotGroupList[[rank]])
    }else{
        message <- paste("After filtering, no taxa remaining wih a", rank,"annotaton. Nothing to show.")
        format(message)  
    }
    if (isAbundPlotGroupUnfilt[[rank]]) {
        plot(abundPlotGroupListUnfilt[[rank]])
    }else{
        message <- paste("No taxa annotated at the", rank,"level. Nothing to show.")
        format(message)
    }
}

Percentage of most abundant taxa

DT::datatable(top20taxaWithPerc)

Clustering (ordination) plots

par(mfrow=c(1,2), las=1)
plot(ordPlotTaxa)
if (isGroupThere) {
    plot(ordPlotSample)
}

Richness plots

par(mfrow=c(1,2), las=1)
plot(plotRichnessBySample)
if (isGroupThere) {
    plot(plotRichnessByGroup)
}

Rarefaction plots

#par(mfrow=c(1,2), las=1)
rarecurve(otu_table(physeqFullObject)@.Data, step=50, cex=0.5, ylab = "OTUs")
#plot(rarefactionPlot2)

Heatmaps

show_pHeatmap()

Two-goups comparison

par(mfrow=c(1,3), las=1)
plot(deseqResults$vPlot)
plot(deseqResults$logPlot)
if (deseqResults$isAllNa) {
    format(deseqResults$isAllNaMsg)
} else {
    plot(deseqResults$pieChart)
}

Full table with all the results to download

write.table(deseqResults$fullTable,diffAbundTableName,row.names = F, col.names = T, quote = F,sep = "\t")
zipped = zipFile(diffAbundTableName)

r zipped

Summary table (top 20 differentially abundant OTUs)

DT::datatable(deseqResults$tableToReport)

SessionInfo

ezSessionInfo()

uzh/ezRun documentation built on April 14, 2024, 5:09 a.m.