下載安裝
git clone https://github.com/shiquan/bamdst.git
cd bamdst/
make
安裝完成之后出現(xiàn)了可執(zhí)行文件
使用
先看看example/里面有些啥
有一個(gè)bed區(qū)域文件和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> -
區(qū)域文件和輸出目錄是必要的，并且需要注意bed的格式，字段與字段之間一定是制表符鹅心，幾個(gè)空格是不對(duì)的展辞！
試一試
bamdst -p MT-RNR1.bed -o ./ test.bam
多了7個(gè)結(jié)果文件键科，分別來(lái)看一下
$ less -S chromosomes.report#該文件中存儲(chǔ)的是從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#信息太多了金吗，我目前覺(jué)得比較重要的有
[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#結(jié)合R可以做出目標(biāo)區(qū)域的深度分布圖
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
#左邊三列分別代表：覆蓋深度十兢，對(duì)應(yīng)該深度的堿基數(shù)，堿基占比（以目標(biāo)區(qū)域的堿基數(shù)為分母）摇庙；
#右邊兩列是高深度向低深度累加的值旱物，分別是堿基數(shù)及其占比，累加到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文件中得到宵呛，沒(méi)有任何限制；
#Rmdup depth去除了PCR重復(fù)的reads, 次優(yōu)比對(duì)reads, 低比對(duì)質(zhì)量的reads(mapQ < 20), 該值與samtools depth的輸出深度類(lèi)似夕凝；
#默認(rèn)使用raw depth來(lái)統(tǒng)計(jì)coverage.report文件中的覆蓋率信息宝穗。 如果要使用rmdup depth計(jì)算覆蓋率，可使用參數(shù)"--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#做出兩個(gè)接頭之間的插入片段大小的圖逮矛，理論上是根據(jù)read1和read2在染色體上的坐標(biāo)信息求得，這時(shí)read1和read2要求比對(duì)到一條染色體上转砖。
$ less region.tsv
#Chr Start Stop Avg depth Median Coverage Coverage(FIX)
chrM 649 1603 0.26 0.0 5.66 5.66
#目標(biāo)區(qū)域文件每一個(gè)區(qū)域的統(tǒng)計(jì)须鼎，其中Coverage以X>0來(lái)統(tǒng)計(jì)
$ less uncover.bed
chrM 672 1603
#--uncover的值默認(rèn)是5，從前面depth_distribution.plot文件中也能看出來(lái)堪藐，深度小于5的堿基數(shù)就是931莉兰；
#包含低覆蓋或者是未覆蓋的挑围，通過(guò)"--uncover"規(guī)定“未覆蓋”的閾值礁竞。
可選參數(shù)
#方括號(hào)中的是程序默認(rèn)的值，即使不加這些參數(shù)杉辙，程序也按這個(gè)來(lái)
-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
側(cè)翼序列：指存在于編碼區(qū)第一個(gè)外顯子和最末一個(gè)外顯子的外側(cè)是一段不被翻譯的核苷酸序列模捂。側(cè)翼序列含有基因調(diào)控順序，如5端含有的啟動(dòng)子蜘矢、 增強(qiáng)子狂男，3端含有的終止子和多聚腺苷酸信號(hào)等，對(duì)基因表達(dá)起重要的調(diào)控作用品腹。
練習(xí)
我用全外顯子數(shù)據(jù)測(cè)試岖食，20G的bam和一個(gè)20kb長(zhǎng)度gene的bed，跑了9min舞吭。
#這是bed文件泡垃，后面的錯(cuò)誤就是由此引起
chr19 1205798 1228434
#下面是coverage.report的報(bào)告析珊，和已知不相符
#臨床中用于基因檢測(cè)的全外測(cè)序深度都是幾十X的，
#并且對(duì)于捕獲區(qū)域的覆蓋度也是八九不離十蔑穴，這里為什么會(huì)這么低忠寻？
[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%
#原來(lái)是我的bed文件編輯錯(cuò)了，應(yīng)該按照外顯子來(lái)分區(qū)存和；
#否則奕剃，那些沒(méi)有捕獲的區(qū)域會(huì)嚴(yán)重拉低計(jì)算的深度和覆蓋度
gene的外顯子信息可以在哪兒獲取捐腿？
NCBI查找基因可以找到外顯子纵朋；
下載參考基因組時(shí)會(huì)有g(shù)ff文件，里面有分區(qū)茄袖；
用annovar注釋變異位點(diǎn)的時(shí)候倡蝙，會(huì)用到很多數(shù)據(jù)庫(kù)文件，hg38_refGene.txt绞佩、hg38_knownGene.txt寺鸥、hg38_ensGene.txt里面也能查找到外顯子信息。
（這是我現(xiàn)在查看外顯子的方法品山，更簡(jiǎn)單的方法是什么呢胆建？）
#重新編輯之后的bed文件，這一次運(yùn)行時(shí)間是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個(gè)位點(diǎn)沒(méi)測(cè)到肘交，有點(diǎn)多鞍试亍！
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
#這個(gè)文件的信息就非常有用涯呻，可以用來(lái)評(píng)估捕獲效率凉驻，
#結(jié)合基因自身的結(jié)構(gòu)還能探究其對(duì)捕獲測(cè)序的影響
#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#沒(méi)有捕獲的區(qū)域
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
