前面我们在第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 跑完之后,对每条染色体都会输出如下文件: over100110789 027065 11730286 22219409 32728526 43251046 53774335 64303971 ~~~~~~~~~~~~~~~~~~~~~~~ 后面省略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%,效果杠杠的!!
加油吧,骚年!!!! 参考: http://www./2014/03/visualize-coverage-exome-targeted-ngs-bedtools.html 文:Jimmy、阿尔的太阳 图文编辑:吃瓜群众
|
|