Last data update: 2014.03.03

R: Get depth of coverage from whole exome sequencing
getcoverageR Documentation

Get depth of coverage from whole exome sequencing

Description

Gets depth of coverage for each exon across all samples from whole exome sequencing files.

Usage

getcoverage(bambedObj, mapqthres)

Arguments

bambedObj

Object returned from getbambed

mapqthres

Mapping quality threshold hold of reads.

Value

Y

Read depth matrix

readlength

Vector of read length for each sample

Author(s)

Yuchao Jiang yuchaoj@wharton.upenn.edu

See Also

getbambed

Examples

library(WES.1KG.WUGSC)
dirPath <- system.file("extdata", package = "WES.1KG.WUGSC")
bamFile <- list.files(dirPath, pattern = '*.bam$')
bamdir <- file.path(dirPath, bamFile)
sampnameFile <- file.path(dirPath, "sampname")
sampname <- as.matrix(read.table(sampnameFile))
chr <- 22
bambedObj <- getbambed(bamdir = bamdir, bedFile = file.path(dirPath, 
    "chr22_400_to_500.bed"), sampname = sampname,
    projectname = "CODEX_demo", chr)
bamdir <- bambedObj$bamdir
sampname <- bambedObj$sampname
ref <- bambedObj$ref
projectname <- bambedObj$projectname
chr <- bambedObj$chr
coverageObj <- getcoverage(bambedObj, mapqthres = 20)
Y <- coverageObj$Y
readlength <- coverageObj$readlength

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(CODEX)
Loading required package: Rsamtools
Loading required package: GenomeInfoDb
Loading required package: stats4
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, cbind, colnames, do.call, duplicated, eval, evalq,
    get, grep, grepl, intersect, is.unsorted, lapply, lengths, mapply,
    match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, rank,
    rbind, rownames, sapply, setdiff, sort, table, tapply, union,
    unique, unsplit

Loading required package: S4Vectors

Attaching package: 'S4Vectors'

The following objects are masked from 'package:base':

    colMeans, colSums, expand.grid, rowMeans, rowSums

Loading required package: IRanges
Loading required package: GenomicRanges
Loading required package: Biostrings
Loading required package: XVector
Loading required package: BSgenome.Hsapiens.UCSC.hg19
Loading required package: BSgenome
Loading required package: rtracklayer

Attaching package: 'CODEX'

The following object is masked from 'package:BiocGenerics':

    normalize

> png(filename="/home/ddbj/snapshot/RGM3/R_BC/result/CODEX/getcoverage.Rd_%03d_medium.png", width=480, height=480)
> ### Name: getcoverage
> ### Title: Get depth of coverage from whole exome sequencing
> ### Aliases: getcoverage
> ### Keywords: package
> 
> ### ** Examples
> 
> library(WES.1KG.WUGSC)
> dirPath <- system.file("extdata", package = "WES.1KG.WUGSC")
> bamFile <- list.files(dirPath, pattern = '*.bam$')
> bamdir <- file.path(dirPath, bamFile)
> sampnameFile <- file.path(dirPath, "sampname")
> sampname <- as.matrix(read.table(sampnameFile))
> chr <- 22
> bambedObj <- getbambed(bamdir = bamdir, bedFile = file.path(dirPath, 
+     "chr22_400_to_500.bed"), sampname = sampname,
+     projectname = "CODEX_demo", chr)
> bamdir <- bambedObj$bamdir
> sampname <- bambedObj$sampname
> ref <- bambedObj$ref
> projectname <- bambedObj$projectname
> chr <- bambedObj$chr
> coverageObj <- getcoverage(bambedObj, mapqthres = 20)
Getting coverage for sample NA06994: read length 100.
Getting coverage for sample NA10847: read length 100.
Getting coverage for sample NA11840: read length 100.
Getting coverage for sample NA12249: read length 100.
Getting coverage for sample NA12716: read length 100.
Getting coverage for sample NA12750: read length 100.
Getting coverage for sample NA12751: read length 100.
Getting coverage for sample NA12760: read length 100.
Getting coverage for sample NA12761: read length 100.
Getting coverage for sample NA12763: read length 100.
Getting coverage for sample NA18966: read length 100.
Getting coverage for sample NA18967: read length 100.
Getting coverage for sample NA18968: read length 100.
Getting coverage for sample NA18969: read length 100.
Getting coverage for sample NA18970: read length 100.
Getting coverage for sample NA18971: read length 100.
Getting coverage for sample NA18972: read length 100.
Getting coverage for sample NA18973: read length 100.
Getting coverage for sample NA18974: read length 100.
Getting coverage for sample NA18975: read length 100.
Getting coverage for sample NA18976: read length 100.
Getting coverage for sample NA18981: read length 100.
Getting coverage for sample NA18987: read length 100.
Getting coverage for sample NA18990: read length 100.
Getting coverage for sample NA18991: read length 100.
Getting coverage for sample NA19098: read length 100.
Getting coverage for sample NA19119: read length 100.
Getting coverage for sample NA19131: read length 100.
Getting coverage for sample NA19137: read length 100.
Getting coverage for sample NA19138: read length 100.
Getting coverage for sample NA19141: read length 100.
Getting coverage for sample NA19143: read length 100.
Getting coverage for sample NA19144: read length 100.
Getting coverage for sample NA19152: read length 100.
Getting coverage for sample NA19153: read length 100.
Getting coverage for sample NA19159: read length 100.
Getting coverage for sample NA19160: read length 100.
Getting coverage for sample NA19171: read length 100.
Getting coverage for sample NA19200: read length 100.
Getting coverage for sample NA19201: read length 100.
Getting coverage for sample NA19204: read length 100.
Getting coverage for sample NA19206: read length 100.
Getting coverage for sample NA19207: read length 100.
Getting coverage for sample NA19209: read length 100.
Getting coverage for sample NA19210: read length 100.
Getting coverage for sample NA19223: read length 100.
> Y <- coverageObj$Y
> readlength <- coverageObj$readlength
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>