【直播】我的基因组77：批量计算每个蛋白编码基因的测序深度及覆盖度

目前我使用的仍然是hg19系统的参考基因组，所以就在gencode数据库里面下载了基于hg19的gtf注释文件，并格式化如下：

head ~/reference/gtf/gencode/protein_coding.hg19.position

我们论坛有专门的教程讲解如何格式化，得到每个基因组的起始终止坐标，就不在此赘述啦（根据gtf格式的基因注释文件得到人所有基因的染色体坐http://www.biotrainee.com/thread-472-1-1.html）

然后我们可以利用大名鼎鼎的bedtools对所有基因批量算出GC含量，如下：

bedtools nuc -fi ~/reference/genome/hg19/hg19.fa -bed   ~/reference/gtf/gencode/protein_coding.hg19.position |cut -f 1-4,5,13 >hg19.protein_coding.gc.txt results.txt

基因长度只是个附带品咯，因为我本来就有基因的起始终止坐标，所以说长度几乎等于是已知的咯。

bedtools的nuc命令还有给出其它信息，我们并不需要，就取第5,13列即可，基本的shell语法大家需要自己学一点，别看了我的直播这么久，还问那些基础问题。

之前我们讲过samtools的depth用法，很容易就可以根据我们拿到的基因起始终止坐标信息来批量依次提取每个基因的被测序的长度，平均测序深度，还有平均测序深度的方差！

脚本如下：

cat ~/reference/gtf/gencode/protein_coding.hg19.position  |while read id
do
arr=($id)
echo ${arr[0]}:${arr[1]}-${arr[2]}
echo -n ${arr[3]} >>results.txt
echo -n -e 't'  >>results.txt
samtools depth -r  ${arr[0]}:${arr[1]}-${arr[2]} ../bamFiles/jmzeng.filter.rmdup.bam | 
perl -alne '{$h{$.}=$F[2];$sum+=$F[2]}END{if($.>0){$mean=$sum/$.;$cum+=($_-$mean)**2 foreach values %h;$cum/=$.;print "$.t$meant$cum"}else{print "0t0t0"}}' >>results.txt
done

这个脚本很简单，主要是对samtools的depth的输入进行简单的统计而已。

我们可以从统计的结果看到有的基因覆盖度极高，但有的基因覆盖度却很低，这是为什么呢？下一讲我们就简单的解析一下蛋白编码基因的测序深度以及覆盖度吧！