前面我们在第8讲提到了公司给我的一个报告的统计表格，有人反映不会做。

本来应该只需要给我6亿条reads的（PE150测序，人30X），但是足足给了我8.9亿条！（但事实上很多paper发表的基因组高于60X的也不少）

表格里面提到了好几个概念，比如duplicate的reads，一般说的PCR造成的duplicate，在找变异的时候需要去除掉。然后是那些比对到了不同染色体的reads pair，虽然只有2.29% ，也是需要重点分析的。（前面我也讲了如何提取以及分别分析它们！）

如果只是想重新前面的这些统计指标，非常很简单，就是samtools工具提供了一个flagstat功能，用法及结果如下：

samtools flagstat P_jmzeng.final.bam

 

899361748 + 0 in total (QC-passed reads + QC-failed reads)

8597742 + 0 secondary

0 + 0 supplementary

132556557 + 0 duplicates

890858540 + 0 mapped (99.05% : N/A)

890764006 + 0 paired in sequencing

445382003 + 0 read1

445382003 + 0 read2

853255862 + 0 properly paired (95.79% : N/A)

881249604 + 0 with itself and mate mapped

1011194 + 0 singletons (0.11% : N/A)

20382234 + 0 with mate mapped to a different chr

12511988 + 0 with mate mapped to a different chr (mapQ>=5)

从结果很明显可以看到，公司的确给了我8.9亿的reads，这个没错，duplicate的情况也的确是132556557/899361748 = 14.73%，而mapping的情况和proper mapping的情况也都显示好了，可以看到公司用的是同样的命令和方法！没有什么神秘的，我们生信工程师，做得就是这个，而且可以做得更好。

后面的探索全基因组区域中碱基覆盖深度不低于多少X的比例，是需要画一个图，有非常多的现成的工具可以使用，包括 BedTools' genomeCov 、 GATK's DepthOfCoverage，还有Picard suite的几个命令。

其实本身原理很简单，就是把全基因组的每个坐标的depth都得到，然后得到depth的频数，然后画图。

我们可以对每条染色体单独来绘图，也可以针对全基因组来绘图，当然，公司给我们的数据是针对全基因组的。 脚本很简单：

samtools mpileup P_jmzeng.final.bam |perl -alne '{if($F[3]>100){$depth{"over100"}++}else{$depth{$F[3]}++}}END{print "$_\t$depth{$_}" foreach sort{$a <=> $b}keys %depth}' 当是需要运行很久，毕竟全基因组的bam文件太大了。

如果分开运行，可以对下面的各条染色体bam文件批量跑一个脚本

脚本如下：

ls P_jmzeng.final.REF*.bam |while read id

do

echo $id

samtools mpileup $id |perl -alne '{if($F[3]>100){$depth{"over100"}++}else{$depth{$F[3]}++}}END{print "$_\t$depth{$_}" foreach sort{$a <=> $b}keys %depth}' >$id.depth.txt

done

跑完之后，对每条染色体都会输出如下文件：

over100 110789

0 27065

1 1730286

2 2219409

3 2728526

4 3251046

5 3774335

6 4303971

~~~~~~~~~~~~~~~~~~~~~~~

后面省略94行，每一行都是两列，第一列是测序深度，第二列是有着该测序深度的位点是多少个！所以行的第二列加起来，就是染色体的长度！！！面这个例子是X染色体的，绘图如下，可以看到X染色体的测序深度其实并不怎么好，全基因组测序深度平均高达44X，可是这个X染色体超过44X的只占极少的比例。

我用excel表格简单画个图如下：（当然，作为一个高级生物信息学工程师，用excel表格是有点low，但是这里只是为了说明一个问题，我们后面还是写程序的，用R语言）

比如我们再看看10号染色体，我随意在R里面画了个图：

a=read.table('P_jmzeng.final.REF_10.bam.depth.txt',stringsAsFactors = F)

plot(a[,2],type = 'l',xaxt="n",lwd=3)

axis(1, at = 1:nrow(a),labels =a[,1] , las=1)

abline(v=44)

可以看到10号染色体的测序深度要显著好于X染色体，大部分的测序深度都超过了30X！！（不过，这里虽然用来R，但是出图很丑，虽然不是专业做可视化的，但是想调整美观一点问题也不大，就是好耗费时间。）

做出depth的累积曲线图了，也是很简单的！

a=read.table('P_jmzeng.final.REF_10.bam.depth.txt',stringsAsFactors = F)

over100=a[1,]

a=a[-1,]

a=rbind(a,over100)

a[,3]= cumsum(a[,2])

a[,4]= 1-a[,3]/sum(a[,2])

plot(a[,4],type = 'l',xaxt="n",lwd=3)

axis(1, at = 1:nrow(a),labels =a[,1] , las=1)

abline(v=10);abline(v=20);abline(v=30)

最后得到的数据如下：第4列，就是不低于多少X的比例

可以看到大于10X测序深度的比例仍然高达92.88%，效果杠杠的！！

这个测序深度累积分布图，就是很多人的重点！！！！（请大家仔细学习我上面的统计脚本和画图方法）

加油吧，骚年！！！！

请扫描以下二维码关注我们，获取直播系列的所有帖子！

参考：

https://www.biostars.org/p/104063/

http://www.gettinggeneticsdone.com/2014/03/visualize-coverage-exome-targeted-ngs-bedtools.html