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

```r

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

```r
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."
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 samlpe
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))

### create plots: 3. richness. Here we use only the unfiltered dataset.
plotTitle <- "Estimate of the community's richness the samples."
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/saturation plots.
abundTableRar <- t(otu_table(physeqFullObject)@.Data)
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)

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)
plot(rarefPlot_1)
#plot(rarefPlot_2)

Heatmaps

show_pHeatmap()

SessionInfo

ezSessionInfo()

uzh/ezRun documentation built on April 24, 2024, 4:01 p.m.