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

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 <- set_names(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)

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

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

sessionInfo()

uzh/ezRun documentation built on May 4, 2024, 3:23 p.m.