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

knitr::opts_chunk$set(echo = TRUE)
require(knitr)
require(kableExtra)
require(plotly)

output <- readRDS("output.rds")
se <- readRDS("deResult.rds")
param <- readRDS("param.rds")


debug <- FALSE

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

FastqScreen_Result {.tabset}

fastqData=fastqData_ppData
fastqDataAdapters=fastqData_rawData

FastqScreen Mapping Rates to Adapters and rRNA

m <- list(l = 80, r = 80, b = 200, t = 100, pad = 0) # more margin space
# Mapping rate
#p <- plot_ly(x=names(fastqData$MappingRate),
#             y=fastqData$MappingRate,
#             type="bar") %>%
#  layout(title = "Overall MappingRate",
#         yaxis = list(title = "MappedReads in %", range = c(0,100)),
#         margin = m)
#p

# MappingRateAdapters
p <- plot_ly(x=names(fastqDataAdapters$MappingRate),
             y=fastqDataAdapters$MappingRate,
             type="bar") %>%
  layout(title = "MappingRate to Adapters - Bowtie2, local Alignment",
         yaxis = list(title = "MappedReads in %", range = c(0,100)),
         margin = m)
p

# Reads
#p <- plot_ly(x=names(fastqData$Reads), y=fastqData$Reads/1e3, type="bar") %>%
#  layout(title = "ProcessedReads",
#         yaxis = list(title = "#Reads in K"),
#         margin = m)
#p

# rRNA Mapping
rRNA_mappingRate = as.data.frame(rRNA_strandInfo/(param[['nReads']]/100))
rRNA_mappingRate$sample <- rownames(rRNA_mappingRate)
p <- plot_ly(rRNA_mappingRate, x=~sample, y=~Sense, 
             type="bar", name="sense") %>%
    add_trace(y = ~Antisense, name = "antisense") %>%
  layout(title = "rRNA Silva Mapping - Bowtie2, end-to-end Alignment",
         yaxis = list(title = "rRNA-Mapping-Rate in %"),
         xaxis = list(title = ""),
         barmode = 'stack', margin = m)
p

# MappingRate
par(mar=c(10.1, 4.1, 4.1, 2.1))
    bplt = barplot(fastqData$MappingRate, las=2, ylim=c(0,100), ylab="MappedReads in %", main="Overall MappingRate", col="royalblue3",
                   names.arg=rep('',length(ezSplitLongLabels(names(fastqData$MappingRate)))))
    if(min(fastqData$MappingRate) < 8){
      #text(y=fastqData$MappingRate+2, font=2, x=bplt, labels=as.character(fastqData$MappingRate), cex= 1, xpd=TRUE)
      text(y=fastqData$MappingRate+2, font=2, x=bplt, srt = 90, adj = 0, labels=as.character(fastqData$MappingRate), cex= 1, xpd=TRUE)
    } else {
      # text(y=fastqData$MappingRate-5, font=2, x=bplt, 
      #      labels=as.character(fastqData$MappingRate), cex= 1.1, col='white', 
      #      xpd=TRUE)
      text(y=fastqData$MappingRate-5, font=2, x=bplt, srt = 90, adj = 1,
           labels=as.character(fastqData$MappingRate), cex= 1.1, col='white', 
           xpd=TRUE)
    }
    text(x = bplt, y = par("usr")[3] - 2, srt = 45, adj = 1, 
         labels = ezSplitLongLabels(names(fastqData$MappingRate)), xpd = TRUE)

# MappingRateAdapters
par(mar=c(10.1, 4.1, 4.1, 2.1))
    bplt = barplot(fastqDataAdapters$MappingRate, las=2, ylim=c(0,100), ylab="MappedReads in %", main="MappingRate to Adapters without trimming", col="royalblue3",
                   names.arg=rep('',length(ezSplitLongLabels(names(fastqDataAdapters$MappingRate)))))
    if(min(fastqDataAdapters$MappingRate) < 8){
      text(y=fastqDataAdapters$MappingRate+2, font=2, x=bplt, srt = 90, adj = 0, labels=as.character(fastqDataAdapters$MappingRate), cex= 1, xpd=TRUE)
    } else {
      text(y=fastqDataAdapters$MappingRate-5, font=2, x=bplt, srt = 90, adj = 1, labels=as.character(fastqDataAdapters$MappingRate), cex= 1.1, col='white', xpd=TRUE)
    }
    text(x = bplt, y = par("usr")[3] - 2, srt = 45, adj = 1, labels = ezSplitLongLabels(names(fastqDataAdapters$MappingRate)), xpd = TRUE)

# Reads
par(mar=c(10.1, 4.1, 4.1, 2.1))
    bplt = barplot(fastqData$Reads/1000, las=2, ylab="#Reads in K", main="ProcessedReads", col="lightblue",
            names.arg=rep('',length(ezSplitLongLabels(names(fastqData$MappingRate)))))
    text(x = bplt, y = par("usr")[3] - 2, srt = 45, adj = 1, labels = ezSplitLongLabels(names(fastqData$MappingRate)), xpd = TRUE)

# rRNA Mapping
par(mar=c(10.1, 4.1, 4.1, 2.1))
    rRNA_mappingRate = t(rRNA_strandInfo/(param[['nReads']]/100))
    bplt = barplot(rRNA_mappingRate, las=2, ylim=c(0,min(max(colSums(rRNA_mappingRate))+20,100)), ylab="rRNA-Mapping-Rate in %", main="rRNA Silva Mapping", col=c("lightblue","darkblue"),
           names.arg=rep('',length(ezSplitLongLabels(rownames(rRNA_strandInfo)))), legend.text = T)
    text(x = bplt, y = par("usr")[3] - min(max(colSums(rRNA_mappingRate))+20,100) * 0.02, srt = 45, adj = 1, labels = ezSplitLongLabels(rownames(rRNA_strandInfo)), xpd = TRUE)

FastqScreen Mapping Per Sample

for (nm in rownames(dataset)){
  par(mar=c(10.1, 4.1, 4.1, 2.1))
  x = fastqData$CommonResults[[nm]]
  if (nrow(x) > 0){
    bplt = barplot(t(x), las=2, ylim=c(0,100), 
                   legend.text=T, ylab="Mapped Reads in %", main=nm, names.arg=rep('', nrow(x)))
      text(x = bplt, y = par("usr")[3] - 2, srt = 45, adj = 1, labels = rownames(x), xpd = TRUE)
  } else {
    plot(1,1, type="n", axes=FALSE, main=nm, xlab="", ylab="", frame=TRUE)
    text(1,1, "no hits found")
  }
}

Mapping to RefSeq mRNA Per Sample

for (nm in rownames(dataset)){
      par(mar=c(10.1, 4.1, 4.1, 2.1))
      x = speciesPercentageTop[[nm]]
      if (is.null(x)) x = matrix(0, 2, 1, dimnames=list(c('UniqueSpeciesHits','MultipleSpeciesHits'),'Misc'))
      bplot = barplot(t(x), col=c("royalblue3", "lightblue"), las=2, ylim=c(0,100),
                      legend.text=T, ylab="Mapped Reads in %", main=nm, names.arg=rep('',nrow(x)) )
      text(y=t(x)[ 1,] + 5, x=bplot, font = 2, labels=t(x)[ 1, ], cex=1.1, col='black')
      text(x = bplot, y = par("usr")[3] - 2, srt = 45, adj = 1, 
           labels = rownames(x), xpd = TRUE)
}

Mapping to Genomes with Kraken

Abundances above 5% are truncated.

for (i in c(1:length(krakenResult))){
      par(mar=c(15.1, 4.1, 4.1, 2.1))
        bplot = barplot(krakenResult[[i]]$readPercentage, names.arg = krakenResult[[i]]$name, col = 'royalblue3', 
                main = names(krakenResult)[i], ylab = 'Mapped Reads in %', las = 2, ylim=c(0, 5))
}

Virus Check

if(param[['virusCheck']]){
  for (nm in rownames(dataset)){ 
    par(mar=c(18.1, 7.1, 2.1, 2.1))
        x = speciesPercentageTopVirus[[nm]]
        if (is.null(x)) x = matrix(0, 2, 1, dimnames=list(c('UniqueSpeciesHits','MultipleSpeciesHits'),'Misc'))
        bplot = barplot(t(x), col=c("royalblue3", "lightblue"), las = 2, ylim = c(0,100),
                        legend.text=T, ylab="Mapped Reads in %", main=nm, names.arg=rep('',nrow(x)) )
        text(y=t(x)[ 1,] + 5, x=bplot, font = 2, labels=t(x)[ 1, ], cex = 1.1, col = 'black')
        text(x = bplot, y = par("usr")[3] - 2, srt = 60, adj = 1, 
             labels = rownames(x), xpd = TRUE)
  }
}

RIN

RIN vs Read Count

if (nrow(dataset) > 1){
  rinCol <- colnames(dataset)[grep('RIN', colnames(dataset))]
  if(length(rinCol) == 1){
    a <- makeScatterplot(dataset, colname1 = rinCol, colname2='Read Count')
    a
    }
}

RIN vs LibConc

if (nrow(dataset) > 1){
  rinCol <- colnames(dataset)[grep('RIN', colnames(dataset))]
  if(length(rinCol) == 1){
    libQuant <- colnames(dataset)[grep('LibConc_100_800bp', colnames(dataset))]  
    if(length(libQuant) == 1){
      b <- makeScatterplot(dataset, colname1 = rinCol, colname2=libQuant)
      b
    }
  }
}

RIN vs rRNA content

if (nrow(dataset) > 1){
  rinCol <- colnames(dataset)[grep('RIN', colnames(dataset))]
  if(length(rinCol) == 1){
    dataset[['rRNA_Content']] <- rRNA_mappingRate$Sense + rRNA_mappingRate$Antisense
    c <- makeScatterplot(dataset, colname1 = rinCol, colname2='rRNA_Content')
    c
  }
}

Settings

getAppVer <- function(appName) { sub("^.+/([^/]+)$", "\\1", Sys.getenv(appName)) }
settings = character()
settings["Configuration File:"] = param$confFile
settings["RefSeq mRNA Reference:"] = REFSEQ_mRNA_REF
settings["FastqScreen Version:"] = getAppVer("FastQScreen")
settings["Bowtie2 Version:"] = getAppVer("Bowtie2")
settings["Bowtie2 Parameters:"] = param$cmdOptions
settings["Minimum AlignmentScore:"] = param$minAlignmentScore
settings["TopSpecies:"] = param$nTopSpecies
kable(as.data.frame(settings), col.names=NA, row.names=TRUE, format="html") %>%
  kable_styling(bootstrap_options = "striped", full_width = F,
                position = "left")

Input Dataset

ezInteractiveTableRmd(values=dataset)

SessionInfo

ezSessionInfo()

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