54.《Bioinformatics Data Skills》之

今天学习GenomicRanges包使用最后一小节，GRanges对象覆盖度计算。

测序简单模拟

测序深度计算公式：

其中C为覆盖度（或深度），N为read数目，L为read长度，G为测序区域长度。

假设read长度为150bp，以小鼠的最短的染色体19号染色体对象，其长度为61431566bp（如下）：

> library(TxDb.Mmusculus.UCSC.mm10.ensGene)
> txdb <- TxDb.Mmusculus.UCSC.mm10.ensGene
> chr19_len <- seqlengths(txdb)["chr19"]
> chr19_len
   chr19
61431566

如果产生5X的覆盖度，那么需要61431566*5/150 = 2047719条read：

> set.seed(0)
> read_len <- 150
> start_pos <- sample(1:(chr19_len-read_len), 2047719, replace = T)
> reads <- GRanges("chr19", IRanges(start_pos, width = read_len))

覆盖度计算

read覆盖度为：

> reads_cov <- coverage(reads)
RleList of length 1
$chr19
integer-Rle of length 61431565 with 3897731 runs
  Lengths: 12  2 23  8 47 13 57  2 23 14 41 ...  7 39 29 36  5 19 15  7 68 36
  Values :  0  1  2  3  4  5  6  5  4  3  4 ...  4  5  4  5  6  5  4  3  2  1

平均深度为5X：

> mean(reads_cov)
chr19
    5

有两种方式确定未覆盖read区域长度，其一是通过table与==函数：

> table(reads_cov == 0)
         FALSE     TRUE
chr19 61016815   414750

其二是使用helper函数，速度更快：

> sum(runLength(reads_cov)[runValue(reads_cov)==0])
 chr19
414750

可知没有read覆盖区域占测序区域比例为：

> 414750/chr19_len
      chr19
0.006751415

约占0.6%，比较有意思的是，未覆盖区域的比例是服从泊松分布的（，即)（如下），这样的现象在随机分布read的鸟枪测序法中常见。

> exp(-5)
[1] 0.006737947