bedtools: flexible tools for genome arithmetic and DNA sequence analysis.

usage:    bedtools <subcommand> [options]

The bedtools sub-commands include:


[ Genome arithmetic ]

    intersect     Find overlapping intervals in various ways.

                  求区域之间的交集，可以用来注释peak，计算reads比对到的基因组区域
                  不同样品的peak之间的peak重叠情况。

    window        Find overlapping intervals within a window around an interval.
    closest       Find the closest, potentially non-overlapping interval.

                  寻找最近但可能不重叠的区域

    coverage      Compute the coverage over defined intervals.

                  计算区域覆盖度

    map           Apply a function to a column for each overlapping interval.
    genomecov     Compute the coverage over an entire genome.
    merge         Combine overlapping/nearby intervals into a single interval.

                  合并重叠或相接的区域

    cluster       Cluster (but don't merge) overlapping/nearby intervals.
    complement    Extract intervals _not_ represented by an interval file.

                  获得互补区域

    subtract      Remove intervals based on overlaps b/w two files.

                  计算区域差集

    slop          Adjust the size of intervals.

                  调整区域大小，如获得转录起始位点上下游3 K的区域

    flank         Create new intervals from the flanks of existing intervals.

    sort          Order the intervals in a file.

                  排序，部分命令需要排序过的bed文件

    random        Generate random intervals in a genome.

                  获得随机区域，作为背景集

    shuffle       Randomly redistrubute intervals in a genome.

                  根据给定的bed文件获得随机区域，作为背景集

    sample        Sample random records from file using reservoir sampling.
    spacing       Report the gap lengths between intervals in a file.
    annotate      Annotate coverage of features from multiple files.

[ Multi-way file comparisons ]

    multiinter    Identifies common intervals among multiple interval files.
    unionbedg     Combines coverage intervals from multiple BEDGRAPH files.

[ Paired-end manipulation ]

    pairtobed     Find pairs that overlap intervals in various ways.
    pairtopair    Find pairs that overlap other pairs in various ways.

[ Format conversion ]

    bamtobed      Convert BAM alignments to BED (& other) formats.
    bedtobam      Convert intervals to BAM records.
    bamtofastq    Convert BAM records to FASTQ records.
    bedpetobam    Convert BEDPE intervals to BAM records.
    bed12tobed6   Breaks BED12 intervals into discrete BED6 intervals.

[ Fasta manipulation ]

    getfasta      Use intervals to extract sequences from a FASTA file.

                  提取给定位置的FASTA序列

    maskfasta     Use intervals to mask sequences from a FASTA file.
    nuc           Profile the nucleotide content of intervals in a FASTA file.

[ BAM focused tools ]

    multicov      Counts coverage from multiple BAMs at specific intervals.
    tag           Tag BAM alignments based on overlaps with interval files.

[ Statistical relationships ]

    jaccard       Calculate the Jaccard statistic b/w two sets of intervals.

                  计算数据集相似性

    reldist       Calculate the distribution of relative distances b/w two files.
    fisher        Calculate Fisher statistic b/w two feature files.

[ Miscellaneous tools ]

    overlap       Computes the amount of overlap from two intervals.
    igv           Create an IGV snapshot batch script.

                  用于生成一个脚本，批量捕获IGV截图

    links         Create a HTML page of links to UCSC locations.

    makewindows   Make interval "windows" across a genome.

                  把给定区域划分成指定大小和间隔的小区间 (bin)

    groupby       Group by common cols. & summarize oth. cols. (~ SQL "groupBy")

                  分组结算，不只可以用于bed文件。

    expand        Replicate lines based on lists of values in columns.
    split         Split a file into multiple files with equal records or base pairs.

ct@ehbio:~$ conda install bedtools

ct@ehbio:~$ mkdir bedtools
ct@ehbio:~$ cd bedtools
ct@ehbio:~$ url=https://s3.amazonaws.com/bedtools-tutorials/web
ct@ehbio:~/bedtools$ curl -O ${url}/maurano.dnaseI.tgz
ct@ehbio:~/bedtools$ curl -O ${url}/cpg.bed
ct@ehbio:~/bedtools$ curl -O ${url}/exons.bed
ct@ehbio:~/bedtools$ curl -O ${url}/gwas.bed
ct@ehbio:~/bedtools$ curl -O ${url}/genome.txt
ct@ehbio:~/bedtools$ curl -O ${url}/hesc.chromHmm.bed

ct@ehbio:~/bedtools$ head -n 3 cpg.bed exons.bed

==> cpg.bed <==

chr1    28735   29810   CpG:_116
chr1    135124  135563  CpG:_30
chr1    327790  328229  CpG:_29



==> exons.bed <==

chr1    11873   12227   NR_046018_exon_0_0_chr1_11874_f 0   +
chr1    12612   12721   NR_046018_exon_1_0_chr1_12613_f 0   +
chr1    13220   14409   NR_046018_exon_2_0_chr1_13221_f 0   +

ct@ehbio:~/bedtools$ bedtools intersect -a cpg.bed -b exons.bed | head -5

chr1    29320   29370   CpG:_116
chr1    135124  135563  CpG:_30
chr1    327790  328229  CpG:_29
chr1    327790  328229  CpG:_29
chr1    327790  328229  CpG:_29

ct@ehbio:~/bedtools$ bedtools intersect -a cpg.bed -b exons.bed -wa -wb > | head -5

ct@ehbio:~/bedtools$ bedtools intersect -a cpg.bed -b exons.bed -wo | head -10

ct@ehbio:~/bedtools$ bedtools intersect -a cpg.bed -b exons.bed -c | head

chr1    28735   29810   CpG:_116    1
chr1    135124  135563  CpG:_30 1
chr1    327790  328229  CpG:_29 3
chr1    437151  438164  CpG:_84 0
chr1    533219  534114  CpG:_94 0
chr1    544738  546649  CpG:_171    0
chr1    713984  714547  CpG:_60 1
chr1    762416  763445  CpG:_115    10
chr1    788863  789211  CpG:_28 9

# -names标注注释来源
# -sorted: 如果使用了这个参数，提供的一定是排序好的bed文件

ct@ehbio:~/bedtools$ bedtools intersect -a exons.bed \
    -b cpg.bed gwas.bed hesc.chromHmm.bed -sorted -wa -wb -names cpg gwas chromhmm \
    | head -10000  | tail -10

ct@ehbio:~/bedtools$ bedtools merge -i exons.bed | head -n 5

chr1    11873   12227
chr1    12612   12721
chr1    13220   14829
chr1    14969   15038
chr1    15795   15947

ct@ehbio:~/bedtools$ bedtools merge -i exons.bed -c 1 -o count | head -n 5

chr1    11873   12227   1
chr1    12612   12721   1
chr1    13220   14829   2
chr1    14969   15038   1
chr1    15795   15947   1

# -c: 指定对哪些列进行操作
# -o: 与-c对应，表示对指定列进行哪些操作
# 这里的用法是对第一列做计数操作，输出这个区域是由几个区域合并来的
# 对第4列做收集操作，记录合并的区域的名字，并逗号分隔显示出来

ct@ehbio:~/bedtools$ bedtools merge -i exons.bed -d 340 -c 1,4 -o count,collapse | head -4

chr1    11873   12227   1   NR_046018_exon_0_0_chr1_11874_f
chr1    12612   12721   1   NR_046018_exon_1_0_chr1_12613_f
chr1    13220   15038   3   NR_046018_exon_2_0_chr1_13221_f,NR_024540_exon_0_0_chr1_14362_r,NR_024540_exon_1_0_chr1_14970_r
chr1    15795   15947   1   NR_024540_exon_2_0_chr1_15796_r

ct@ehbio:~/bedtools$ head genome.txt 

chr1    249250621
chr10   135534747
chr11   135006516
chr11_gl000202_random   40103
chr12   133851895
chr13   115169878
chr14   107349540
chr15   102531392



ct@ehbio:~/bedtools$ bedtools complement -i exons.bed -g genome.txt | head -n 5

chr1    0   11873
chr1    12227   12612
chr1    12721   13220
chr1    14829   14969
chr1    15038   15795

# 对单行FASTA，可如此计算
# 如果是多行FASTA，则需要累加

ct@ehbio:~/bedtools$ awk 'BEGIN{OFS=FS="\t"}{\
    if($0~/>/) {seq_name=$0;sub(">","",seq_name);} \
    else {print seq_name,length;} }' ../bio/genome.fa | tee ../bio/genome.txt 

chr1    60001
chr2    54001
chr3    54001
chr4    60001

ct@ehbio:~/bedtools$ bedtools genomecov -ibam ../bio/map.sortP.bam -bga \
    -g ../bio/genome.txt | head

# 这个warning很有意思，因为BAM中已经有这个信息了，就不需要提供了

*****
*****WARNING: Genome (-g) files are ignored when BAM input is provided. 
*****

# bedgraph文件，前3列与bed相同，最后一列表示前3列指定的区域的覆盖度。

chr1    0   11  0
chr1    11  17  1
chr1    17  20  2
chr1    20  31  3
chr1    31  36  4
chr1    36  43  6
chr1    43  44  7
chr1    44  46  8
chr1    46  48  9
chr1    48  54  10

ct@ehbio:~/bedtools$ bedtools jaccard \
    -a fHeart-DS16621.hotspot.twopass.fdr0.05.merge.bed \
    -b fHeart-DS15839.hotspot.twopass.fdr0.05.merge.bed

intersection    union-intersection  jaccard n_intersections
81269248    160493950   0.50637 130852

ct@ehbio:~/bedtools$ for i in *.merge.bed; do \
    for j in *.merge.bed; do \
    bedtools jaccard -a $i -b $j | cut -f3 | tail -n +2 | sed "s/^/$i\t$j\t/"; \
    done; done >total.similarity

root@ehbio:~# yum install parallel.noarch
# parallel 后面双引号("")内的内容为希望用parallel执行的命令，
# 整体写法与Linux下命令写法一致。
# 双引号后面的 三个相邻冒号 (:::)默认用来传递参数的，可多个连写。
# 每个三冒号后面的参数会被循环调用，而在命令中的引用则是根据其出现的位置，分别用{1}, {2}
# 表示第一个三冒号后的参数，第二个三冒号后的参数。
#
# 这个命令可以替换原文档里面的整合和替换, 相比于原文命令生成多个文件，这里对每个输出结果
# 先进行了比对信息的增加，最后结果可以输入一个文件中。
#

ct@ehbio:~/bedtools$ parallel "bedtools jaccard -a {1} -b {2} | awk 'NR> | cut -f 3 \
    | sed 's/^/{1}\t{2}\t/'" ::: `ls *.merge.bed` ::: `ls *.merge.bed`  >totalSimilarity.2

# 上面的命令也有个小隐患，并行计算时的输出冲突问题，可以修改为输出到单个文件,再cat到一起
ct@ehbio:~/bedtools$ parallel "bedtools jaccard -a {1} -b {2} | awk 'NR> | cut -f 3 \
    | sed 's/^/{1}\t{2}\t/' >{1}.{2}.totalSimilarity_tmp" ::: `ls *.merge.bed` ::: `ls *.merge.bed` 
ct@ehbio:~/bedtools$ cat *.totalSimilarity_tmp >totalSimilarity.2



# 替换掉无关信息

ct@ehbio:~/bedtools$ sed -i -e 's/.hotspot.twopass.fdr0.05.merge.bed//' \
    -e 's/.hg19//' totalSimilarity.2  

fMusle_leg-DS19115	fMusle_back-DS18454	0.55
fMusle_leg-DS19115	fHeart-DS15643	0.4
fHeart-DS15643	fHeart-DS16621	0.8
fHeart-DS15643	fHeart-DS15839	0.7
fHeart-DS16621	fHeart-DS15839	0.7

# 这里面b和c可以用一个，因为是一个对称阵
# 如果2和3列内容不同，此脚本也可用
ct@ehbio:~/bedtools$ awk 'BEGIN{OFS=FS="\t"}{a[$1, $2]=$3; b[$1]=1; c[$2]=1;}END\
	{printf("ID"); for(i in c) printf("\t%s",  i); \
		for (i in b) {printf("%s", i); for(j in c) {printf("\t%s", a[i, j]);} print "";}}

ct@ehbio:~/bedtools$ parallel "bedtools jaccard -a {1} -b {2} \
      | awk 'NR>1' | cut -f 3 \
      > {1}.{2}.jaccard" \
     ::: `ls *.merge.bed` ::: `ls *.merge.bed`

find . \

    | grep jaccard \
    | xargs grep "" \
    | sed -e s"/\.\///" \
    | perl -pi -e "s/.bed./.bed\t/" \
    | perl -pi -e "s/.jaccard:/\t/" \
    > pairwise.dnase.txt

cat pairwise.dnase.txt \
| sed -e 's/.hotspot.twopass.fdr0.05.merge.bed//g' \
| sed -e 's/.hg19//g' \
> pairwise.dnase.shortnames.txt

awk 'NF==3' pairwise.dnase.shortnames.txt \
| awk '$1 ~ /^f/ && $2 ~ /^f/' \
| python make-matrix.py \
> dnase.shortnames.distance.matrix