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

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

knitr::opts_chunk$set(echo = TRUE)
## This report requires the object resultList, dataset, param
require(knitr)
require(kableExtra)

## debug:
#require(ezRun)
#load("/home/gtan/analysis/gtan/p2438-RNABamStats/resultList.RData")
#load("/home/gtan/analysis/gtan/p2438-RNABamStats/dataParam.rda")
## end of debug

conds = ezConditionsFromDataset(dataset, param=param)
samples = rownames(dataset)
sampleColors = getSampleColors(conds, samples)
bamFiles = dataset$BAM

if (param$writeIgvSessionLink){
#     titles[["Genome Browser"]] = "Genome Browser"
#     addTitle(doc, titles[[length(titles)]], 2, id=titles[[length(titles)]]) ## NOTEP: igv links are not yet added to the report.
  if (length(bamFiles) > 4){
    idx = which(!duplicated(conds))
    idx = idx[1:min(4, length(idx))]
  } else {
    idx = 1:length(bamFiles)
  }
  for (each in idx){
    writeIgvSession(genome=getIgvGenome(param), refBuild=param$ezRef["refBuild"],
                    file=basename(sub(".bam", "-igv.xml", bamFiles[each])),
                    bamUrls=paste(PROJECT_BASE_URL, bamFiles[each], sep="/"))
    writeIgvJnlp(jnlpFile=basename(sub(".bam", "-igv.jnlp", bamFiles[each])),
                 projectId=param$projectId,
                 sessionUrl=paste(PROJECT_BASE_URL, sub(".bam", "-igv.xml", 
                                                        bamFiles[each]), 
                                  sep="/"))
  }
}

DNA_QC_Statistics {.tabset}

Overview

is.numeric0 <- function(x) {
  identical(x, numeric(0))
}

basicStats = ezMatrix(0, rows=samples, cols=c('MappingRate', 'DuplicationRate', 'avgCoverage', 'errorRate', 'insertionRate', 'deletionRate'))
for (sm in samples){
    if(any(sapply(resultList[[sm]], is.numeric0))){
      resultList[[sm]][sapply(resultList[[sm]], is.numeric0)] = 0
    }
    basicStats[sm, 'MappingRate'] = resultList[[sm]]$mappingRate
    basicStats[sm, 'DuplicationRate'] = resultList[[sm]]$dupRate
    basicStats[sm, 'avgCoverage'] = resultList[[sm]]$avgCoverage
    basicStats[sm, 'errorRate'] = resultList[[sm]]$errorRate
    basicStats[sm, 'insertionRate'] = resultList[[sm]]$insertRate
    basicStats[sm, 'deletionRate'] = resultList[[sm]]$delRate
}

ezInteractiveTableRmd(basicStats, digits=3)

Multi-Matching Reported in Bam File

The plot holds for each sample the number of reads in Millions that have X matches in the target and are reported in the file.

mmValues = integer()
for (sm in samples){
  mmValues = union(mmValues, as.integer(names(resultList[[sm]]$multiMatchInFileTable)))
}
mmCounts = ezMatrix(0, rows=samples, cols=sort(mmValues))
for (sm in samples){
  mm = resultList[[sm]]$multiMatchInFileTable
  mmCounts[sm, names(mm)] = mm
}

alignmentCountBarPlot(mmCounts, relative=FALSE,
                      file="multiMatchInFile-barplot.txt")
alignmentCountBarPlot(mmCounts, relative=TRUE)

txtFile <- "read-alignment-statistics.txt"
ezWrite.table(mmCounts, file=txtFile, head="Sample")

r txtFile

InsertSize

plots <- list.files(".","fragmentSize_.*png$")
   for(i in 1:length(plots)){
      filename <- plots[i]
      cat("![](",filename,")")
   }

Library Complexity

plots <- list.files(".","libComplexity_.*png$")
   for(i in 1:length(plots)){
      filename <- plots[i]
      cat("![](",filename,")")
   }

Plots about DuplicationRate

Coverage

Qualimap

Multi-Sample Report

plots <- list.files(".","QualiMapStats_.*png$")
   for(i in 1:length(plots)){
      filename <- plots[i]
      cat("![](",filename,")")
   }

Per Sample Reports

require(DT)
qualimapResults = data.frame(SampleName = samples, ReportLink = file.path(samples, 'qualimapReport.html'), stringsAsFactors = F)
qualimapResults$ReportLink =  paste0('<a href="',qualimapResults$ReportLink,'">',qualimapResults$ReportLink,'</a>') 
datatable(qualimapResults, escape = F,
                                caption='', colnames=colnames(qualimapResults),
                                options=list(dom = 'Bfrtip', pageLength=25, autoWidth=TRUE))

Input Dataset

ezInteractiveTableRmd(dataset)

SessionInfo

ezSessionInfo()

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