下载安装
git clone https://github.com/shiquan/bamdst.git
cd bamdst/
make
安装完成之后出现了可执行文件
使用
先看看example/里面有些啥
有一个bed区域文件和bam文件：MT-RNR1.bed, test.bam
再看看软件的使用方法
bamdst --help
USAGE : bamdst [OPTION] -p <probe.bed> -o <output_dir> [in1.bam [in2.bam ... ]]
or : bamdst [OPTION] -p <probe.bed> -o <output_dir> -
区域文件和输出目录是必要的，并且需要注意bed的格式，字段与字段之间一定是制表符，几个空格是不对的！
试一试
bamdst -p MT-RNR1.bed -o ./ test.bam
多了7个结果文件，分别来看一下
$ 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"规定“未覆盖”的阈值。
可选参数
#方括号中的是程序默认的值，即使不加这些参数，程序也按这个来
-f, --flank [200] flank n bp of each region
-q [20] map quality cutoff value, greater or equal to the value will be count
--maxdepth [0] set the max depth to stat the cumu distribution.
--cutoffdepth [0] list the coverage of above depths
--isize [2000] stat the inferred insert size under this value
--uncover [5] region will included in uncover file if below it
--bamout BAMFILE target reads will be exported to this bam file
-1 begin position of bed file is 1-based
-h, --help print this help info
侧翼序列：指存在于编码区第一个外显子和最末一个外显子的外侧是一段不被翻译的核苷酸序列。侧翼序列含有基因调控顺序，如5端含有的启动子、 增强子，3端含有的终止子和多聚腺苷酸信号等，对基因表达起重要的调控作用。
练习
我用全外显子数据测试，20G的bam和一个20kb长度gene的bed，跑了9min。
#这是bed文件，后面的错误就是由此引起
chr19 1205798 1228434
#下面是coverage.report的报告，和已知不相符
#临床中用于基因检测的全外测序深度都是几十X的，
#并且对于捕获区域的覆盖度也是八九不离十，这里为什么会这么低？
[Target] Average depth 9.48
[Target] Average depth(rmdup) 8.12
[Target] Coverage (>0x) 46.70%
[Target] Coverage (>=4x) 21.99%
[Target] Coverage (>=10x) 14.74%
[Target] Coverage (>=30x) 10.55%
[Target] Coverage (>=100x) 2.45%
#原来是我的bed文件编辑错了，应该按照外显子来分区；
#否则，那些没有捕获的区域会严重拉低计算的深度和覆盖度
gene的外显子信息可以在哪儿获取？
NCBI查找基因可以找到外显子；
下载参考基因组时会有gff文件，里面有分区；
用annovar注释变异位点的时候，会用到很多数据库文件，hg38_refGene.txt、hg38_knownGene.txt、hg38_ensGene.txt里面也能查找到外显子信息。
（这是我现在查看外显子的方法，更简单的方法是什么呢？）
#重新编辑之后的bed文件，这一次运行时间是8min
chr19 1205798 1207203
chr19 1218416 1218500
chr19 1219323 1219413
chr19 1220372 1220505
chr19 1220580 1220717
chr19 1221212 1221340
chr19 1221948 1222006
chr19 1222984 1223172
chr19 1226453 1226663
chr19 1227592 1228435
$ less coverage.report#可以看到有明显的提高了
[Target] Average depth 36.37
[Target] Average depth(rmdup) 31.26
[Target] Coverage (>0x) 64.32%
[Target] Coverage (>=4x) 48.63%
[Target] Coverage (>=10x) 43.96%
[Target] Coverage (>=30x) 40.60%
[Target] Coverage (>=100x) 13.16%
$ less depth_distribution.plot#有1169个位点没测到，有点多啊！
0 1169 0.356838 2107 0.643162
1 336 0.102564 1771 0.540598
2 86 0.026252 1685 0.514347
3 92 0.028083 1593 0.486264
4 98 0.029915 1495 0.456349
5 18 0.005495 1477 0.450855
$ zcat region.tsv.gz | lsx
#这个文件的信息就非常有用，可以用来评估捕获效率，
#结合基因自身的结构还能探究其对捕获测序的影响
#Chr Start Stop Avg depth Median Coverage Coverage(FIX)
chr19 1205798 1207203 22.53 1.0 52.10 52.10
chr19 1218416 1218500 100.05 102.0 100.00 100.00
chr19 1219323 1219413 115.11 121.5 100.00 100.00
chr19 1220372 1220505 76.17 78.0 100.00 100.00
chr19 1220580 1220717 40.50 39.0 100.00 100.00
chr19 1221212 1221340 113.42 111.0 100.00 100.00
chr19 1221948 1222006 46.52 43.0 100.00 100.00
chr19 1222984 1223172 142.72 153.5 100.00 100.00
chr19 1226453 1226663 38.43 39.0 100.00 100.00
chr19 1227592 1228435 1.05 0.0 41.16 41.16
$ less uncover.bed#没有捕获的区域
chr19 1205798 1206763
chr19 1227592 1227638
chr19 1227647 1227654
chr19 1227664 1227680
chr19 1227688 1228435
reference
https://github.com/shiquan/bamdst
https://blog.csdn.net/biubiuv/article/details/40348135
