1.bedtools 文檔
$ bedtools --version
bedtools v2.31.1coverage Compute the coverage over defined intervals.
Usage:bedtools coverage [OPTIONS] -a <FILE> \-b <FILE1, FILE2, ..., FILEN>(or):coverageBed [OPTIONS] -a <FILE> \-b <FILE1, FILE2, ..., FILEN>注意: As of version 2.24.0, the coverage tool has changed such that the coverage is computed for the A file, not the B file. This changes the command line interface to be consistent with the other tools.
從2.24.0版本開始，計算的是A文件的覆蓋度，而不是B文件。
重要參數:
- https://bedtools.readthedocs.io/en/latest/content/tools/coverage.html
## bedtools coverage -a positions.bed -b aligned_reads.bam -d > position_depth.txt
-a BAM/BED/GFF/VCF file “A”. Each feature in A is compared to B in search of overlaps. Use “stdin” if passing A with a UNIX pipe.
-b One or more BAM/BED/GFF/VCF file(s) “B”. Use “stdin” if passing B with a UNIX pipe. NEW!!!: -b may be followed with multiple databases and/or wildcard (*) character(s).
-d Report the depth at each position in each A feature. Positions reported are one based. Each position and depth follow the complete A feature.
2.測試
(1) 準備bed文件
$ cd /home/wangjl/tmp找到一個IGV的峰圖范圍
$ cat positions.bed
chr1 246766959 246768109
(2)運行
$ bedtools coverage -a positions.bed -b /datapool/wangjl/scPolyA-seq2/chenxi/PBMC/cluster_bam/RNA/RNA_cluster0.pA_reads.sorted.bam -d > position_depth.txt
# 20:09->20:11
警告:***** WARNING: File /datapool/wangjl/scPolyA-seq2/chenxi/PBMC/cluster_bam/RNA/RNA_cluster0.pA_reads.sorted.bam has inconsistent naming convention for record:GL000008.2 47002 47037 A01758:191:HGG7FDSX5:3:1561:3830:23844_AGTAAACC_TCCACCGTAT 255 -***** WARNING: File /datapool/wangjl/scPolyA-seq2/chenxi/PBMC/cluster_bam/RNA/RNA_cluster0.pA_reads.sorted.bam has inconsistent naming convention for record:GL000008.2 47002 47037 A01758:191:HGG7FDSX5:3:1561:3830:23844_AGTAAACC_TCCACCGTAT 255 -
(3)查看輸出文件
可見，前三列和輸入的-a bed文件一致，第四列是逐個堿基位置編號，第五列是該位置的測序深度。
246768109-246766959=1150, 可見，這是一個半開半閉區間。比如 1-2 只算1個堿基，到底哪邊開呢？一般是 [)。在 (5)中我們做測試。
$ wc position_depth.txt1150 5750 36843 position_depth.txt
$ head position_depth.txt
chr1 246766959 246768109 1 17
chr1 246766959 246768109 2 17
chr1 246766959 246768109 3 17
...
$ tail position_depth.txt
...
chr1 246766959 246768109 1148 19
chr1 246766959 246768109 1149 19
chr1 246766959 246768109 1150 19
(4)R繪圖，和IGV峰圖比較
按照x=第4列,y=5th畫圖
$ R
> dat=read.table("position_depth.txt")
> plot(dat$V4, dat$V5, type="o", cex=0.5)
縮小到長寬基本一致，看起來還是挺像IGV峰圖的。
圖1：上圖是IGV圖，下圖是逐位點的測序深度連線圖。
(5)驗證 bedtools 對bed是取左閉右開區間
在IGV中找到最高峰附近的1bp
目測: 8 上有reads，9上沒有。
$ cat positions2.bed
chr1 246767468 246767469$ bedtools coverage -a positions2.bed -b /datapool/wangjl/scPolyA-seq2/chenxi/PBMC/cluster_bam/RNA/RNA_cluster0.pA_reads.sorted.bam -d > position_depth2.txt$ cat position_depth2.txt
chr1 246767468 246767469 1 15
說明1.單個堿基，2.有reads。支持[)區間。
2)再次測試：chr1 246767468 246767470
輸出:
$ cat position_depth2.txt
chr1 246767468 246767470 1 15
chr1 246767468 246767470 2 15
失敗，竟然都有reads?!
難道，是取了滑動平均值？取平均就不是按位點取測序深度了。
3)再試，爭取測試到有reads和無reads的分界線
發現往后延申坐標不行，都是15。而往前延申坐標可以！
$ cat positions2.bed
chr1 246767458 246767470
輸出:
$ cat position_depth2.txt
chr1 246767458 246767470 1 80
chr1 246767458 246767470 2 78
chr1 246767458 246767470 3 60
chr1 246767458 246767470 4 53
chr1 246767458 246767470 5 53
chr1 246767458 246767470 6 53
chr1 246767458 246767470 7 53
chr1 246767458 246767470 8 53
chr1 246767458 246767470 9 46
chr1 246767458 246767470 10 46 #這個位置應該是58+10-1=67
chr1 246767458 246767470 11 15
chr1 246767458 246767470 12 15
推測：15就是本該是0的背景噪音了？
結論:
- i) bedtools取bed區間時，按照左閉右開區間。
- ii) bed記錄的坐標是0-based，就是igv(1-based)看到的坐標減去1。
推理: 我想要某個點的 pos=8 的reads 深度，該怎么設置bed？
- i)bedtools 就需要設置 [pos, pos+1)位置；
- ii)寫成bed格式要坐標-1，就是 chr pos-1 pos，也就是 chr 7 8
驗證1: 取pos=8的深度，igv有
$ cat positions2.bedchr1 246767467 246767468輸出:$ cat position_depth2.txtchr1 246767467 246767468 1 46
驗證2: 取pos=9的深度，igv看沒
$ cat positions2.bedchr1 246767467 246767469輸出:$ cat position_depth2.txtchr1 246767467 246767469 1 46chr1 246767467 246767469 2 15
圖2：使用的樣本為 cluster0。IGV坐標中 pos=8有reads，pos=9沒有reads。
結論：IGV和位點深度量化相對大小是一致的，絕對大小不一致。
- IGV 最高點是38，最低點是0。圖1中IGV最高也才到66，還是低于(5)-3中的最高值80。
- IGV繼續放大pos=9，還是有零星2個reads的，也比bedtools給出的低。
說明：IGV可能做了取子集？只是展示了bam的一部分？
)
)
)
