统计测序深度和覆盖度的工具

1.GATK DepthOfCoverage

cat test.sh 

#!/usr/bin/bash
HG38=/path/hg38/hg38.fa
GATK=/path/biosoft/gatk3.7/GenomeAnalysisTK.jar

java -jar -Xmx10g $GATK  -T DepthOfCoverage 
-R $HG38 
-o ./test 
-I /path2bam/03_bam_processing/03_base_recal/test.sorted_MarkDuplicates_BQSR.bam 
-L /path2bed/target.bed 
--omitDepthOutputAtEachBase --omitIntervalStatistics 
-ct 1 -ct 5 -ct 10 -ct 20 -ct 30 -ct 50 -ct 100

结果生成四个文件

 4096 Jan 13 21:39 ./
 4096 Jan 13 21:05 ../
 6417 Jan 13 21:39 test.sample_cumulative_coverage_counts
 GenekVIP 6412 Jan 13 21:39 test.sample_cumulative_coverage_proportions
 9362 Jan 13 21:39 test.sample_statistics
 291 Jan 13 21:39 test.sample_summary

　less test.sample_sumary 结果不是很好

sample_id       total   mean    granular_third_quartile granular_median granular_first_quartile %_bases_above_1 %_bases_above_5 %_bases_above_10        %_bases_above_20        %_bases_above_3
test    2025198 17.50   1       1       1       5.1     2.7     2.6     2.5     2.5     2.4     2.3
Total   2025198 17.50   N/A     N/A     N/A

#

2.bamdst: bam文件深度统计

#

软件安装：

git clone https://github.com/shiquan/bamdst.git
cd bamdst/
make

测试：

cat  test.sh

#!/usr/bin/bash
/path2bamdst/biosoft/bamdst/bamdst -p /path2bed/target.bed -o ./ ../03_base_recal/test.sorted.markdup.realign.BQSR.bam
echo "run status $?"

结果文件：

   1630 Jan 13 16:58 chromosomes.report
   4331 Jan 13 16:58 coverage.report
  64188 Jan 13 16:58 depth_distribution.plot
 838584 Jan 13 16:58 depth.tsv.gz
  18119 Jan 13 16:58 insertsize.plot
   9519 Jan 13 16:58 region.tsv.gz
     0  Jan 13 16:58 uncover.bed

看一下内容：

$ less -S chromosomes.report#该文件中存储的是从bam文件中获取的染色体深度、覆盖度信息
Chromosome        DATA(%)       Avg depth          Median       Coverage%        Cov 4x %       Cov 10x %       Cov 30x %
chrM       100.00            0.26              0.0           5.66            2.83            0.00            0.00
$ less coverage.report#信息太多了，我目前觉得比较重要的有
[Total] Raw Reads (All reads)    15
[Total] Mapped Reads    15
[Total] Properly paired    0#Paired reads with properly insert size
[Total] Fraction of Properly paired    0.00%
[Total] Read and mate paired    0#read1 and read2 paired.
[Total] Fraction of Read and mate paired    0.00%
[Total] Map quality cutoff value    20
[Total] MapQuality above cutoff reads    10
[Total] Fraction of MapQ reads in all reads    66.67%
[Total] Fraction of MapQ reads in mapped reads    66.67%
[Target] Average depth    0.26
[Target] Average depth(rmdup)    0.06
[Target] Coverage (>0x)    5.66%
[Target] Coverage (>=4x)    2.83%
[Target] Coverage (>=10x)    0.00%
[Target] Coverage (>=30x)    0.00%
$ less depth_distribution.plot#结合R可以做出目标区域的深度分布图
0       900     0.943396        54      0.056604
1       0       0.000000        54      0.056604
2       0       0.000000        54      0.056604
3       27      0.028302        27      0.028302
4       4       0.004193        23      0.024109
#左边三列分别代表：覆盖深度，对应该深度的碱基数，碱基占比（以目标区域的碱基数为分母）；
#右边两列是高深度向低深度累加的值，分别是碱基数及其占比，累加到X=1为止。
$ less depth.tsv
#Chr    Pos     Raw Depth       Rmdup depth     Cover depth
chrM    650     8       6       8
chrM    651     8       6       8
chrM    652     8       6       8
chrM    653     9       6       9
chrM    654     9       6       9
#Raw Depth从输入bam文件中得到，没有任何限制；
#Rmdup depth去除了PCR重复的reads, 次优比对reads, 低比对质量的reads(mapQ < 20), 该值与samtools depth的输出深度类似；
#默认使用raw depth来统计coverage.report文件中的覆盖率信息。 如果要使用rmdup depth计算覆盖率，可使用参数"--use_rmdup"；
#The coverage depth is the raw depth with consider of deletion region, so this value should be equal to or greated than the raw depth.
$ less insertsize.plot#做出两个接头之间的插入片段大小的图，理论上是根据read1和read2在染色体上的坐标信息求得，这时read1和read2要求比对到一条染色体上。
$ less region.tsv
#Chr    Start   Stop    Avg depth       Median  Coverage        Coverage(FIX)
chrM    649     1603    0.26    0.0     5.66    5.66
#目标区域文件每一个区域的统计，其中Coverage以X>0来统计
$ less uncover.bed
chrM    672     1603
#--uncover的值默认是5，从前面depth_distribution.plot文件中也能看出来，深度小于5的碱基数就是931；
#包含低覆盖或者是未覆盖的，通过"--uncover"规定“未覆盖”的阈值。

REFERENCE   简书作者：hsy_hzauer 链接：https://www.jianshu.com/p/ac7b8e4d76e4