Comment on page
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文件格式的具体定义请参考:
​https://genome.ucsc.edu/FAQ/FAQformat.html#format1​
​https://bedtools.readthedocs.io/en/latest/content/general-usage.html​
​bedtools是由犹他⼤学Quinlan实验室开发的一套工具，实现了各种在genomic intervals之间进行比较的功能，是一个非常实用的工具。虽然名字叫做bedtools，但是它对于其他文件格式，如bam，vcf等等也提供了支持。
有很多其他的工具，如R包GenomicRanges等,也提供了类似的功能，请有兴趣的同学自行了解。
二代测序数据处理与sam/bam格式
在mapping一节我们提到，在传统的二代测序的数据分析中，第一步往往都是把测到的reads mapping回基因组。sam(Sequence Alignment/Map)文件是最早由bwa的作者Heng Li设计的一种记录reads mapping信息的文件格式，后来被广泛接受，得到了更普遍的应用。目前各种mapping的工具的结果基本上都是以bam/sam文件格式作为输出的。
sam是一种纯文本文件，bam文件是压缩过的，二进制版本的sam文件。在真实的生物信息分析中，为了节约存储，我们通常以bam文件的形式存储mapping的结果。
sam文件由header部分和数据部分组成。header部分存储了参考序列，reads排列顺序，mapping所使用的命令等信息。数据部分存储了每一条read mapping结果的详细信息。希望具体了解的同学请阅读https://samtools.github.io/hts-specs/SAMv1.pdf​
​samtools是最早由Heng Li开发的，用于处理sam/bam文件的一个工具，目前在二代测序的数据分析被极为广泛的使用。在后续的维护中，samtools中用到的函数又被封装到了C程序库htslib中，方便开发者使用。
除了samtools之外，还有其他很多工具也提供了相似的功能，例如bamtools，GATK中的很多工具，等等。在脚本语言python和R中，也有相应的package用于解析sam/bam文件，如pysam和Rsamtools。有兴趣的同学可以自行了解。
1) 示例文件和环境设置
1a) 示例文件
示例文件请从清华云下载。
genes.bed: E.coli K-12基因在参考基因组NC_000913.3上的位置
#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中提供的，根据clip-seq确定的调控因子CsrA的binding peaks
#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的长度
#cat chrom.bed
NC_000913.3     0       4641652
SRR1573434.sampled.bam: GSE61327中提供的E.coli K-12一个野生型样本用bowtie2 mapping后的结果
1b)环境配置
docker
我们还是交互式的进入先前的docker容器:
docker exec -it bioinfo_tsinghua bash
mkdir samtools-bedtools && cd samtools-bedtools
本节只依赖于bedtools和samtools两个软件。
容器内已经装好了samtools
bedtools我们可以下载预编译好的版本:
#当前目录为: /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
如果网络状况不好，可以在宿主机从清华云下载预编译的版本，通过共享目录在容器内部访问。
如果不想每次都输bedtools前面的一串路径，可以在~/.bashrc中将bedtools所在的目录加入环境变量:
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适用的算法。
如上图所示，还有很多参数可以对输出进行控制，例如：
-wa：相当于从A文件中过滤出和B文件有重叠的interval。
-wa -wb参数：对于每一个overlap，既输出第一个文件的interval的原始坐标，也输出第二个文件的interval的原始坐标。
-v：相当于从A文件中过滤出和B文件没有重叠的interval(类似grep的-v参数)。
-wo: write overlap，在每一行结果末尾统计这个重叠区域有多长
还有其他很多参数，文档有详细介绍，请自行参考
我们从示例文件出发提供下面几个例子:
#-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    .       +
# -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命令。
我们需要注意，bedtools merge强制要求输入需要按chromsome和start position排序。对于没有这样sort过的bed文件，可以用gnu sort先进行排序:
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:
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的例子：
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落在每一个基因上:
# -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列:
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上面:
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文件的一种变体(请参考https://genome.ucsc.edu/goldenPath/help/bedgraph.html)，也是genomic coverage的一种常见的存储形式。bedgraph可以直接load到igv上进行可视化。
利用bedtools genomecov命令，我们可以方便的从bam文件出发计算出genome coverage:
# -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
# 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一节中，我们用到了一个用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，等等。
博德研究所的这个网页对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等等。
# 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进行排序。
#按默认的染色体位置进行排序,-@是线程数
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个数以及所占的比例进行统计
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
作业所需文件可从清华云下载
（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参数。）
1.1 Genome Browser
Next - Part III. NGS DATA ANALYSES
2.RNA-seq
Last modified 8mo ago