# 1.2 bedtools and samtools ## 0) 背景 ### genomic interval的比较 * 用染色体名称，起始位置坐标，终止位置坐标(还可以包括正负链的信息)，我们可以定义基因组上的一段区域(genomic interval/genomic range)。 * 生物信息中用到的基因组数据基本上都可以描述成基因组上的一段区域，以及注释到这一段区域的数据。我们在学习linux文本处理的基本命令时介绍的用于基因组注释的gtf/gff文件，就是这样的一个例子。 * 在很多数据分析的任务中，我们经常需要在不同的genomic intervals之间进行比较，例如： * 取两组genomic intervals的交集 * 将一组genomic intervals中有重合的区域合并到一起 * 从一组genomic intervals中减去另外一组genomic intervals(差集) * 统计一组genomic intervals有多少个落在另外一组genomic intervals上 * 对于一组genomic intervals，统计落在基因组每一个位置的genomic intervals的数量（即计算它们基因组上的coverage） * ... * 高效的实现上述的一些操作需要依赖特殊的数据结构，例如interval tree等。实现一个interval tree一类的数据结构来可靠的支持上述操作不是一件trivial的事情，但是好在我们并不需要自己实现这些功能，用很多经过长期测试的工具可以直接拿来使用。 * 如果genomic intervals是事先按染色体和start position排过序的，上述的一些操作可以更高效的实现。 * 和gtf/gff文件类似，bed文件格式是一种存储基因组的一段区域，以及相关的信息的另外一种非常常见的文件格式。和gtf/gff文件格式不同，bed文件的坐标是从0开始的开区间。 * bed文件格式的具体定义请参考: * * * [bedtools](https://bedtools.readthedocs.io)是由犹他⼤学Quinlan实验室开发的一套工具，实现了各种在genomic intervals之间进行比较的功能，是一个非常实用的工具。虽然名字叫做bedtools，但是它对于其他文件格式，如bam，vcf等等也提供了支持。 * 有很多其他的工具，如R包[GenomicRanges](https://bioconductor.org/packages/GenomicRanges/)等,也提供了类似的功能，请有兴趣的同学自行了解。 ### 二代测序数据处理与sam/bam格式 * 在[mapping](https://book.ncrnalab.org/teaching/part-iii.-ngs-data-analyses/1.mapping)一节我们提到，在传统的二代测序的数据分析中，第一步往往都是把测到的reads mapping回基因组。sam(Sequence Alignment/Map)文件是最早由[bwa](http://bio-bwa.sourceforge.net/)的作者[Heng Li](http://www.liheng.org/)设计的一种记录reads mapping信息的文件格式，后来被广泛接受，得到了更普遍的应用。目前各种mapping的工具的结果基本上都是以bam/sam文件格式作为输出的。 * sam是一种纯文本文件，bam文件是压缩过的，二进制版本的sam文件。在真实的生物信息分析中，为了节约存储，我们通常以bam文件的形式存储mapping的结果。 * sam文件由header部分和数据部分组成。header部分存储了参考序列，reads排列顺序，mapping所使用的命令等信息。数据部分存储了每一条read mapping结果的详细信息。希望具体了解的同学请阅读 * [samtools](http://www.htslib.org/)是最早由Heng Li开发的，用于处理sam/bam文件的一个工具，目前在二代测序的数据分析被极为广泛的使用。在后续的维护中，samtools中用到的函数又被封装到了C程序库[htslib](https://github.com/samtools/htslib)中，方便开发者使用。 * 除了samtools之外，还有其他很多工具也提供了相似的功能，例如[bamtools](https://github.com/pezmaster31/bamtools)，[GATK](https://gatk.broadinstitute.org/hc/en-us/categories/360002369672)中的很多工具，等等。在脚本语言python和R中，也有相应的package用于解析sam/bam文件，如[pysam](https://pysam.readthedocs.io/en/latest/api.html)和[Rsamtools](https://bioconductor.org/packages/Rsamtools)。有兴趣的同学可以自行了解。 ## 1) 示例文件和环境设置 ### 1a) 示例文件 * 示例文件请从[清华云](https://cloud.tsinghua.edu.cn/d/ad22768345664924b202/?p=%2FFiles%2FPART_III%2F1.mapping%2Fsamtools%26bedtools%2Fexample\&mode=list)下载。 * `genes.bed`: E.coli K-12基因在参考基因组NC\_000913.3上的位置 ```bash #head genes.bed NC_000913.3 189 255 thrL . + NC_000913.3 336 2799 thrA . + NC_000913.3 2800 3733 thrB . + NC_000913.3 3733 5020 thrC . + NC_000913.3 5233 5530 yaaX . + NC_000913.3 5682 6459 yaaA . - NC_000913.3 6528 7959 yaaJ . - NC_000913.3 8237 9191 talB . + NC_000913.3 9305 9893 mog . + NC_000913.3 9927 10494 satP . - ``` * `CsrA.peaks.bed`: [GSE102380](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE102380)中提供的，根据clip-seq确定的调控因子CsrA的binding peaks ```bash #head CsrA.peaks.bed NC_000913.3 311 336 . . + NC_000913.3 430 468 . . - NC_000913.3 551 577 . . + NC_000913.3 1001 1040 . . + NC_000913.3 1058 1080 . . + NC_000913.3 1100 1180 . . - NC_000913.3 1302 1329 . . + NC_000913.3 2077 2104 . . - NC_000913.3 2226 2269 . . - NC_000913.3 2262 2286 . . + ``` * `chrom.bed`: 参考基因组NC\_000913.3的长度 ```bash #cat chrom.bed NC_000913.3 0 4641652 ``` * `SRR1573434.sampled.bam`: [GSE61327](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE61327)中提供的E.coli K-12一个野生型样本用bowtie2 mapping后的结果 ### 1b)环境配置 #### docker * 我们还是交互式的进入先前的docker容器: ```bash docker exec -it bioinfo_tsinghua bash mkdir samtools-bedtools && cd samtools-bedtools ``` * 本节只依赖于`bedtools`和`samtools`两个软件。 * 容器内已经装好了`samtools` * `bedtools`我们可以下载预编译好的版本: ```bash #当前目录为: /home/test/samtools-bedtools mkdir bin wget -O bin/bedtools https://github.com/arq5x/bedtools2/releases/download/v2.30.0/bedtools.static.binary chmod u+x bin/bedtools ``` * 如果网络状况不好，可以在宿主机从[清华云](https://cloud.tsinghua.edu.cn/d/ad22768345664924b202/?p=%2FFiles%2FPART_III%2F1.mapping%2Fsamtools%26bedtools\&mode=list)下载预编译的版本，通过共享目录在容器内部访问。 * 如果不想每次都输bedtools前面的一串路径，可以在`~/.bashrc`中将bedtools所在的目录加入环境变量: ```bash echo "export PATH=$PATH:/home/test/samtools-bedtools/bin" >> ~/.bashrc source ~/.bashrc which bedtools # /home/test/samtools-bedtools/bin/bedtools ``` * 请从清华云下载所需的数据，然后在容器内通过共享目录访问。 ## 2) bedtools * 我们接下来介绍怎样用bedtools实现前面提到的一些对于genomic intervals进行比较的任务 ### 2.1) 取交集 * `bedtools intersect`比较在一组genomic intevals和另外一组(或多组genomic inteval)的重叠情况，然后找它们重叠的区域 * 在默认参数下，bedtools的提供的包括intersect在内的多数工具都不考虑strand的信息，输出gemomic interval的重叠部分(至少有1nt重合就认为存在重叠) * `bedtools intersect`提供了很多可选的参数，来对软件的行为进行更精细的控制，例如: * `-s`: 当两个genomic inteval不仅位置上存在重叠，而且还在同一个strand上时，才认为存在重叠 * `-f`: 当两个genomic inteval重叠的比例大于一个给定值时，才认为存在重叠 * `-sorted`: `bedtools intersect`对于输入的genomic interval的次序没有强制的要求，但是如果我们的输入已经按染色体和坐标sort过，就可以通过这个参数告知bedtools使用效率更高，但只对sorted过的interval适用的算法。 ![](https://4115668567-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LPVsf5VZbQ7h14X29qW%2F-Mjwl7PMCtH34muKBlGl%2F-Mjwl_JWBhk1jSZFF2UM%2Fbedtools_1.png?alt=media\&token=0668d9b2-269f-4dc4-9c5c-e46f0b62ae27) * 如上图所示，还有很多参数可以对输出进行控制，例如： * `-wa`：相当于从A文件中过滤出和B文件有重叠的interval。 * `-wa -wb`参数：对于每一个overlap，既输出第一个文件的interval的原始坐标，也输出第二个文件的interval的原始坐标。 * `-v`：相当于从A文件中过滤出和B文件没有重叠的interval(类似grep的-v参数)。 * `-wo`: write overlap，在每一行结果末尾统计这个重叠区域有多长 * 还有其他很多参数，文档有详细介绍，请自行参考 * 我们从示例文件出发提供下面几个例子: ```bash #-s: require intervals to reside in same strand bedtools intersect -s -a genes.bed -b CsrA.peaks.bed | head # NC_000913.3 551 577 thrA . + # NC_000913.3 1001 1040 thrA . + # NC_000913.3 1058 1080 thrA . + # NC_000913.3 1302 1329 thrA . + # NC_000913.3 2262 2286 thrA . + # NC_000913.3 2390 2449 thrA . + # NC_000913.3 2721 2799 thrA . + # NC_000913.3 2800 2805 thrB . + # NC_000913.3 3054 3093 thrB . + # NC_000913.3 4626 4678 thrC . + ``` ```bash # -s: require intervals to reside in same strand # -wa -wb: for each overlap output original entry in file a and file b bedtools intersect -s -a genes.bed -b CsrA.peaks.bed -wa -wb | head # NC_000913.3 336 2799 thrA . + NC_000913.3 551 577 . . + # NC_000913.3 336 2799 thrA . + NC_000913.3 1001 1040 . . + # NC_000913.3 336 2799 thrA . + NC_000913.3 1058 1080 . . + # NC_000913.3 336 2799 thrA . + NC_000913.3 1302 1329 . . + # NC_000913.3 336 2799 thrA . + NC_000913.3 2262 2286 . . + # NC_000913.3 336 2799 thrA . + NC_000913.3 2390 2449 . . + # NC_000913.3 336 2799 thrA . + NC_000913.3 2721 2805 . . + # NC_000913.3 2800 3733 thrB . + NC_000913.3 2721 2805 . . + # NC_000913.3 2800 3733 thrB . + NC_000913.3 3054 3093 . . + # NC_000913.3 3733 5020 thrC . + NC_000913.3 4626 4678 . . + ``` ### 2.2) 合并重叠或相邻的基因组区域 * 有的时候我们会希望把有重合的，或者是相邻的genomic inteval merge到一起。这时候就可以用到`bedtools merge`命令。 ![](https://4115668567-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LPVsf5VZbQ7h14X29qW%2F-Mjwl7PMCtH34muKBlGl%2F-MjwlgqAxncEGzBriBkP%2Fbedtools_2.png?alt=media\&token=015bbfc4-1e8f-45a5-87ed-75c903966733) * 我们需要注意，`bedtools merge`强制要求输入需要按chromsome和start position排序。对于没有这样sort过的bed文件，可以用gnu sort先进行排序: ```bash sort -k1,1 -k2,2n input.bed > input.sorted.bed ``` * `bedtools merge`也提供了很多可选的参数，来对软件的行为进行更精细的控制 * `-d`: distance，只要两个interval距离小于d，就进行merge。 * `-c`结合`-o`参数: 对于需要merge的intervals，将注释到这些intervals的数据按指定的方式整合后作为结果输出。bedtools为这里所谓“指定的方式”提供了很多选项，如sum, min, max, mean,collapse,distinct,count等等。`-c`后面可以指定多列(逗号隔开)，`-o`后面可以对每一列指定相应的操作。所以`-c`结合`-o`参数可以实现一些非常灵活的功能。 * 在下面这个例子中，我们把同一个strand上相邻的基因merge到一起，并统计merge后的每一个interval来源于哪几个基因，来源的基因数量，以及所在的strand: ```bash bedtools merge -s -c 4,1,6 -o collapse,count,distinct -i genes.bed | head # NC_000913.3 189 255 thrL 1 + # NC_000913.3 336 2799 thrA 1 + # NC_000913.3 2800 5020 thrB,thrC 2 + # NC_000913.3 5233 5530 yaaX 1 + # NC_000913.3 5682 6459 yaaA 1 - # NC_000913.3 6528 7959 yaaJ 1 - # NC_000913.3 8237 9191 talB 1 + # NC_000913.3 9305 9893 mog 1 + # NC_000913.3 9927 10494 satP 1 - # NC_000913.3 10642 11356 yaaW 1 - ``` ### 2.3) 差集 * 使用`bedtools subtract`可以从一组genomic intervals中减去另外一组genomic intervals。该命令的很多参数和intersect和merge都非常相似，不再赘述。 * 我们给出一个计算E.coli K-12 integenic region的例子： ```bash bedtools subtract -a chrom.bed -b genes.bed | head # NC_000913.3 0 189 # NC_000913.3 255 336 # NC_000913.3 2799 2800 # NC_000913.3 5020 5233 # NC_000913.3 5530 5682 # NC_000913.3 6459 6528 # NC_000913.3 7959 8237 # NC_000913.3 9191 9305 # NC_000913.3 9893 9927 # NC_000913.3 10494 10642 ``` ### 2.4) 计算count * 不难看出，统计一组genomic intervals有多少个落在另外一组genomic intervals上，实际上是在两组intervals之间做intersection的一个特例。所以这还是可以通过`bedtools intersect`来实现。 * 在下面这个例子中，我们统计有多少个peaks落在每一个基因上: ```bash # -c: For each entry in A, report the number of overlaps with B. bedtools intersect -s -a genes.bed -b CsrA.peaks.bed -c | head # NC_000913.3 189 255 thrL . + 0 # NC_000913.3 336 2799 thrA . + 7 # NC_000913.3 2800 3733 thrB . + 2 # NC_000913.3 3733 5020 thrC . + 1 # NC_000913.3 5233 5530 yaaX . + 0 # NC_000913.3 5682 6459 yaaA . - 1 # NC_000913.3 6528 7959 yaaJ . - 1 # NC_000913.3 8237 9191 talB . + 3 # NC_000913.3 9305 9893 mog . + 2 # NC_000913.3 9927 10494 satP . - 1 ``` * 针对这一目的，bedtools还专门方封装了一个`bedtools coverage`命令，默认会多输出3列: ```bash bedtools coverage -s -a genes.bed -b CsrA.peaks.bed | head # NC_000913.3 189 255 thrL . + 0 0 66 0.0000000 # NC_000913.3 336 2799 thrA . + 7 275 2463 0.1116525 # NC_000913.3 2800 3733 thrB . + 2 44 933 0.0471597 # NC_000913.3 3733 5020 thrC . + 1 52 1287 0.0404040 # NC_000913.3 5233 5530 yaaX . + 0 0 297 0.0000000 # NC_000913.3 5682 6459 yaaA . - 1 44 777 0.0566281 # NC_000913.3 6528 7959 yaaJ . - 1 18 1431 0.0125786 # NC_000913.3 8237 9191 talB . + 3 183 954 0.1918239 # NC_000913.3 9305 9893 mog . + 2 84 588 0.1428571 # NC_000913.3 9927 10494 satP . - 1 28 567 0.0493827 # Meaning of last four columns: # 1) The number of features in B that overlapped the A interval. # 2) The number of bases in A that had non-zero coverage. # 3) The length of the entry in A. # 4) The fraction of bases in A that had non-zero coverage. ``` * bedtools intersect还可以直接从bam文件出发统计有多少个Reads落在每一个指定的interval上面: ```bash bedtools intersect -a genes.bed -b SRR1573434.sampled.bam -c > gene.counts.txt head gene.counts.txt # NC_000913.3 189 255 thrL . + 140 # NC_000913.3 336 2799 thrA . + 925 # NC_000913.3 2800 3733 thrB . + 266 # NC_000913.3 3733 5020 thrC . + 411 # NC_000913.3 5233 5530 yaaX . + 7 # NC_000913.3 5682 6459 yaaA . - 32 # NC_000913.3 6528 7959 yaaJ . - 5 # NC_000913.3 8237 9191 talB . + 489 # NC_000913.3 9305 9893 mog . + 15 # NC_000913.3 9927 10494 satP . - 11 ``` ### 2.5) 计算genome coverage * 对于一组genomic intervals，统计落在基因组每一个位置的genomic intervals的数量（即计算它们基因组上的coverage） * bedgraph格式是bed文件的一种变体(请参考)，也是genomic coverage的一种常见的存储形式。bedgraph可以直接load到igv上进行可视化。 * 利用`bedtools genomecov`命令，我们可以方便的从bam文件出发计算出genome coverage: ```bash # -bg: print output in bedgraph format # -ibam: input is in bam format. input must be sorted by coordinate !!! samtools sort SRR1573434.sampled.bam > SRR1573434.sampled.sorted.bam samtools index SRR1573434.sampled.sorted.bam bedtools genomecov -ibam SRR1573434.sampled.sorted.bam -bg > SRR1573434.coverage.bedgraph head SRR1573434.coverage.bedgraph # NC_000913.3 71 102 1 # NC_000913.3 130 157 2 # NC_000913.3 157 159 11 # NC_000913.3 159 161 12 # NC_000913.3 161 163 11 # NC_000913.3 163 164 21 # NC_000913.3 164 165 52 # NC_000913.3 165 166 57 # NC_000913.3 166 168 60 # NC_000913.3 168 169 62 ``` * 对于真核生物的RNA-seq，我们通常会加上一个`-split`参数，对于intron spanning reads,不把spanning intron的部分计入coverage中。 * 需要注意的是，输入的bam文件**一定应当是按coordinate sort过的**。如果没有sort过，bedtools genomecov并不会报错，但结果会出现问题。对于怎么sort bam文件，我们在samtools一节会详细介绍。 * 从bed文件出发，也可以计算genome coverage ```bash # convert bam to bed file bedtools bamtobed -i SRR1573434.sampled.sorted.bam > SRR1573434.sampled.sorted.bed # get chromosome size cut -f 1,3 chrom.bed > chrom.size # get coverage bedtools genomecov -i SRR1573434.sampled.sorted.bed -g chrom.size -bg > SRR1573434.coverage.2.bedgraph ``` * 在[mapping](https://book.ncrnalab.org/teaching/part-iii.-ngs-data-analyses/1.mapping)一节中，我们用到了一个用perl语言实现的把bam文件里面aligned reads对应的interval转换成bed文件的脚本。这里`bedtools bamtobed`的作用是基本相同的。有些工具只接受bed文件，不接受bam文件，这时bamtobed的转换就能派上用场。当然我们这里只是为了举例，`bedtools genomecov`两种文件格式都可以接受的。 * 从bed文件出发计算genome coverage还需要提供一个记录每一条染色体长度的文件。之所以从bam文件出发不需要这个输入，是因为bam的header里面就有染色体长度的信息。 ## 3) samtools * 本章中我们介绍samtools最常用的几个命令: * view * sort * index * flagstat ### 3.1 ) view * view命令的主要功能是在sam和bam文件之间进行转换，以及根据每一个alignment record的flag，位置，指定tag，mapping质量等多种指标对bam文件的alignment record进行filter。它还提供了对bam文件进行sub sampling等功能。 (1) sam/bam文件互相转换 * bam文件是二进制文件，如果我们用`head`一类的命令查看，看到的东西会像是一堆乱码。如果我们想查看bam文件数据部分的前几行，正确的做法是`samtools view input.bam | head`。 * 下面给出了其他一些例子: ``` # bam to sam # without -h parameter, header will not be printed samtools view -h SRR1573434.sampled.bam > SRR1573434.sampled.sam # sam back to bam # samtools view -b -S SRR1573434.sampled.sam > SRR1573434.sampled.bam ``` (2) sam/bam filtering * bam文件中每一个record的flag用一组二进制数记录了每一个read的mapping状况，例如是否unmapped，mapping到forward strand还是reverse strand；对于paired end reads，当前record对应第一个mate还是第二个mate，于当前alignment的read paired的那条read是否是unmapped，等等。 * 博德研究所的[这个网页](https://broadinstitute.github.io/picard/explain-flags.html)对sam/bam文件flag的定义给出了一个非常好的总结，供大家参考。 * `samtools view`的`-F`,`-f`和`-G`三个参数允许我们根据flag对bam文件进行各种各样的filtering。 * `-f`: INT, only include reads with all of the FLAGs in INT present \[0] * `-F`: INT, only include reads with none of the FLAGS in INT present \[0] * `-G`: INT, only EXCLUDE reads with all of the FLAGs in INT present \[0] * 例如: ``` # discard unmapped reads # flag 4 indicates unmapped samtools view -bF 4 SRR1573434.sampled.bam > SRR1573434.sampled.F.4.bam # extract unmapped reads samtools view -bf 4 SRR1573434.sampled.bam > SRR1573434.sampled.f.4.bam # discard reads pairs with unmapped read # flag 4 indicates unmapped, flag 8 indicates mate unmapped # flag 4+8=12 indicates read itself is unmapped or has unmapped mate samtools view -bF 12 SRR1573434.sampled.bam > SRR1573434.sampled.F.12.bam # extract properly mapped reads samtools view -bf 2 SRR1573434.sampled.bam > SRR1573434.sampled.properly.mapped.bam ``` * samtools view还允许根据除了flag之外的多种指标进行filter，例如mapping quality, bam record的tag等等。 ```bash # extract reads with mapping quality >=4 samtools view -q 4 -O bam SRR1573434.sampled.bam > SRR1573434.sampled.MAQ.gt.4.bam ``` ### 3.2 ) sort * sort对bam文件中的alignment records进行排序。 ```bash #按默认的染色体位置进行排序,-@是线程数 samtools sort -@ 8 SRR1573434.sampled.bam > SRR1573434.sampled.sorted.by.coord.bam -n参数则是根据read名进行排序 samtools sort -n -@ 8 SRR1573434.sampled.bam > SRR1573434.sampled.sorted.by.name.bam ``` * 有很多工具都要求输入的bam文件根据染色体位置进行排序(例如igv，以及我们前面提到的`bedtools genomecov`)，以便根据位置信息进行对mapping到特定位置的reads进行快速访问。 ### 3.3 ) index * index是为了根据基因组的位置对bam文件中mapping到相应位置的alignment records进行快速访问而设计出的一种数据结构。它大致可以理解成基因组上的位置与磁盘相应alignment records所在的chunk对应的offset之间的一个look up table。所以只有当bam文件根据coordinates排序后，才能利用`samtools index`建立这种index。index文件通常有一个`.bai`后缀。 ``` samtools index SRR1573434.sampled.sorted.by.coord.bam ls SRR1573434.sampled.sorted.by.coord.bam.bai ``` ### 3.4 ) flagstat * `flagstat`对不同flag的reads个数以及所占的比例进行统计 ```bash samtools flagstat SRR1573434.sampled.bam 364728 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 secondary 0 + 0 supplementary 0 + 0 duplicates 351059 + 0 mapped (96.25% : N/A) 364728 + 0 paired in sequencing 182364 + 0 read1 182364 + 0 read2 349528 + 0 properly paired (95.83% : N/A) 349988 + 0 with itself and mate mapped 1071 + 0 singletons (0.29% : N/A) 0 + 0 with mate mapped to a different chr 0 + 0 with mate mapped to a different chr (mapQ>=5) ``` ## 4) Homework 作业所需文件可从[清华云](https://cloud.tsinghua.edu.cn/d/ad22768345664924b202/?p=%2FFiles%2FPART_III%2F1.mapping%2Fsamtools%26bedtools%2Fhomework\&mode=list)下载 * （1）我们提供的bam文件`COAD.ACTB.bam`是单端测序分析的结果还是双端测序分析的结果？为什么？(提示：可以使用`samtools flagstat`） * （2）查阅资料回答什么叫做"secondary alignment"？并统计提供的bam文件中，有多少条记录属于"secondary alignment?" （提示：可以使用samtools view -f 获得对应secondary alignment的records进行统计） * （3）请根据`hg38.ACTB.gff`计算出在ACTB基因的每一条转录本中都被注释成intron的区域，以bed格式输出。并提取`COAD.ACTB.bam`中比对到ACTB基因intron区域的bam信息，后将bam转换为fastq文件。 > 提示： > > * 写脚本把ACTB在gff中第三列为"gene"的interval放在一个bed文件中，第三列为"exon"的intervals放在另外一个bed文件中，再使用`bedtools subtract`。 > * 请注意bed文件使用的是0-based coordinate，gff文件使用的是1-based coordinate。 > * 鼓励其他实现方法，描述清楚过程即可 * (4) 利用`COAD.ACTB.bam`计算出reads在ACTB基因对应的genomic interval上的coverage，以bedgraph格式输出。 （提示：对于真核生物转录组测序向基因组mapping得到的bam文件，bedtools genomecov有必要加-split参数。）