前面搜集了一些计算测序覆盖度和深度的方法，这次详细介绍下samtools如何解决了我的问题。
覆盖度和深度
测序深度是指测序得到的总碱基数与待测基因组大小的比值，可以理解为基因组中每个碱基被测序到的平均次数。测序深度 = reads长度 × 比对的reads数目 / 参考序列长度。假设一个基因大小为2M，测序深度为10X，那么获得的总数据量为20M。
覆盖度是指测序获得的序列占整个基因组的比例。指的是基因组上至少被检测到1次的区域，占整个基因组的比例。当然，有些文章中也会将测序深度称为Coverage，容易给我们带来混淆。
参考文章：NGS基础概念-depth and coverage
问题
我们有500多个SNP位点，交给公司设计panel，样本测序后想检测这500多个SNP位点是否被测到，以及这些位点的测序的深度。目前只知道SNP位点，并不知道设计的panel具体区间是什么。
尝试1：Qualimap
Qualimap是功能比较全的一款质控软件，可以对bam文件，RNA-seq，Counts数据质控，也支持比对数据，counts数据和表观数据的比较和聚类，并提供了GUI界面和命令行界面。
- Qualimap：http://qualimap.bioinfo.cipf.es/doc_html/index.html
- 命令行文档：http://qualimap.bioinfo.cipf.es/doc_html/command_line.html#command-line
下载安装Qualimap
wget https://bitbucket.org/kokonech/qualimap/downloads/qualimap_v2.2.1.zip
unzip https://bitbucket.org/kokonech/qualimap/downloads/qualimap_v2.2.1.zip
cd qualimap_v2.2.1
./qualimap
# 安装相关的R包
# NOISeq，Repitools, Rsamtools, GenomicFeatures, rtracklayer
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("NOISeq", version = "3.8")
BiocManager::install("Repitools", version = "3.8")
BiocManager::install("Rsamtools", version = "3.8")
BiocManager::install("GenomicFeatures", version = "3.8")
BiocManager::install("rtracklayer", version = "3.8")
使用
qualimap bamqc -bam file.bam -outdir qualimap_results -outformat pdf
我先是用的上面的命令，算出来结果只有每条染色体的平均coverage，而且非常低。
再一想上面用的命令是不对的，那是针对全基因组，参考序列长度是基于全基因组的长度，所以得出的coverage才如此低。而我的数据是靶向测序，只是全基因组上一些片段，所以我必须知道靶向的区域，才能计算coverage。问题是我只知道关心的SNP位点，并不知道公司设计的靶向区域是什么。如果知道靶向区域，是可以做个bed文件，然后加上
-gff参数。
- -gff,--feature-file <arg> Feature file with regions of interest in
GFF/GTF or BED format
那么有没有什么方法可以统计单个位点的测序深度呢？
尝试2：DepthOfCoverage
看了下GATK的DepthOfCoverage有个-L参数是可以提供靶向区域，除此之外还有个omitDepthOutputAtEachBase参数，还是解决不了我的问题。
java -Xmx4g -jar GenomeAnalysisTK.jar
-R ref_genome.fasta
-T DepthOfCoverage
--omitDepthOutputAtEachBase
--omitIntervalStatistics
--omitLocusTable
-ct 0 -ct 1 -ct 5 -ct 10 -ct 15 -ct 20 -ct 25 -ct 30
--nBins 99 --start 1 --stop 100
-I sample.bam
-o sample.depth
尝试3：samtools mpileup/depth
samtools mpileup/depth在直播我的基因组里面讲过，我最初尝试的也是这个方法，但是统计方向错了，一是我明明关注的是单个位点的情况，我却统计每条染色体的覆盖度；二是用错了参考序列的长度，最终导致算的覆盖度很低，刚开始没明白为什么。
samtools depth是可以直接统计每个位点的覆盖深度，输出结果包括三列:
- reference name
- postion
- coverage depth
使用
samtools depth final.bam > all_position_depth
all_position_depth
然后从all_position_depth文件提取出所关心的SNP位点即可：
final_depth
总结
全基因组、全外显子组或者靶向测序计算覆盖度和深度时一定要注意参考序列长度的选择。这一点在直播我的基因组系列教程（ 「直播」我的基因组（21）：为什么我算出的染色体覆盖度与公司差异甚大）也有提到，健明师兄最初算的染色体长度是全部碱基的长度，包含碱基N，而公司统计的是去除碱基N的，还有就是他们对Y染色体的重复序列在比对分析的时候也mask掉了。
而我的是靶向测序，最初统计覆盖度时没有考虑到参考序列长度，实际计算的结果是以全长基因组为参考序列长度，所以结果就出奇的低了。
另外如果是只关心某些位点是否被测到，以及测到的深度，可以直接用samtools depth。
走的这么多弯路归根结底还是对概念、原理理解的不透彻，今后要加强这方面的学习。
网友评论