Platform for Drug Discovery


BEDtools


概要


  • 目的の異なるいくつかのBEDに関係する処理を行うプログラム群です。
  • BEDファイルの処理というよりは、ゲノム上にあるアノテーション情報の和をとったり、差をとったり、積をとったり、配列を切り出したりします。
  • BEDtoolsという名称ですが、殆どのプログラムはBED以外にGFF,VCFも同様に扱えるようです。

注意


  • BEDtoolsは経験上BED6の情報しか考慮していません。BED12フォーマットでは、exon-intron構造や、CDSの開始/終了を定義できますが、これらの情報は無視されます。以下確認済みの
    • subtractBedでは7カラム目以降は全く修正されず、入力そのままが出力されます。
    • fastaFromBedでは開始座標から終了座標が切り出されます。ブロックの開始位置、ブロックの長さは無視されます。

仕様の詳細確認


プログラムリスト


  • #スクリプト名機能概要応用例/利用シーンなど
    1annotateBed複数のBED等のファイルと対象となるBEDファイルを入力として受け取り、対象BEDに対するcoverageのdepth, breadthやcountなどを同時に算出します。
    2bamToBedbamをbedに変換します。
    3bed12ToBed6bed12をbed6に変換します。ブロックサイズが複数ある場合は、複数行に分割されます。
    4bedToBambedをbamに変換します。
    5bedToIgv
    6closestBed領域Aの近傍の領域Bを探します。BED形式のChIPピークの近傍の遺伝子を探索するなど。
    7complementBedゲノムを全領域とした領域Aの補集合を作る?
    8coverageBedcoverageを算出します。bamを変換したbedを、遺伝子領域のbedでカウントしたり、exomeのターゲット領域のcoverageを算出する等利用する機会が多いです。
    9fastaFromBedbedで定義された開始座標、終了座標を元に、fasta配列を切り出します。遺伝子領域の塩基配列や、ChIP-seqでピークコールされた領域の塩基配列を纏めて切り出すなど。
    10flankBed入力のBEDの上流、下流の一定領域情報を返します。元の領域は含まず、上流、下流がそれぞれ別のエントリーとなります。
    11genomeCoverageBedゲノムに対するcoverageを算出します。whole genomeのリシーケンシング用途など。
    12intersectBed領域A基準に領域Bとの重複領域を返します。
    積集合になると思います。領域が2つに分かれる場合は行も2行になります。
    13linksBed
    14maskFastaFromBed入力FASTAをBEDで与えられた領域でマスクして出力します。ゲノム配列にリピート領域の座標を与えてマスクするなど。
    15mergeBed複数のbedで指定された領域の重複を纏めて一つにします。
    16multiBamCov
    17nucBedターゲットとなるfastaの配列ファイルの、bedで指定された領域の塩基組成の統計量を返します。ゲノムの特定領域の塩基組成を調べる場合など。
    18overlapあるファイルに含まれる1行ごとに存在する
    2領域(A,B)の開始(A)、終了(A)、開始(B)、終了(B)をもとに、
    重複塩基数(正の数)か距離(負の数)を返す。
    pair endのBAMファイル中で、insertサイズを算出するなど。
    19pairToBed
    20pairToPair
    21shuffleBed
    22slopBed入力のBEDの上流、下流の一定領域情報を返します。出力には元の領域を含む点で、flankBedと異なります。
    23sortBed領域のサイズや、染色体番号ごとにサイズや、スコアで並べ替える等。
    24subtractBed領域Aから領域Bと重なる部分を除きます。2つ以上に分断される場合は行が分かれます。
    25tagBam
    26unionBedGraphs
    27windowBed

使い方


coverageBedの使い方


  • coverageBedは、BEDフォーマットで記述されたゲノム上の領域情報に対して、coverageを算出します。
  • toolPath/BEDTools-Version-2.12.0/bin/coverageBed -a [カウント対象].bed -b [領域].bed
    • bedフォーマットのマッピング結果等のカバー率を計測する場合
  • toolPath/BEDTools-Version-2.12.0/bin/coverageBed -abam [カウント対象].bam -b [領域].bed
    • bamフォーマットのマッピング結果等のカバー率を計測する場合
  • toolPath/BEDTools-Version-2.12.0/bin/coverageBed -abam [カウント対象].bam -b [領域].bed -s
    • ストランドを考慮。
  • デフォルトオプションでの出力内容は以下
    • #説明備考
      0Bのファイル内容がそのまま出力されます。chr1 16762998 16812569 NBPF1 1 -
      1B領域に重なるAの件数。4955
      2Aのカバー率がゼロでないB領域の塩基数14242
      3B領域の長さ49571
      4Aのカバー率がゼロでないB領域の塩基の割合。0.2873051=14242/49571
  • "-d"オプションの効果
    • 使い方
      toolPath/BEDTools-Version-2.12.0/bin/coverageBed -abam [カウント対象].bam -b [領域].bed -d
    • 塩基単位で被覆数を出力します。
    • 出力例
      chr1 67051160 67163158 WDR78 1 - 1 0
      
      chr1 67051160 67163158 WDR78 1 - 2 0
      
      chr1 67051160 67163158 WDR78 1 - 3 0
      
      :
      
      chr1 8335050 8800286 RERE 1 - 465234 0
      
      chr1 8335050 8800286 RERE 1 - 465235 0
      
      chr1 8335050 8800286 RERE 1 - 465236 0
      
      chr1 16762998 16812569 NBPF1 1 - 1 76
      
      chr1 16762998 16812569 NBPF1 1 - 2 76
      
      chr1 16762998 16812569 NBPF1 1 - 3 76
      
      chr1 16762998 16812569 NBPF1 1 - 4 78
      
      chr1 16762998 16812569 NBPF1 1 - 5 78
      
      chr1 16762998 16812569 NBPF1 1 - 6 79
      
      :
      
      chr1 16762998 16812569 NBPF1 1 - 49569 44
      
      chr1 16762998 16812569 NBPF1 1 - 49570 41
      
      chr1 16762998 16812569 NBPF1 1 - 49571 34
      
      chr1 41748270 42157083 HIVEP3 1 - 1 1
      
      chr1 41748270 42157083 HIVEP3 1 - 2 0
      
      chr1 41748270 42157083 HIVEP3 1 - 3 0
      
      :
      
      chr1 184532027 184550311 PRG4 1 + 18282 0
      
      chr1 184532027 184550311 PRG4 1 + 18283 0
      
      chr1 184532027 184550311 PRG4 1 + 18284 0
    • Image
      • 図は、"NBPF1"の被覆率をRのbarplotで描画したもの。長さが短ければ、Excelの棒グラフでも描画は可能です。※32,000がExcel棒グラフの限度。
  • "-hist"オプションの効果
    • 使い方
      toolPath/BEDTools-Version-2.12.0/bin/coverageBed -abam [カウント対象].bam -b [領域].bed -hist
    • 出力例
      chr1 67051160 67163158 WDR78 1 - 0 105033 111998 0.9378114
      
      chr1 67051160 67163158 WDR78 1 - 1 3413 111998 0.0304738
      
      chr1 67051160 67163158 WDR78 1 - 2 1541 111998 0.0137592
      
      :
      
      chr1 67075871 67163158 WDR78 1 - 27 17 87287 0.0001948
      
      chr1 67075871 67163158 WDR78 1 - 28 4 87287 0.0000458
      
      chr1 67075871 67163158 WDR78 1 - 29 3 87287 0.0000344
      
      chr1 8335050 8406334 RERE 1 - 0 56798 71284 0.7967847
      
      chr1 8335050 8406334 RERE 1 - 1 5194 71284 0.0728635
      
      :
      
      chr1 8335050 8800286 RERE 1 - 105 1 465236 0.0000021
      
      chr1 8335050 8800286 RERE 1 - 107 16 465236 0.0000344
      
      chr1 8335050 8800286 RERE 1 - 108 5 465236 0.0000107
      
      :
      
      chr1 184532027 184550311 PRG4 1 + 85 8 18284 0.0004375
      
      chr1 184532027 184550311 PRG4 1 + 86 12 18284 0.0006563
      
      all 0 3492689 3736822 0.9346683
      
      all 1 114083 3736822 0.0305294
      
      all 2 46056 3736822 0.0123249
      
      :
      
      all 951 1 3736822 0.0000003
      
      all 957 17 3736822 0.0000045
      
      all 959 6 3736822 0.0000016
      
      all 960 17 3736822 0.0000045
    • 以下のように、ある特定のB領域(遺伝子など)ごとに被覆数と頻度を示す情報が出力されます。
    • #説明備考
      0~4上記と同じ
      5被覆数0
      6頻度35329=49571-14242 被覆数が0の塩基数
      7全数49571
      8相対頻度0.7126949=35329/49571
    • Image
      • 実際の出力では、1塩基ごとに出力されます。0件の場合は行が出力されないため、後半は被覆数がとびとびになります。

fastaFromBedの使い方


説明

  • fastaFromBedは、BEDフォーマットで記述されたゲノム上のアノテーション情報をもとに、配列を切り出す処理を行います。
  • BED12形式でアノテーション情報を指定してもexon/intron構造は考慮されず、全長(exon+intron)が抜き出されます。
  • -nameオプションの効果
    • 入力のBED
      chr20   62551017        62551133        MIR941-1##NR_030637     0       +       62551133  62551133        0       2       32,57,  0,59,
    • -nameあり
      >MIR941-1##NR_030637
      
      CCCGGCTGTGTGGACATGTGCCCAGGGCCCGGGACAGCGCCACGGAAGAGGACGCACAGGACAGCGCCACGGAAGAGGACGCACCCGGCTGTGTGCACATGTGCCCAGGGCCCGGG
    • -nameなし
      >chr20:62551017-62551133
      
      CCCGGCTGTGTGGACATGTGCCCAGGGCCCGGGACAGCGCCACGGAAGAGGACGCACAGGACAGCGCCACGGAAGAGGACGCACCCGGCTGTGTGCACATGTGCCCAGGGCCCGGG
  • -sオプションの効果
    • -sオプションでstrandの考慮をしますが、BEDファイルにstrandのカラム(第6カラム)が無くてもエラーは出ず適切に処理されます。
    • 入力のBED
      chr20   62551017        62551133        MIR941-1##NR_030637     0       +       62551133  62551133        0       2       32,57,  0,59,
      
      chr12   45580873        45581252        MIR1975##NR_031739      0       -       45581252  45581252        0       2       26,20,  0,359,
    • -sなし
      >MIR941-1##NR_030637
      
      CCCGGCTGTGTGGACATGTGCCCAGGGCCCGGGACAGCGCCACGGAAGAGGACGCACAGGACAGCGCCACGGAAGAGGACGCACCCGGCTGTGTGCACATGTGCCCAGGGCCCGGG
      
      >MIR1975##NR_031739
      
      caatgttaaatcagcttaacaataacccaTGTTgggcacagtggcttacgcctataatcccagcactttgggaggccgaggtgggtggatcgcctaaggtcaggagttccagaccagcctgaccaacatggtgaaaccccgtctctactaaaaatacaaaaattagccaggcgtggtggtgtgcacctgtagtcacagctactcaggaggctgagacaggagaattgcttgaacccaggaggcagaggttgcagcgatccaagattgtgccattgcactccagcctgggcgacagagcaaggctccatctcaaaacaaacaaacaaacgaacaaaaaaGCAAAACAAAAAAACAATAACCCACAACACTCGGACCAACt
    • -sあり
      >MIR941-1##NR_030637
      
      CCCGGCTGTGTGGACATGTGCCCAGGGCCCGGGACAGCGCCACGGAAGAGGACGCACAGGACAGCGCCACGGAAGAGGACGCACCCGGCTGTGTGCACATGTGCCCAGGGCCCGGG
      
      >MIR1975##NR_031739
      
      aGTTGGTCCGAGTGTTGTGGGTTATTGTTTTTTTGTTTTGCttttttgttcgtttgtttgtttgttttgagatggagccttgctctgtcgcccaggctggagtgcaatggcacaatcttggatcgctgcaacctctgcctcctgggttcaagcaattctcctgtctcagcctcctgagtagctgtgactacaggtgcacaccaccacgcctggctaatttttgtatttttagtagagacggggtttcaccatgttggtcaggctggtctggaactcctgaccttaggcgatccacccacctcggcctcccaaagtgctgggattataggcgtaagccactgtgcccAACAtgggttattgttaagctgatttaacattg
      • 第6カラムのstrandが'-'のMIR1975##NR_031739だけがcomplementになります。
    • -sオプションとは関係ありませんが、上記のexon/intronが考慮されていないことは上記の例でも確認できます。
    • アノテーション名全長exon長fastaの配列長
      MIR941-1##NR_03063711689116
      MIR1975##NR_03173937956379

実行例

  • 入力データ1:ゲノムfasta + インデックスファイル(.fai)
    • ここでは、インデックスファイルを予め用意していますが、無ければ実行時に自動生成されます。
    • index用の.faiファイルは、samtoolsで作成したものを流用可能です。
    • $ ls -1 /dataDir/genomes/hg19withSamtoolIndex
      
      hg19.fa
      
      hg19.fa.fai
      
      $ head -n 3 /dataDir/genomes/hg19withSamtoolIndex/hg19.fa
      
      >gt;chr1
      
      NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
      
      NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
  • 入力データ2:遺伝子アノテーション情報(refFLat.bed)
    • refFlat.bedは以下の2つの方法で入手可能です。方法はこちらを参照くださいgenePred2bed
      • UCSCの"refFlat.txt"をgenePred2bedで変換します。
      • UCSC Table BrowserでrefSeq GeneをBEDフォーマットでダウンロードします。
    • $ head -n 2 refFlat.bed
      
      chr19   49534925        49535044        NR_024244       0       -       49535044        49535044  0       1       119,    0,
      
      chr19   50643458        50643577        NR_024243       0       -       50643577        50643577  0       1       119,    0,
    • ゲノムの配列がchr1~chr22+chrX,chrYのみの場合は、そのままだと"chr17_ctg5_hap1"のようなハプロタイプへのアノテーションがマッピングできないのでエラーとなります。
      • unable to find FASTA index entry for 'chr17_ctg5_hap1'
    • 以下のスクリプトでハプロタイプをアノテーションを除外しておくとうまくいきます。
      • $ awk '$1~/^chr([0-9]+|X|Y)$/{print $0}' refFlat.bed > refFlat.only24chrom.bed
  • コマンド実行
    • $ toolDir/BEDTools-Version-2.12.0/bin/fastaFromBed -name -fi dataDir/genomes/hg19withSamtoolIndex/hg19.fa -bed refFlat.only24chrom.bed -fo refGenes.fa
  • 出力結果
    • トランスクリプトのfastaファイルが得られます。
    • head -n 4 refGenes.fa
      
      >NR_024244
      
      AAAAAAATTTTTTTCTCAACCAATGTGGACCCGGTTGGCCTCGAACTCGTACCCTCGAACCCTCCCCTCCCTGAGGGCCCGAGGGCACGCGCAACCGGTCGGAGCCACAATAGCTCGGG
      
      >NR_024243
      
      AAAAAAATTTTTTTTTCAACCCATGTGGACCAGGTTGGCCTTGAACTCGTACCCTCGAACCCCTGCCTCCCTGAGGGCCCGAGGGCAGGCGCATCCAGCCGGAGCCACAGTGGCTCCGG

overlapの使い方


  • BEDtoolsに同梱されていますが、xxBedという命名からも予想されるように、BEDファイルに対する処理を想定していないように見受けられます。
  • 実例としては、pair-endのマッピングを行ったBAMファイルから、insertサイズを抽出する処理を例に、ツールの動作を確認します。
  • 入力に使用するBAMファイルはSRR005720(external link)からダウンロードしたfastqをBWAのデフォルトオプションでペアエンドマッピングした結果です。
  • 入力の例はこちら
    $ samtools view SRR005720.bam | head -n 3
    
    SRR005720.1150114       163     chr1    9991    37      6M1I29M =       10407   452     GT..途中略..MD:Z:7C27
    
    SRR005720.3956281       113     chr1    10000   0       36M     chr17   81195003        0 AT..途中略..59995,36M,2;
    
    SRR005720.4118717       89      chr1    10009   0       36M     =       10009   0       AC..途中略..MD:Z:36
    • 3列目がリード1のマップ箇所の染色体番号
    • 4列目がそれぞれリード1のマップ箇所の開始座標。
    • 7列目がリード2のマップ箇所の染色体番号。リード1と同一の場合は"="。
    • 8列目ががそれぞれリード2のマップ箇所の開始座標。
    • 9列目が総断片長(insertサイズ+両リード長)。
    • 10列目が配列。
  • 前処理
    $ samtools view SRR005720.bam | awk '$7=="="{print $1"\t"$4"\t"$4+length($10)"\t"$8"\t"$8+length($10)"\t"$9}' | head
    
    -n 3
    
    SRR005720.1150114       9991    10027   10407   10443   452
    
    SRR005720.4118717       10009   10045   10009   10045   0
    
    SRR005720.4118717       10009   10045   10009   10045   0
    • 前処理としてawkコマンドを使用して、同じ染色体に乗っているものだけ対象にするため、7列目が"="のものだけ抽出する処理を挟みます。
    • 加えて、配列名、リード1開始位置、リード1開始位置+リード1長、リード2開始位置、リード2開始位置+リード1長、insert長を出力させます。
    • ここでは、リード1開始位置+リード1長をリード1終了位置の概算、リード2開始位置+リード1長をリード2終了位置の概算として使用します。
      • 正確にはInDelの有無によって、終了位置がずれる場合があるので、正しくはBEDして、ペアを対応付けるなどの対応が必要かと思われます。
  • 本処理
    $ samtools view SRR005720.bam | awk '$7=="="{print $1"\t"$4"\t"$4+length($10)"\t"$8"\t"$8+length($10)"\t"$9}' | /toolDir/BEDTools-Version-2.13.2/bin/overlap -cols 2,3,4,5 | head -n 3
    
    SRR005720.1150114       9991    10027   10407   10443   452     -380
    
    SRR005720.4118717       10009   10045   10009   10045   0       36
    
    SRR005720.4118717       10009   10045   10009   10045   0       36
    • この内、"overlap"ツールでは、2,3列、4,5列を"-cols"引数で渡しています。
    • 出力の最後に"-380","36"などの数値が表れていると思いますが、これはこのペアのリードのinsert長が10370であることを意味しています。
    • 念のため計算して確かめてみます。
      リード1開始位置リード1終了位置リード2開始位置リード2終了位置リード2開始位-リード1終了位置overlapの結果リード2終了位置-リード1開始位置総断片長
      19991100271040710443380-380452452
      210009100451000910045-3636360
    • 少なくともoverlapの結果は入力から考えると期待通りの結果になっています。

flankBedの使い方


  • flankBedは、入力となるある領域の上流域、下流域を取り出します。
    • もとの遺伝子領域は結果には含まれません。
    • 遺伝子アノテーション情報から、その上流部位だけ切りだすなどの使用方法が考えられます。
  • この例では、refSeq geneの上流部位等を切り出す方法を記載します。
    • 入力データは、UCSC Table browserから、RefGeneをBED形式でダウンロードしたものを使用しています。
    • 入力の例
      $ head -n 3 refGene.bed
      
      chr1    67051159        67163158        NM_024763       0       -       67052400        67163102  0       17      1292,157,227,99,122,158,152,87,203,195,156,140,157,113,185,175,226,       0,9472,13931,14923,20696,21102,22737,24821,27580,34595,49258,58481,61892,78265,80340,92312,111773,
      
      chr1    67075869        67163158        NM_207014       0       -       67075923        67163102  0       10      198,203,195,156,140,157,113,185,175,226,        0,2870,9885,24548,33771,37182,53555,55630,67602,87063,
      
      chr1    16762998        16812569        NM_017940       0       -       16763024        16791103  0       29      270,112,177,53,173,52,164,52,206,73,215,103,210,212,73,215,103,210,212,73,215,103,210,155,127,70,79,272,150,      0,890,1716,2490,3262,4062,5156,9225,10570,11240,12350,13400,15276,16828,17503,18627,19677,21572,23133,23808,24932,25982,27930,28242,29524,30675,31014,44591,49421,
    • 入力の例2:"-g"オプションに渡す染色体情報ファイル。
      $ head -n 3 /dataDir/hg19/chromInfo/chromInfo.txt
      
      chr1    249250621       /gbdb/hg19/hg19.2bit
      
      chr2    243199373       /gbdb/hg19/hg19.2bit
      
      chr3    198022430       /gbdb/hg19/hg19.2bit
      • この例では、3列目がありますが、不要です。1列目に染色体番号、タブで区切り、2列名に染色体の長さが入ります。この情報を用い、切りだされた領域が、染色体長を超えることを防いでいるものと考えられます。3列目以降はあっても動作上支障は無いようです。
    • 実行例1:遺伝子両端100塩基を切り出します。
      $ /toolDir/BEDTools-Version-2.13.2/bin/flankBed -i refGene.bed -g /dataDir/hg19/chromInfo/chromInfo.txt -b 100 | head -n 10
      
      chr1    67051059        67051159        NM_024763       0       -       67052400        67163102  0       17      1292,157,227,99,122,158,152,87,203,195,156,140,157,113,185,175,226,       0,9472,13931,14923,20696,21102,22737,24821,27580,34595,49258,58481,61892,78265,80340,92312,111773,
      
      chr1    67163158        67163258        NM_024763       0       -       67052400        67163102  0       17      1292,157,227,99,122,158,152,87,203,195,156,140,157,113,185,175,226,       0,9472,13931,14923,20696,21102,22737,24821,27580,34595,49258,58481,61892,78265,80340,92312,111773,
      
      chr1    67075769        67075869        NM_207014       0       -       67075923        67163102  0       10      198,203,195,156,140,157,113,185,175,226,        0,2870,9885,24548,33771,37182,53555,55630,67602,87063,
      
      chr1    67163158        67163258        NM_207014       0       -       67075923        67163102  0       10      198,203,195,156,140,157,113,185,175,226,        0,2870,9885,24548,33771,37182,53555,55630,67602,87063,
      
      chr1    16762898        16762998        NM_017940       0       -       16763024        16791103  0       29      270,112,177,53,173,52,164,52,206,73,215,103,210,212,73,215,103,210,212,73,215,103,210,155,127,70,79,272,150,      0,890,1716,2490,3262,4062,5156,9225,10570,11240,12350,13400,15276,16828,17503,18627,19677,21572,23133,23808,24932,25982,27930,28242,29524,30675,31014,44591,49421,
      
      chr1    16812569        16812669        NM_017940       0       -       16763024        16791103  0       29      270,112,177,53,173,52,164,52,206,73,215,103,210,212,73,215,103,210,212,73,215,103,210,155,127,70,79,272,150,      0,890,1716,2490,3262,4062,5156,9225,10570,11240,12350,13400,15276,16828,17503,18627,19677,21572,23133,23808,24932,25982,27930,28242,29524,30675,31014,44591,49421,
      
      chr1    50286172        50286272        NM_001144777    0       +       50286423        50439437  0       7       169,241,104,154,226,39,1018,    0,96943,129075,145751,147547,149415,152837,
      
      chr1    50440127        50440227        NM_001144777    0       +       50286423        50439437  0       7       169,241,104,154,226,39,1018,    0,96943,129075,145751,147547,149415,152837,
      
      chr1    8334950 8335050 NM_001042681    0       -       8337733 8638943 0       23      2717,181,147,721,223,1379,114,162,200,93,163,81,99,100,125,49,105,97,106,126,71,469,481,  0,3015,3696,5792,7360,7708,9359,10279,11652,12342,13408,70323,113521,142659,145001,156222,188809,204070,205013,262156,271905,303568,464755,
      
      chr1    8800286 8800386 NM_001042681    0       -       8337733 8638943 0       23      2717,181,147,721,223,1379,114,162,200,93,163,81,99,100,125,49,105,97,106,126,71,469,481,  0,3015,3696,5792,7360,7708,9359,10279,11652,12342,13408,70323,113521,142659,145001,156222,188809,204070,205013,262156,271905,303568,464755,
      • BED12特有の7列目以降は考慮されていませんので、、7列目以降は除去する方が無難です。
      • 上流の切り出しと、下流の切り出しが、別の行になるので、結果として、入力ファイルの2倍の行数が出力されます。
    • 実行例2:ストランドを考慮せずに上流100塩基、下流200塩基を切り出します。
      $ /toolDir/BEDTools-Version-2.13.2/bin/flankBed -i refGene.bed -g /dataDir/hg19/chromInfo/chromInfo.txt -l 100 -r 200 | head -n 10
      
      chr1    67051059        67051159        NM_024763       0       -       67052400        67163102  0       17      1292,157,227,99,122,158,152,87,203,195,156,140,157,113,185,175,226,       0,9472,13931,14923,20696,21102,22737,24821,27580,34595,49258,58481,61892,78265,80340,92312,111773,
      
      chr1    67163158        67163358        NM_024763       0       -       67052400        67163102  0       17      1292,157,227,99,122,158,152,87,203,195,156,140,157,113,185,175,226,       0,9472,13931,14923,20696,21102,22737,24821,27580,34595,49258,58481,61892,78265,80340,92312,111773,
      
      chr1    67075769        67075869        NM_207014       0       -       67075923        67163102  0       10      198,203,195,156,140,157,113,185,175,226,        0,2870,9885,24548,33771,37182,53555,55630,67602,87063,
      
      chr1    67163158        67163358        NM_207014       0       -       67075923        67163102  0       10      198,203,195,156,140,157,113,185,175,226,        0,2870,9885,24548,33771,37182,53555,55630,67602,87063,
      
      chr1    16762898        16762998        NM_017940       0       -       16763024        16791103  0       29      270,112,177,53,173,52,164,52,206,73,215,103,210,212,73,215,103,210,212,73,215,103,210,155,127,70,79,272,150,      0,890,1716,2490,3262,4062,5156,9225,10570,11240,12350,13400,15276,16828,17503,18627,19677,21572,23133,23808,24932,25982,27930,28242,29524,30675,31014,44591,49421,
      
      chr1    16812569        16812769        NM_017940       0       -       16763024        16791103  0       29      270,112,177,53,173,52,164,52,206,73,215,103,210,212,73,215,103,210,212,73,215,103,210,155,127,70,79,272,150,      0,890,1716,2490,3262,4062,5156,9225,10570,11240,12350,13400,15276,16828,17503,18627,19677,21572,23133,23808,24932,25982,27930,28242,29524,30675,31014,44591,49421,
      
      chr1    50286172        50286272        NM_001144777    0       +       50286423        50439437  0       7       169,241,104,154,226,39,1018,    0,96943,129075,145751,147547,149415,152837,
      
      chr1    50440127        50440327        NM_001144777    0       +       50286423        50439437  0       7       169,241,104,154,226,39,1018,    0,96943,129075,145751,147547,149415,152837,
      
      chr1    8334950 8335050 NM_001042681    0       -       8337733 8638943 0       23      2717,181,147,721,223,1379,114,162,200,93,163,81,99,100,125,49,105,97,106,126,71,469,481,  0,3015,3696,5792,7360,7708,9359,10279,11652,12342,13408,70323,113521,142659,145001,156222,188809,204070,205013,262156,271905,303568,464755,
      
      chr1    8800286 8800486 NM_001042681    0       -       8337733 8638943 0       23      2717,181,147,721,223,1379,114,162,200,93,163,81,99,100,125,49,105,97,106,126,71,469,481,  0,3015,3696,5792,7360,7708,9359,10279,11652,12342,13408,70323,113521,142659,145001,156222,188809,204070,205013,262156,271905,303568,464755,
      • ストランドを考慮しないので、逆鎖のアノテーションであっても、純鎖にあるアノテーションと同じく、ゲノム座標が小さい方が上流、大きい方が上流となっています。
    • 実行例3:ストランドを考慮して上流100塩基、下流200塩基を切り出します。
      $ /toolDir/BEDTools-Version-2.13.2/bin/flankBed -s -i refGene.bed -g /dataDir/hg19/chromInfo/chromInfo.txt -l 100 -r 200 | head -n 10
      
      chr1    67050959        67051159        NM_024763       0       -       67052400        67163102  0       17      1292,157,227,99,122,158,152,87,203,195,156,140,157,113,185,175,226,       0,9472,13931,14923,20696,21102,22737,24821,27580,34595,49258,58481,61892,78265,80340,92312,111773,
      
      chr1    67163158        67163258        NM_024763       0       -       67052400        67163102  0       17      1292,157,227,99,122,158,152,87,203,195,156,140,157,113,185,175,226,       0,9472,13931,14923,20696,21102,22737,24821,27580,34595,49258,58481,61892,78265,80340,92312,111773,
      
      chr1    67075669        67075869        NM_207014       0       -       67075923        67163102  0       10      198,203,195,156,140,157,113,185,175,226,        0,2870,9885,24548,33771,37182,53555,55630,67602,87063,
      
      chr1    67163158        67163258        NM_207014       0       -       67075923        67163102  0       10      198,203,195,156,140,157,113,185,175,226,        0,2870,9885,24548,33771,37182,53555,55630,67602,87063,
      
      chr1    16762798        16762998        NM_017940       0       -       16763024        16791103  0       29      270,112,177,53,173,52,164,52,206,73,215,103,210,212,73,215,103,210,212,73,215,103,210,155,127,70,79,272,150,      0,890,1716,2490,3262,4062,5156,9225,10570,11240,12350,13400,15276,16828,17503,18627,19677,21572,23133,23808,24932,25982,27930,28242,29524,30675,31014,44591,49421,
      
      chr1    16812569        16812669        NM_017940       0       -       16763024        16791103  0       29      270,112,177,53,173,52,164,52,206,73,215,103,210,212,73,215,103,210,212,73,215,103,210,155,127,70,79,272,150,      0,890,1716,2490,3262,4062,5156,9225,10570,11240,12350,13400,15276,16828,17503,18627,19677,21572,23133,23808,24932,25982,27930,28242,29524,30675,31014,44591,49421,
      
      chr1    50286172        50286272        NM_001144777    0       +       50286423        50439437  0       7       169,241,104,154,226,39,1018,    0,96943,129075,145751,147547,149415,152837,
      
      chr1    50440127        50440327        NM_001144777    0       +       50286423        50439437  0       7       169,241,104,154,226,39,1018,    0,96943,129075,145751,147547,149415,152837,
      
      chr1    8334850 8335050 NM_001042681    0       -       8337733 8638943 0       23      2717,181,147,721,223,1379,114,162,200,93,163,81,99,100,125,49,105,97,106,126,71,469,481,  0,3015,3696,5792,7360,7708,9359,10279,11652,12342,13408,70323,113521,142659,145001,156222,188809,204070,205013,262156,271905,303568,464755,
      
      chr1    8800286 8800386 NM_001042681    0       -       8337733 8638943 0       23      2717,181,147,721,223,1379,114,162,200,93,163,81,99,100,125,49,105,97,106,126,71,469,481,  0,3015,3696,5792,7360,7708,9359,10279,11652,12342,13408,70323,113521,142659,145001,156222,188809,204070,205013,262156,271905,303568,464755,
      • 順鎖の"+"ストランドアノテーションは、上流がゲノムの座標が小さい方向へ100塩基、下流はゲノムの座標が大きい方向へ200塩基切りだされているのに対し、
      • 逆鎖の"-"ストランドアノテーションは、上流がゲノムの座標が大きい方向へ100塩基、下流はゲノムの座標が小さい方向へ200塩基切りだされているのが分かるかと思います。
    • 実行例4:ストランドを考慮して上流100塩基のみ切り出す。
      $ /panasas/vol01/CIPpipeline/tool/BEDTools-Version-2.13.2/bin/flankBed -s -i /ftp01/inoue_group/data/DS00066004/refGene.bed -g /panasas/vol01/CIPdata/genomes/hg19/chromInfo/chromInfo.txt -l 100 -r 0 | head -n 5
      
      chr1    67163158        67163258        NM_024763       0       -       67052400        67163102  0       17      1292,157,227,99,122,158,152,87,203,195,156,140,157,113,185,175,226,       0,9472,13931,14923,20696,21102,22737,24821,27580,34595,49258,58481,61892,78265,80340,92312,111773,
      
      chr1    67163158        67163258        NM_207014       0       -       67075923        67163102  0       10      198,203,195,156,140,157,113,185,175,226,        0,2870,9885,24548,33771,37182,53555,55630,67602,87063,
      
      chr1    16812569        16812669        NM_017940       0       -       16763024        16791103  0       29      270,112,177,53,173,52,164,52,206,73,215,103,210,212,73,215,103,210,212,73,215,103,210,155,127,70,79,272,150,      0,890,1716,2490,3262,4062,5156,9225,10570,11240,12350,13400,15276,16828,17503,18627,19677,21572,23133,23808,24932,25982,27930,28242,29524,30675,31014,44591,49421,
      
      chr1    50286172        50286272        NM_001144777    0       +       50286423        50439437  0       7       169,241,104,154,226,39,1018,    0,96943,129075,145751,147547,149415,152837,
      
      chr1    8800286 8800386 NM_001042681    0       -       8337733 8638943 0       23      2717,181,147,721,223,1379,114,162,200,93,163,81,99,100,125,49,105,97,106,126,71,469,481,  0,3015,3696,5792,7360,7708,9359,10279,11652,12342,13408,70323,113521,142659,145001,156222,188809,204070,205013,262156,271905,303568,464755,
      • 上流塩基のみ切り出しているため、出力行は、入力ファイルと同数になります。

slopBedの使い方とflankBedとの違い


  • slopBedは、入力となる領域を広げ、上流域、下流域を付け加えます。
    • 操作としては、flankBedと同一ですが、もとの入力領域を含めるか否かで結果が違っています。
    • 遺伝子アノテーション情報から、その上流部位、下流部位を含めた領域を取得するのに役立ちます。
  • この例では、refSeq geneの上流部位、下流部位を付け加える例を示しています。
    • 入力データはflankBedと同じく、UCSC Table browserから、RefGeneをBED形式でダウンロードしたものを使用しています。
    • 入力の例
      $ head -n 3 refGene.bed
      
      chr1    67051159        67163158        NM_024763       0       -       67052400        67163102  0       17      1292,157,227,99,122,158,152,87,203,195,156,140,157,113,185,175,226,       0,9472,13931,14923,20696,21102,22737,24821,27580,34595,49258,58481,61892,78265,80340,92312,111773,
      
      chr1    67075869        67163158        NM_207014       0       -       67075923        67163102  0       10      198,203,195,156,140,157,113,185,175,226,        0,2870,9885,24548,33771,37182,53555,55630,67602,87063,
      
      chr1    16762998        16812569        NM_017940       0       -       16763024        16791103  0       29      270,112,177,53,173,52,164,52,206,73,215,103,210,212,73,215,103,210,212,73,215,103,210,155,127,70,79,272,150,      0,890,1716,2490,3262,4062,5156,9225,10570,11240,12350,13400,15276,16828,17503,18627,19677,21572,23133,23808,24932,25982,27930,28242,29524,30675,31014,44591,49421,
    • 実行例1:遺伝子両端100塩基を付け加えます。
      $ /toolDir/BEDTools-Version-2.13.2/bin/slopBed -i refGene.bed -g /dataDir/hg19/chromInfo/chromInfo.txt -b 100 | head -n 5
      
      chr1    67051059        67163258        NM_024763       0       -       67052400        67163102  0       17      1292,157,227,99,122,158,152,87,203,195,156,140,157,113,185,175,226,       0,9472,13931,14923,20696,21102,22737,24821,27580,34595,49258,58481,61892,78265,80340,92312,111773,
      
      chr1    67075769        67163258        NM_207014       0       -       67075923        67163102  0       10      198,203,195,156,140,157,113,185,175,226,        0,2870,9885,24548,33771,37182,53555,55630,67602,87063,
      
      chr1    16762898        16812669        NM_017940       0       -       16763024        16791103  0       29      270,112,177,53,173,52,164,52,206,73,215,103,210,212,73,215,103,210,212,73,215,103,210,155,127,70,79,272,150,      0,890,1716,2490,3262,4062,5156,9225,10570,11240,12350,13400,15276,16828,17503,18627,19677,21572,23133,23808,24932,25982,27930,28242,29524,30675,31014,44591,49421,
      
      chr1    50286172        50440227        NM_001144777    0       +       50286423        50439437  0       7       169,241,104,154,226,39,1018,    0,96943,129075,145751,147547,149415,152837,
      
      chr1    8334950 8800386 NM_001042681    0       -       8337733 8638943 0       23      2717,181,147,721,223,1379,114,162,200,93,163,81,99,100,125,49,105,97,106,126,71,469,481,  0,3015,3696,5792,7360,7708,9359,10279,11652,12342,13408,70323,113521,142659,145001,156222,188809,204070,205013,262156,271905,303568,464755,
      • BED12特有の7列目以降には影響を与えませんので、7列目以降は除去する方が無難です。
    • 実行例2:ストランドを考慮せずに上流100塩基、下流200塩基を付け加えます。
      $ /toolDir/BEDTools-Version-2.13.2/bin/slopBed -i refGene.bed -g /dataDir/hg19/chromInfo/chromInfo.txt -l 100 -r 200 | head -n 5
      
      chr1    67051059        67163358        NM_024763       0       -       67052400        67163102  0       17      1292,157,227,99,122,158,152,87,203,195,156,140,157,113,185,175,226,       0,9472,13931,14923,20696,21102,22737,24821,27580,34595,49258,58481,61892,78265,80340,92312,111773,
      
      chr1    67075769        67163358        NM_207014       0       -       67075923        67163102  0       10      198,203,195,156,140,157,113,185,175,226,        0,2870,9885,24548,33771,37182,53555,55630,67602,87063,
      
      chr1    16762898        16812769        NM_017940       0       -       16763024        16791103  0       29      270,112,177,53,173,52,164,52,206,73,215,103,210,212,73,215,103,210,212,73,215,103,210,155,127,70,79,272,150,      0,890,1716,2490,3262,4062,5156,9225,10570,11240,12350,13400,15276,16828,17503,18627,19677,21572,23133,23808,24932,25982,27930,28242,29524,30675,31014,44591,49421,
      
      chr1    50286172        50440327        NM_001144777    0       +       50286423        50439437  0       7       169,241,104,154,226,39,1018,    0,96943,129075,145751,147547,149415,152837,
      
      chr1    8334950 8800486 NM_001042681    0       -       8337733 8638943 0       23      2717,181,147,721,223,1379,114,162,200,93,163,81,99,100,125,49,105,97,106,126,71,469,481,  0,3015,3696,5792,7360,7708,9359,10279,11652,12342,13408,70323,113521,142659,145001,156222,188809,204070,205013,262156,271905,303568,464755,
      • ストランドを考慮しないので、逆鎖のアノテーションであっても、純鎖にあるアノテーションと同じく、ゲノム座標が小さい方が上流、大きい方が上流となっています。
    • 実行例3:ストランドを考慮して上流100塩基、下流200塩基を付け加えます。
      $ /toolDir/BEDTools-Version-2.13.2/bin/slopBed -s -i refGene.bed -g /dataDir/hg19/chromInfo/chromInfo.txt -l 100 -r 200 | head -n 5
      
      chr1    67050959        67163258        NM_024763       0       -       67052400        67163102  0       17      1292,157,227,99,122,158,152,87,203,195,156,140,157,113,185,175,226,       0,9472,13931,14923,20696,21102,22737,24821,27580,34595,49258,58481,61892,78265,80340,92312,111773,
      
      chr1    67075669        67163258        NM_207014       0       -       67075923        67163102  0       10      198,203,195,156,140,157,113,185,175,226,        0,2870,9885,24548,33771,37182,53555,55630,67602,87063,
      
      chr1    16762798        16812669        NM_017940       0       -       16763024        16791103  0       29      270,112,177,53,173,52,164,52,206,73,215,103,210,212,73,215,103,210,212,73,215,103,210,155,127,70,79,272,150,      0,890,1716,2490,3262,4062,5156,9225,10570,11240,12350,13400,15276,16828,17503,18627,19677,21572,23133,23808,24932,25982,27930,28242,29524,30675,31014,44591,49421,
      
      chr1    50286172        50440327        NM_001144777    0       +       50286423        50439437  0       7       169,241,104,154,226,39,1018,    0,96943,129075,145751,147547,149415,152837,
      
      chr1    8334850 8800386 NM_001042681    0       -       8337733 8638943 0       23      2717,181,147,721,223,1379,114,162,200,93,163,81,99,100,125,49,105,97,106,126,71,469,481,  0,3015,3696,5792,7360,7708,9359,10279,11652,12342,13408,70323,113521,142659,145001,156222,188809,204070,205013,262156,271905,303568,464755,
      • 順鎖の"+"ストランドアノテーションは、上流がゲノムの座標が小さい方向へ100塩基、下流はゲノムの座標が大きい方向へ200塩基付け加えられるのに対し、
      • 逆鎖の"-"ストランドアノテーションは、上流がゲノムの座標が大きい方向へ100塩基、下流はゲノムの座標が小さい方向へ200塩基付け加えられるのが分かるかと思います。
    • 実行例4:ストランドを考慮して上流100塩基のみ付け加えます。
      $ /toolDir/BEDTools-Version-2.13.2/bin/slopBed -s -i refGene.bed -g /dataDir/hg19/chromInfo/chromInfo.txt -l 100 -r 0 | head -n 5
      
      chr1    67051159        67163258        NM_024763       0       -       67052400        67163102  0       17      1292,157,227,99,122,158,152,87,203,195,156,140,157,113,185,175,226,       0,9472,13931,14923,20696,21102,22737,24821,27580,34595,49258,58481,61892,78265,80340,92312,111773,
      
      chr1    67075869        67163258        NM_207014       0       -       67075923        67163102  0       10      198,203,195,156,140,157,113,185,175,226,        0,2870,9885,24548,33771,37182,53555,55630,67602,87063,
      
      chr1    16762998        16812669        NM_017940       0       -       16763024        16791103  0       29      270,112,177,53,173,52,164,52,206,73,215,103,210,212,73,215,103,210,212,73,215,103,210,155,127,70,79,272,150,      0,890,1716,2490,3262,4062,5156,9225,10570,11240,12350,13400,15276,16828,17503,18627,19677,21572,23133,23808,24932,25982,27930,28242,29524,30675,31014,44591,49421,
      
      chr1    50286172        50440127        NM_001144777    0       +       50286423        50439437  0       7       169,241,104,154,226,39,1018,    0,96943,129075,145751,147547,149415,152837,
      
      chr1    8335050 8800386 NM_001042681    0       -       8337733 8638943 0       23      2717,181,147,721,223,1379,114,162,200,93,163,81,99,100,125,49,105,97,106,126,71,469,481,  0,3015,3696,5792,7360,7708,9359,10279,11652,12342,13408,70323,113521,142659,145001,156222,188809,204070,205013,262156,271905,303568,464755,

complementBedの使い方


  • intergenic regionのBEDをカスタムで定義したいシチュエーションなどで、全ゲノムから、遺伝子領域とその周辺領域を除外したBEDを取得する例をcomplementBedを利用して実行してみます。
  • 手順1:slopBedを利用して、遺伝子領域とその近傍前後5Kの領域を切り出します。
  • 手順2:complementBedを使用して、手順1のBEDの補集合を切り出します。結果として、遺伝子近傍5K以外の領域をintergenic regionとして切りだすことができます。
  • 用意するファイル1:遺伝子領域のbedファイル。UCSCのtable browserなどを利用し、hg19のrefFlatのbedを取得します。
  • 用意するファイル2:染色体長の情報
    • $ head -n 3 /dataDir/hg19/chromInfo/chromInfo.txt
      
      chr1    249250621       /gbdb/hg19/hg19.2bit
      
      chr2    243199373       /gbdb/hg19/hg19.2bit
      
      chr3    198022430       /gbdb/hg19/hg19.2bit
    • この例では、3列目がありますが、不要です。1列目に染色体番号、タブで区切り、2列名に染色体の長さが入ります。
  • 実行コマンド
    $ /toolDir/BEDTools-Version-2.13.2/bin/slopBed -i data/refFlat.hg19.bed -g /dataDir/hg19/chromInfo/chromInfo.txt -b 5000 | /toolDir/BEDTools-Version-2.13.2/bin/complementBed -g /dataDir/hg19/chromInfo/chromInfo.txt > intergenic.bed
    
    $ head -n 3 intergenic.bed
    
    chr1    0       6873
    
    chr1    41081   64090
    
    chr1    75008   129772
  • ちなみに、染色体番号、ゲノム開始座標順に並べた、slopBedの結果が以下のようになっているため、遺伝子と周辺領域を除外した領域になっていることが確認できるかと思います。
    • 染色体番号開始座標終了座標備考
      chr1687319409
      chr1936134370上の領域と部分重複
      chr12961041081上の領域と部分重複
      chr12961041081上の領域と部分重複
      chr16409075008
      chr1129772145566
    • $ /panasas/vol01/CIPpipeline/tool/BEDTools-Version-2.13.2/bin/slopBed -i data/refFlat.hg19.bed -g /panasas/vol01/CIPdata/genomes/hg19/chromInfo/chromInfo.txt -b 5000 | sort -k 1,1 -k 2,2n | head -n 6
      
      chr1    6873    19409   DDX11L1 0       +       14409   14409   0       3       354,109,1189,     0,739,1347,
      
      chr1    9361    34370   WASH7P  0       -       29370   29370   0       11      468,69,152,159,198,136,137,147,99,154,50, 0,608,1434,2245,2496,2871,3244,3553,3906,10376,14959,
      
      chr1    29610   41081   FAM138A 0       -       36081   36081   0       3       564,205,361,      0,666,1110,
      
      chr1    29610   41081   FAM138F 0       -       36081   36081   0       3       564,205,361,      0,666,1110,
      
      chr1    64090   75008   OR4F5   0       +       69090   70008   0       1       918,    0,
      
      chr1    129772  145566  LOC729737       0       -       140566  140566  0       3       4924,58,492,      0,5017,5302,

intersecBedの使い方


インストール


    1. ソースを本家サイトのリンク先の以下からダウンロードする。
      • バージョンなどは適宜最新のもの等選んでください。
    http://code.google.com/p/bedtools/downloads/detail?name=BEDTools.v2.10.1.tar.gz

    1. tarの中身のBEDTools-Version-2.10.1/READMEの末尾の"Installation"を参考に以下のように実行して問題なく完了。
    $ tar xvzf BEDTools.v2.10.1.tar.gz
    
    $ cd BEDTools-Version-2.10.1/
    
    $ make clean
    
    Cleaning up.
    
    $ make all
    
    $ make all
    
    Building BEDTools:
    
    =========================================================
    
    - Building in src/utils/lineFileUtilities
    
      * compiling lineFileUtilities.cpp
    
    :
    
    (中略)
    
    :
    
    - Building in src/windowBed
    
      * compiling windowMain.cpp
    
      * compiling windowBed.cpp
    
      * linking windowBed

    • ※警告は出たものの、問題なくインストール完。

参考情報


Contact us
Copyright © 2009-2017 National Institute of Genetics  [Site Policy] [Privacy Policy]