Platform for Drug Discovery


Bismark


論文


特徴


  • bowtieをマッピングのエンジンとしています。
  • C->T変換、G->A変換後のゲノムに、変換後リードをマッピングし、結果を統合します。
  • 本体がperlなので、インストールが簡単で詰まる可能性も低いように思います。
  • インデックスの作成ステップでは2.9GHzのCPUのマシンで1.5時間要しました。
  • ヒトゲノムで16Kバイトのfastaファイルに対して、インデックスのサイズが合計14Gバイト。
  • version 0.6以降でBowtie2をサポートし、defaultの出力形式も独自書式からSAMに変更。インデックスも共通では無くなっている。

出力フォーマット


1. シングルエンドの場合


  • SAM形式の標準カラムに加え、オプショナルフィールドでXX、XM、XR、XGが存在する。
  • #名称説明
    12NM-tagリファレンス(参照配列)からの編集距離
    13XX-tag塩基対塩基のミスマッチ情報
    14XM-tagメチレーション文字列
    15XR-tagアライメント時のリードの変換状態(C->T or G-A)
    16XG-tagアライメント時のゲノムの変換状態(C->T or G-A)

2. ペアエンドの場合


  • シングルエンドの場合に説明文が同じ。

3. シングルエンドの場合(0.5.4以前の標準、0.6以降では'--vanilla'オプションを指定した場合)


  • #意味凡例
    1配列ID SRR020138.15024317 SALK_2029:7:100:1672:902 length=86
    2strand(-|+)-
    3染色体 chr1
    4開始座標 57798677
    5終了座標 57798726
    6リードの配列 ATTTGGGGATGATTTTTTATTATTTTTAGGATTTATGGGATGGGAAAGAA
    7対応するゲノムの配列 ACTTGGGGATGACCTCTCATTACCCCTAGGATCTACGGGATGGGAAAGAAAC
    8メチル化状況[.zZxXhH].h..........hh.h.h....hhhh......h..z..............
    9リードの変換(CT|GA)CT
    10ゲノムの変換(CT|GA)GA
    11リードクオリティ ACCB88<=BCACCCCCCCCCCCCCCCCB6=ACCCB@)8+?@+=<@CA:@@

ペアエンドの場合(0.5.4以前の標準、0.6以降では'--vanilla'オプションを指定した場合)


  • #意味凡例
    1配列ID HWUSI-EAS611_100205:2:1:13:1732#0
    2strand(-|+)+
    3染色体 14
    4開始座標 62880539
    5終了座標 62880652
    6リードの配列1 CGGGATTTCGCGGAGTACGGGTGATCGTGTGGAATATAGA
    7対応するゲノムの配列1 CGGGACTCCGCGGAGCACGGGTGACCGTGTGGAATACAGAGT
    8メチル化状況1[.zZxXhH]Z....h.xZ.Z....h.Z......xZ..........x...
    9リードの配列2 CAACTATCTAAAACTAAAAATAACGCCGCCCAAAAACTCT
    10 対応するゲノムの配列2 TCCGGCTGTCTGGAGCTGAAGATGGCGCCGCCCAGAAGCTCT
    11メチル化状況2[.zZxXhH].zx..x...xh.h..x..h..hh.Z..Z....x..h....
    12リードの変換(CT|GA)CT
    13ゲノムの変換(CT|GA)CT
    14リードクオリティ1 IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
    15リードクオリティ2 IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

4. メチル化状況の記号の意味


  • #記号意味
    1zCpGにおける非メチル化C
    2ZCpGにおけるメチル化C
    3xCHGにおける非メチル化C
    4XCHGにおけるメチル化C
    5hCHHにおける非メチル化C
    6HCHHにおけるメチル化C


インストールログ


  • perlにパスが通っていれば、解凍するだけで(少なくともLinux環境では)インストールせずに実行可能です。
  • 試していないですが、bowtieさえ入っていればwindows環境でも動くのではないかと思います。
  • #wget でダウンロード
    $ wget http://www.bioinformatics.bbsrc.ac.uk/projects/bismark/bismark_v0.5.4.tar.gz
    
    #tarで解凍
    $ tar xvzf bismark_v0.5.4.tar.gz
    
    #ついでにテストデータもダウンロード
    $ wget http://www.bioinformatics.bbsrc.ac.uk/projects/bismark/test_data.fastq

実行ログ


1. インデックスの作成


  • この実行環境ではbowtieにパスが通っているので"which bowtie"でパスが出ます。
    • $ which bowtie
      /toolDir/bowtie-0.12.7/bowtie
      
      $ which bowtie-build
      /toolDir/bowtie-0.12.7/bowtie-build
  • パスが通っていない場合は"--path_to_bowtie <bowtieの場所>"のオプションを付加して実行する必要があると思います。
  • "./ref/"にインデックス作成対象のゲノム配列"hg19.fa"が入っています。
  • 実行後に"Bisulfite_Genome"というディレクトリが作成され、その中に"CT_conversion"と"GA_conversion"用のbowtie indexファイルがそれぞれ作成されました。
  • $ ls ref/
    
    hg19.fa
    
    $ ./bismark_v0.5.4/bismark_genome_preparation ./ref/
    
    ..長いので割愛..
    
    $ ls ref/Bisulfite_Genome/*
    
    ref/Bisulfite_Genome/CT_conversion:
    
    BS_CT.1.ebwt  BS_CT.3.ebwt  BS_CT.rev.1.ebwt  genome_mfa.CT_conversion.fa
    
    BS_CT.2.ebwt  BS_CT.4.ebwt  BS_CT.rev.2.ebwt
    
    
    
    ref/Bisulfite_Genome/GA_conversion:
    
    BS_GA.1.ebwt  BS_GA.3.ebwt  BS_GA.rev.1.ebwt  genome_mfa.GA_conversion.fa
    
    BS_GA.2.ebwt  BS_GA.4.ebwt  BS_GA.rev.2.ebwt
  • 実行時間はXeon(R)2.93GHz1コアで3時間程。メモリ利用は軽微。

2. テストラン(シングルエンド)


  • $ ./bismark_v0.5.4/bismark -q --phred33-quals -n 1 -l 50 --directional ./ref/ test_data.fastq
    
    $ ls
    
    test_data.fastq_Bismark_mapping_report.txt
    
    test_data.fastq_bismark.txt
  • テスト結果の中身

2.1. INPUT_Bismark_mapping_report.txtの中身


  • 以下のようなテキスト形式のレポートが作成されます。
  • $ cat test_data.fastq_Bismark_mapping_report.txt
    
    Bismark report for: test_data.fastq
    
    Option '--directional' specified: alignments to complementary strands will be ignored
    
    Bowtie was run against the bisulfite genome of /panasas/vol01/home/nobfujii/tmpBismark/ref/ with the specified options: -q --phred33-quals -n 1 -l 50 -k 2 --best --chunkmbs 512
    
    Final Alignment report
    
    ======================
    
    Sequences analysed in total:    10000
    
    Number of alignments with a unique best hit from the different alignments:      4795
    
    Mapping efficiency:     48.0%
    
    Sequences with no alignments under any condition:       4437
    
    Sequences did not map uniquely: 759
    
    Number of sequences with unique best (first) alignment came from the bowtie output:
    
    CT/CT:  2402    ((converted) top strand)
    
    CT/GA:  2393    ((converted) bottom strand)
    
    GA/CT:  0       (complementary to (converted) top strand)
    
    GA/GA:  0       (complementary to (converted) bottom strand)
    
    Number of alignments to (merely theoretical) complementary strands being rejected in total:       9
    
    Final Cytosine Methylation Report
    
    =================================
    
    Total number of C's analysed:   39596
    
    Total methylated C's in CpG context:     1323
    
    Total methylated C's in CHG context:    21
    
    Total methylated C's in CHH context:    100
    
    Total C to T conversions in CpG context:        665
    
    Total C to T conversions in CHG context:        9884
    
    Total C to T conversions in CHH context:        27603
    
    C methylated in CpG context:    66.5%
    
    C methylated in CHG context:    0.2%
    
    C methylated in CHH context:    0.4%

2.2. INPUT_bismark.txtの中身


  • $ head test_data.fastq_bismark.txt
    
    Bismark version: v0.5.4
    
    SRR020138.15024317 SALK_2029:7:100:1672:902 length=86   -       chr1    57798677        57798726  ATTTGGGGATGATTTTTTATTATTTTTAGGATTTATGGGATGGGAAAGAA      ACTTGGGGATGACCTCTCATTACCCCTAGGATCTACGGGATGGGAAAGAAAC      .h..........hh.h.h....hhhh......h..z..............      CT        GA      ACCB88<=BCACCCCCCCCCCCCCCCCB6=ACCCB@)8+?@+=<@CA:@@
    
    SRR020138.15024319 SALK_2029:7:100:1672:31 length=86    +       chr2    10166575        10166624  ATTTTGTTATAGAGTGGGGTATTTTCGGGAAGAAGGAGGAGGAGTGTATT      ACTCTGTTACAGAGTGGGGCATCCTCGGGAAGAAGGAGGAGGAGTGCACTGT      .h.x.....x.........h..hh.Z....................h.x.      CT        CT      BCCCCBCCCCA?:=ACCBCABCCCCCBCCA??5=9@4BB@;??B@BABBA
    
    SRR020138.15024320 SALK_2029:7:100:1672:1164 length=86  -       chr5    28344472        28344521  TATAATTGATTTTTGAATAATGTGGGGTTTAGGGGTGTTGATATTTTGTG      TACAATTAACCCTTGAACAATGTGGGGCTTAGGGGTGCTGACACTCTGTGAA      ..h......hhh.....h.........h.........x...h.h.x....      CT        GA      BACBACA:>CBBBB@@BC??CABB:=?@BC<:2::?<@A=9A@ABBB@@@
    
    SRR020138.15024321 SALK_2029:7:100:1672:433 length=86   +       chr14   38711099        38711148  TTTTGAGTAGAGAAGTTAGTATTTTAGGGAATTTTTGATTTTTTTAAGTT      CCTTGAGTAGAGAAGCCAGTATTCCAGGGAACCTCTGACCTCCTCAAGTTAC      hh.............hx......hx......hh.x...hh.hh.h.....      CT        CT      BCCBB?B@@A>@-4BBB:7@BBBCBBC@@=A@BCACA;BCBBCBB@@@BB
    
    SRR020138.15024323 SALK_2029:7:100:1672:1628 length=86  +       chr13   22989135        22989184  TAATTGAGGTATAATTAATGAATTATTGGAGATTAGATGTGGGTAAATTT      TAACTGAGGCATAATTAATGAATTACTGGAGACTAGATGTGGGCAAATCTGA      ...x.....h...............x......h..........h....x.      CT        CT      BBBCCCB@9@BCBBCCBBCCABCCBCCCB?B?CB?>=B>>>1B:?@BBBA
    
    SRR020138.15024326 SALK_2029:7:100:1672:1418 length=86  -       chr5    126218343       126218392 AGATATTAGTGATGTAAAGGTAAGTTTGAATGGGAAGATAGATATATAAA      AGACACCAGTGATGTAAAGGTAAGTTTGAATGGGAAGACAGACACACAAACA      ...h.hx...............................x...h.h.h...      CT        GA      >BBCBCCBA>CBCC;A@A30>A@AACCA:BC?3@?AB>B@CABABBBB@@
    
    SRR020138.15024328 SALK_2029:7:100:1672:110 length=86   -       chr11   123316320       123316369 TATTTTGTGGTAGAGAGGGGGTTTTATTATGTTGGTTATGTTGGTTTTAA      TATTTTGTGGTAGAGAGGGGGTTTCACCATGTTGGCCAGGCTGGTCTCAAAC      ........................h.hh.......hx...x....h.h..      CT        GA      B=BCCB:?23:=+8(;A@/3:?CBA/AA2@,=?)+:=+%)7?*-;BC=03
    
    SRR020138.15024329 SALK_2029:7:100:1672:668 length=86   -       chr1    211068652       211068701 TTATTATTTTATTTATATTGTGTGTTATATATGTGTGAAGATGTGAAAAT      TTATTATTCCATCTACATTGTGTGCCATACATGTGTGAAGATGTGAAAATGA      ........hh..h..h........hh...h....................      CT        GA      BBABBABCBB@BBB@CBCB=B@A<@B@C?B@@6>2?4>;35@:(2<@@>B
    
    SRR020138.15024331 SALK_2029:7:100:1673:1282 length=86  +       chr11   77736285        77736334  AGTTCGAGTTTTAAAAATGAAGAATTTGGAGTTTGATGTTTGAAGGTAGG      AGTCCGAGTTCCAAAAATGAAGAACTTGGAGTCTGATGTTTGAAGGCAGGAA      ...xZ.....hh............h.......x.............x...      CT        CT      B-CCACCACCBCBCBBCCBCB7BBBCCCC@ACCCBCCBACCCBB=BBB@B

3. テストラン(ペアエンド)


  • 入力データはデータサイズが小さいこともあり、ターゲットbisulfite-SeqのSRR039649のデータの先頭100エントリを抜いて使用しました。
    $ head -n 400 SRR039649_1.fastq > SRR039649_h100_1.fastq
    
    $ head -n 400 SRR039649_2.fastq > SRR039649_h100_2.fastq
  • 実行コマンドは以下
  • $ /toolDir/bismark_v0.5.4/bismark --path_to_bowtie /toolDir/bowtie-0.12.7/ /dataDir/genomes/hg19multi/bismark_v0.5.4/ -1 data/pairTest/SRR039649_h100_1.fastq -2 data/pairTest/SRR039649_h100_2.fastq
  • 出力のリスト
    $ ls -1 data/pairTest/
    
    SRR039649_h100_1.fastq
    
    SRR039649_h100_1.fastq_Bismark_paired-end_mapping_report.txt
    
    SRR039649_h100_1.fastq_bismark_pe.txt
    
    SRR039649_h100_2.fastq

3.1. INPUT_Bismark_paired-end_mapping_report.txtの中身


  • $ cat data/pairTest/SRR039649_h100_1.fastq_Bismark_paired-end_mapping_report.txt
    
    Bismark report for: data/pairTest/SRR039649_h100_1.fastq and data/pairTest/SRR039649_h100_2.fastq
    
    Bowtie was run against the bisulfite genome of /panasas/vol01/CIPdata/genomes/hg19multi/bismark_v0.5.4/ with the specified options: -q -k 2 --best --chunkmbs 512
    
    Final Alignment report
    
    ======================
    
    Sequences analysed in total:    100
    
    Number of paired-end alignments with a unique best hit: 60
    
    Mapping efficiency:     60.0%
    
    Sequences with no alignments under any condition:       31
    
    Sequences did not map uniquely: 9
    
    Number of sequences with unique best (first) alignment came from the bowtie output:
    
    CT/GA/CT:       15      ((converted) top strand)
    
    GA/CT/CT:       18      (complementary to (converted) top strand)
    
    GA/CT/GA:       10      (complementary to (converted) bottom strand)
    
    CT/GA/GA:       17      ((converted) bottom strand)
    
    Final Cytosine Methylation Report
    
    =================================
    
    Total number of C's analysed:   1314
    
    Total methylated C's in CpG context:    152
    
    Total methylated C's in CHG context:    1
    
    Total methylated C's in CHH context:    7
    
    Total C to T conversions in CpG context:        160
    
    Total C to T conversions in CHG context:        334
    
    Total C to T conversions in CHH context:        660
    
    C methylated in CpG context:    48.7%
    
    C methylated in CHG context:    0.3%
    
    C methylated in CHH context:    1.0%

3.2. INPUT_Bismark_paired-end_mapping_report.txtの中身


  • $ head -n 5 data/pairTest/SRR039649_h100_1.fastq_bismark_pe.txt
    
    Bismark version: v0.5.4
    
    SRR039649.1 HWUSI-EAS291:6:1:1:584      +       chr2    27334672        27334775        GTTGGTAAGCGTCGCGCGTAGTTTCGGTNNNNNNNC    GCCGGCAAGCGCCGCGCGCAGCCCCGGTTCGCGCTCGC        .xz..h...Z.xZ.Z.Z.x..hhxZ..........Z    CCTTTCTATACTAACCAAAACCGAAACCAACTCCTT    CCCCTTTCTGTACTGGCCAAGGCCGGAGCCGGCTCCTT        .......x....xh....hh..Zx.h..zx......    CT      CT      BCCCCBCBBB<?<*@BBBABBA%%%%%%!!!!!!!%    BCCBCBBBBAABBCBBCBBCBCABCCBBBCBCCCBA
    
    SRR039649.2 HWUSI-EAS291:6:1:1:896      +       chr20   44685815        44685927        GTTCAATGTGCGGCGTATGTATATGGTCNNNNNNNT    GTCCAACGTGCGGCGCATGCACACGGCCGTGCGGCTGA        ..hH..z...Z..Z.h...h.h.z..xZ........    TACAATTACGAAAAAACCCAAACATATTAAACAAAA    ATTGCGGTTGCGGGGAGGCCCAGGCATGTTGAGCAAAA        .h.zx..h.Zxhh.hh....xh...h..h.h.....    CT      CT      BCCCCCCBACB77BCBCCCBCCBC%%%%!!!!!!!%    BCBABA3A@>?B??BBABBBBBBCBBCCBBBBBBA<
    
    SRR039649.3 HWUSI-EAS291:6:1:1:1082     -       chr20   34204064        34204163        GTTTTTGTTGAGAGTTTTTTGTTTATGTNNNNNNNT    GCTTCTGCTGAGAGCTTCCTGCCCATGTGGTTCCTCCG        .h..x..x......h..hx..hhh...........x    CCCAAAAACAAAAAAACCCAATACCCAAACTACAAT    AGCCCGAGGGCAGAAGAGCCCGGGGCCCGAGCTGCGGT        ...z.hhh..x..h.h...zx.h...z.h..x.zx.    CT      GA      BBCCCCBBCCB<=1=CCCCBAACB%%%%!!!!!!!%    BCCBCBBBBBBBABABBBCCBCCCCCB?ACBCCCCC
    
    SRR039649.4 HWUSI-EAS291:6:1:1:1425     -       chr13   113973375       113973509       GTATCGTTTATTATATTTTAAGTTTCGTNNNNNNNT    GCACCGCCCACCACACCCCAAGCCCCGTGGCTGGACAG        .h.xZ.hhh.hh.h.hhhh...hhxZ.........x    AACCGACCACGCAACACAACTTCCTACGTCCACTCC    GTGGCCGGCCACGCAGCACAGCTTCCTGCGTCCACTCC        hh..Zx....Z..x....x......x.Z........    CT      GA      BBCCCCBCCBCCBBABCCCBBB@B%%%%!!!!!!!%    BCCCBBBBBBBBBBCCBBCBBBCCBCCCBCBBBBBC

4. 出力のsam変換(シングルエンド)


  • sam変換には、header付加のために染色体の長さ情報ファイルが必要なようです。以下のような書式で用意するか、bamのヘッダにある情報を抜き出す等すると良いかと思います。
  • $ head /dataDir/genomes/hg19/chromInfo/chromInfo.txt 
    
    chr1    249250621
    chr2    243199373
    chr3    198022430
    chr4    191154276
    chr5    180915260
    chr6    171115067
    chr7    159138663
    chrX    155270560
    chr8    146364022
    chr9    141213431
    :
  • 実行は以下のようにすれば動作が確認できました。
  • $ ./bin/bismark2SAM_v5_xm.pl -c  /dataDir/genomes/hg19/chromInfo/chromInfo.txt -i data/singleTest/test_data.fastq_bismark.txt
  • 出力例は以下のようになります。
    • "@"で始まるヘッダ行が邪魔なので、"grep -v "@""で除外しています。
  • $ grep -v "@" data/singleTest/test_data.fastq_bismark.txt.sam | head  -n 5
    
    SRR020138.15024379 SALK_2029:7:100:1673:1570 length=86  16      chr4    149800847       255     50M     *       0       0       ATCTCATCTAAAAAATAACAAATCCTTTCTTCCTATAAACCTATAAAATC    BBCCCCBCCBBABBBCBAA;B8BB?=/)?CCCB?BBBBBCC=CBBBBB:B      NM:i:10 MD:Z:A8A1A4AA3A14A1A3A4A2     XM:Z:..h....x...h.h..............h...hh....h.x........h
    
    SRR020138.15024520 SALK_2029:7:100:1675:1949 length=86  16      chr4    162629397       255     50M     *       0       0       CAAATCTCTTTACTTTTAAACTTTTATAACATAAATTATTCCGCTACCCT    BBBC:>CCCC>BCBCBCCBCBCCBCBCCCBACCAABBB?ABCBBB?CBCC      NM:i:3  MD:Z:2A16A25A4  XM:Z:....x..Z......................h................x..
    
    SRR020138.15024544 SALK_2029:7:100:1675:484 length=86   0       chr13   75419792        255     50M     *       0       0       ATTATTATATTGAAAATAGATTGAGTTTATAGGTTTTTTAGTTTTGATTG    BCCCCCBCBCB:AB?BCB:;BC=;:ACCBB;49=ACBBA9+>CBBBBBB=      NM:i:9  MD:Z:2T1T2T8T3T15T1T3TT6      XM:Z:..h.h..h........x...x...............h.x...hx......
    
    SRR020138.15024569 SALK_2029:7:100:1676:419 length=86   16      chr8    23164539        255     50M     *       0       0       GAAATAACGTCAATTAAACTAAAAATAACATTTTTTCCCATTCCCACTAA    BCB?<;4=;AB1*/<AB?>ABABCBCCCCBB<BCB:<A?.7>?>AABBBB      NM:i:7  MD:Z:2AA2A13AA4AA22   XM:Z:......................hh....hx...........Z.h..hh.Z
    
    SRR020138.15024647 SALK_2029:7:100:1677:612 length=86   16      chr3    152504928       255     50M     *       0       0       AAATATAAACCCATATATCAAAATTACAATATTTTCTTAAAACATATTTT    BCCBCCCCBCCCAA?BCBACBCB?BBBCCCB;>B?B<B112<CBABBCCC      NM:i:12 MD:Z:A5AAA5A1A3A7A1A7AA1A8    XM:Z:........h.hh.......h.x.......x...h.h.....hhh.....x

5. 出力のsam変換(ペアエンド)


  • 実行は以下のようにすれば動作が確認できました。
  • $ ./bin/bismark2SAM_v5_xm.pl -p -c /dataDir/genomes/hg19/chromInfo/chromInfo.txt -i data/pairTest/SRR039649_h100_1.fastq_bismark_pe.txt  -o data/pairTest/SRR039649_h100_1.fastq_bismark_pe.txt.sam
  • 出力例は以下のようになります。
    • "@"で始まるヘッダ行が邪魔なので、"grep -v "@""で除外しています。
  • $ grep -v "@" data/pairTest/SRR039649_h100_1.fastq_bismark_pe.txt.sam | head -n 5
    
    SRR039649.1 HWUSI-EAS291:6:1:1:582      131     chr2    27334740        255     36M     chr2      27334672        -103    AAGGAGTTGGTTTCGGTTTTGGTTAGTATAGAAAGG    BCCBCBBBBAABBCBBCBBCBCABCCBBBCBCCCBA      NM:i:9  MD:Z:6TT2T1T3TT4TT4T7   XM:Z:.......x....xh....hh..Zx.h..zx......
    
    SRR039649.2 HWUSI-EAS291:6:1:1:896      67      chr20   44685815        255     36M     chr20     44685892        112     GTTCAATGTGCGGCGTATGTATATGGTCNNNNNNNT    BCCCCCCBACB77BCBCCCBCCBC%%%%!!!!!!!%      NM:i:14 MD:Z:2T3T8T3T1T1T2T1NNNNNNN1    XM:Z:..hH..z...Z..Z.h...h.h.z..xZ........
    
    SRR039649.3 HWUSI-EAS291:6:1:1:1082     115     chr20   34204128        255     36M     chr20     34204064        -99     ANNNNNNNACATAAACAAAAAACTCTCAACAAAAAC    BBCCCCBBCCB<=1=CCCCBAACB%%%%!!!!!!!%      NM:i:17 MD:Z:ANNNNNNN4AAA2AA2A6A2A2A1   XM:Z:.h..x..x......h..hx..hhh...........x
    
    SRR039649.3 HWUSI-EAS291:6:1:1:1082     179     chr20   34204064        255     36M     chr20     34204128        99      CCCAAAAACAAAAAAACCCAATACCCAAACTACAAT    BCCBCBBBBBBBABABBBCCBCCCCCB?ACBCCCCC      NM:i:16 MD:Z:3A1AAA2A2A1A3AATA3A1A2A1AA1        XM:Z:...z.hhh..x..h.h...zx.h...z.h..x.zx.
    
    SRR039649.4 HWUSI-EAS291:6:1:1:1422     179     chr13   113973375       255     36M     chr13     113973474       134     AACCGACCACGCAACACAACTTCCTACGTCCACTCC    BCCCBBBBBBBBBBCCBBCBBBCCBCCCBCBBBBBC      NM:i:6  MD:Z:AA3A7A4A6A10       XM:Z:hh..Zx....Z..x....x......x.Z........

sam出力からメチル化情報を取り出す際の注意


  • 以下のように、paired-endの出力の"XM:Z:"タグに含まれるメチル化情報は、逆鎖をとる必要があったり無かったりするので、使用時には注意が必要です。
    • 旧バージョンのバージョン依存のバグのようです。少なくともSAM出力default対応したv.0.7.7では問題が解決しています。

比較的新しい情報 v0.7.7


  • 公共データSRR094905(ペアエンド)をbismark 0.7.7でdefaultのSAM形式でマップした出力例。
    
    #first strandが順鎖、second strandが逆鎖にヒットした例
    
    SRR094905.487_HWI-BRUNOP20X:0617:2:1:10579:1235_length=75/1     67      chr15   39211433  255     75M     =       39211541        183     TTTCGGGTTTAAGTGATTTTTTAGTTTNNNNNTTTTAAGTAGNNGGANNNNTAGGGGTTTATTATTATGATTAGG       HHHHHHHHHHHHHHHEHHHHHHFHF54!!!!!4544455555!!555!!!!455544354@<F>GFF>F/F?BFC       NM:i:30 XX:Z:1CC6C8C1CC2CC1CAGCC1CCC6CT3ACTA6CCC1CC1CC5C3 XM:Z:.hxZ.....h........h.hx..hh.......hhh.....................hhh.hh.hh.....x...        XR:Z:CT   XG:Z:CT
    
    SRR094905.487_HWI-BRUNOP20X:0617:2:1:10579:1235_length=75/2     131     chr15   39211541  255     75M     =       39211433        -183    GTTTTATTATATTGGTTAGGTTGGTTTCGAATTTTTGATTTTGTGATTTGTTTATTTTATTTTTTTAAAGTGTTG       GEHHHHHGHHFGHHHHGHHHHGHHGHFHHHGEHHHHHHHHHHHHHGGHHGBGHGGGGGGGHHHHHHHHHHHHHHH       NM:i:26 XX:Z:4C1C2C5CC3C4C5C1CC3CC7C2CCC1CC1C1CCC1CCC6C2  XM:Z:....h.h..h.....hx...x....h.Z...h.hx...hh.......x..hhh.hh.h.hhh.hhh......x..        XR:Z:GA   XG:Z:CT
    
    #上記の例から、順にメチル化情報(XMタグ)、リード配列、ゲノム配列を順に抜き出したもの。全てゲノム配列に対する順鎖情報に対応しているようです。
    
    First strand methylation state_: .hxZ.....h........h.hx..hh.......hhh.....................hhh.hh.hh.....x...
    First strand read sequence_____: TTTCGGGTTTAAGTGATTTTTTAGTTTNNNNNTTTTAAGTAGNNGGANNNNTAGGGGTTTATTATTATGATTAGG
    First strand genome sequence___: TCCCGGGTTCAAGTGATTCTCCAGCCTCAGCCTCCCAAGTAGCTGGAACTATAGGGGCCCACCACCATGATCAGG
    
    Second strand methylation state: ....h.h..h.....hx...x....h.Z...h.hx...hh.......x..hhh.hh.h.hhh.hhh......x..
    Second strand read sequence____: GTTTTATTATATTGGTTAGGTTGGTTTCGAATTTTTGATTTTGTGATTTGTTTATTTTATTTTTTTAAAGTGTTG
    Second strand genome sequence__: GTTTCACTACATTGGCCAGGCTGGTCTCGAACTCCTGACCTTGTGATCTGCCCACCTCACCCTCCCAAAGTGCTG
    
    #first strandが逆鎖、second strandが順鎖にヒットした例
    
    SRR094905.499_HWI-BRUNOP20X:0617:2:1:12572:1236_length=75/1     115     chr1    212584290 255     75M     =       212584205       -160    ACCATCTTTACAAAATATACTCTCNNNNTATNNTCACATAATCNNNNNCTTAAAAAACACCAACCATATTAAATT       :AA=9F?FFB==@=<+45944454!!!!455!!4434455554!!!!!55GGGHHHHHHHHHHHHHHHHHHHHHH       NM:i:26 XX:Z:G12GG1G1G5TCTG1G1CT6GG2ATCGT3G2GG6G4G2GG3  XM:Z:x............hh.h.h..........h.........hh..........h..hh......x....h..hh...  XR:Z:CT XG:Z:GA
    
    SRR094905.499_HWI-BRUNOP20X:0617:2:1:12572:1236_length=75/2     179     chr1    212584205 255     75M     =       212584290       160     CATACCTCTCTCCTAACTTCTAATAACAACCAACAATCCTTAAAATTCGTTAACTCATAAAAACATCCAATCTCC       HHHHHHHHHHHHHHHHHHHHHHHGHHHHIHHHHHHHIHHGHHHHHEHHHHEHHHHHGHHFHHCEG;HHHHGHHHH       NM:i:17 XX:Z:3G10GG5GG1G3G3G8GG1G6GG3G2G1GG12   XM:Z:...h..........xh.....xh.h...x...x........hh.h...Z..hh...z..h.hh............  XR:Z:GA XG:Z:GA
    
    #上記の例から、順にメチル化情報(XMタグ)、リード配列、ゲノム配列を順に抜き出したもの。逆鎖の場合も全てゲノム配列に対する順鎖情報に対応しているようです。
    
    First strand methylation state_: x............hh.h.h..........h.........hh..........h..hh......x....h..hh...
    First strand read sequence_____: ACCATCTTTACAAAATATACTCTCNNNNTATNNTCACATAATCNNNNNCTTAAAAAACACCAACCATATTAAATT
    First strand genome sequence___: GCCATCTTTACAAGGTGTGCTCTCTCTGTGTCTTCACATGGTCATCGTCTTGAAGGACACCAGCCATGTTGGATT
    
    Second strand methylation state: ...h..........xh.....xh.h...x...x........hh.h...Z..hh...z..h.hh............
    Second strand read sequence____: CATACCTCTCTCCTAACTTCTAATAACAACCAACAATCCTTAAAATTCGTTAACTCATAAAAACATCCAATCTCC
    Second strand genome sequence__: CATGCCTCTCTCCTGGCTTCTGGTGACAGCCAGCAATCCTTGGAGTTCGTTGGCTCGTAGAGGCATCCAATCTCC
    
    #insertionのある例 順鎖
    
    SRR094905.1352_HWI-BRUNOP20X:0617:2:1:1707:1306_length=75/1     67      chr8    46980085        255     44M1I30M        =       46980236        226     TTTATAGAATTAAGTTTTATTTTTGATTTAGTTGTTTGTAAATATTTTTTTTTTGTGGAATTTGGAAAGGGATGT     HHHHHIHHHHHHHHHHHHHHGGGGGHHHHEFHHFHHHHHHHHHHHHFHGGGGGG#####################     NM:i:14 XX:Z:C1C1C9CC3C8C9C3C1NC15C10C2 XM:Z:h.h.x.........hh...h........x.........h...h..h...............x..........z..        XR:Z:CT XG:Z:CT
    
    Forward with insertion methylation state_: h.h.x.........hh...h........x.........h...h..h...............x..........z..
    Forward with insertion read sequence_____: TTTATAGAATTAAGTTTTATTTTTGATTTAGTTGTTTGTAAATATTTTTTTTTTGTGGAATTTGGAAAGGGATGT
    Forward with insertion genome sequence___: CTCACAGAATTAAGCCTTACTTTTGATTCAGTTGTTTGCAAACA-CTTTTTTTTGTGGAATCTGGAAAGGGACGTT
    Forward with insertion CIGAR_____________: <====================44====================>I<=============30=============>
    Forward with insertion position__________: 12345678901234567890123456789012345678901234 123456789012345678901234567890
    
    #insertionのある例 逆鎖
    
    SRR094905.2599_HWI-BRUNOP20X:0617:2:1:3844:1407_length=75/1     115     chr3    30878631        255     26M3I46M        =       30878523        -180    TATTTTTAAAAAAAAATAAACATCAATTTTTTAATATTTAAAAAAAACTTTAAAATACTCTAAAAAATTTCAAAC     BHGEGGEEEEEEEEHHGHFHHGFHHHFHGHHHHHFHHHCEEHHHHHHHHHHHHHHHHFHFHHHHHHHHHHHHHHH     NM:i:14 XX:Z:12G1G10GNNN10G3G2G15G2G9   XM:Z:............h.h..........x.............h...h..h...............h..h.........        XR:Z:CT XG:Z:GA
    
    Reverse with insertion methylation state_: ............h.h..........x.............h...h..h...............h..h.........
    Reverse with insertion read sequence_____: TATTTTTAAAAAAAAATAAACATCAATTTTTTAATATTTAAAAAAAACTTTAAAATACTCTAAAAAATTTCAAAC
    Reverse with insertion genome sequence___: TATTTTTAAAAAGAGATAAACATCAG---TTTAATATTTGAAAGAAGCTTTAAAATACTCTAGAAGATTTCAAACAGG
    Reverse with insertion CIGAR_____________: <============26==========>III<====================46======================>
    Reverse with insertion position__________: 12345678901234567890123456   1234567890123456789012345678901234567890123456
    
    
    #deletionのある例 順鎖
    
    SRR094905.2185_HWI-BRUNOP20X:0617:2:1:4781:1374_length=75/1     67      chr6    3598762 255     36M1D39M        =       3598829 142     ATTTAAATTGTTTTTGGTTTTATGTAAATTTTAGGATTTTTTTTTTATTTTTATAAAAAATGATATTGGAATTTT     HHHHGHHHHHHHHEEGFHHHHHHHHHHDHHHDHHHEHHHHHHGGGGGGGGGGGEGHHHG;CH?H=@#########     NM:i:10 XX:Z:7C3C7CC3C19C5C12C5T5       XM:Z:.......x...h.......hh...h...................h.....h............h...........        XR:Z:CT XG:Z:CT
    
    Forward with deletion methylation state_: .......x...h.......hh...h...........-........h.....h............h...........
    Forward with deletion read sequence_____: ATTTAAATTGTTTTTGGTTTTATGTAAATTTTAGGA-TTTTTTTTTTATTTTTATAAAAAATGATATTGGAATTTT
    Forward with deletion genome sequence___: ATTTAAACTGTCTTTGGTTCCATGCAAATTTTAGGATTTTTTTTTCTATTTCTATAAAAAATGACATTGGTATTTT
    Forward with deletion CIGAR_____________: <================36================>D<===================39================>
    Forward with deletion position__________: 123456789012345678901234567890123456 123456789012345678901234567890123456789
    
    
    #deletionのある例 逆鎖
    
    SRR094905.1216_HWI-BRUNOP20X:0617:2:1:7029:1295_length=75/2     179     chr2    113452479       255     39M3D36M        =       113452617       213     TAACTAAAATTACAAACACCTACCACCATACCAAAATAATTTTTTTTTTTTTTTTAAAACGAAATCTCGCTCTAT     HHHHHHHHHHHHHHHHHHHHHHHHHGHHHHHHHHFHHDGEHHHFHEGEEGEEBCEEDEGEBEFF###########     NM:i:16 XX:Z:2G2GG7G14G3GG2G17G1G3G1G9G1        XM:Z:..h..xh.......x..............h...xh..h.................h.h..Zx.h....Z....x.        XR:Z:GA XG:Z:GA
    
    
    Reverse with deletion methylation state_: ..h..xh.......x..............h...xh..h.---................h.h..Zx.h....Z....x.
    Reverse with deletion read sequence_____: TAACTAAAATTACAAACACCTACCACCATACCAAAATAA---TTTTTTTTTTTTTTTTAAAACGAAATCTCGCTCTAT
    Reverse with deletion genome sequence___: TAGCTGGAATTACAGACACCTACCACCATGCCAGGATGATTTTTTTTTTTTTTTTTTTGAGACGGAGTCTCGCTCTGT
    Reverse with deletion CIGAR_____________: <===================39================>DDD<================36================>
    Reverse with deletion position__________: 123456789012345678901234567890123456789   123456789012345678901234567890123456
  • 余談:ゲノム配列を切りだす方法
    $ chr=chr1;pos=39211433;len=75;samtools faidx hg19.fa $chr:$pos-`expr $pos + $len - 1` | tr a-z A-Z
    
    >CHR1:39211433-39211507
    
    CGAGGTATCTGAGCCATTCTAGCCTGCTGGCTCACATTGGCCCAAGGAACAGCAGTTGTA
    
    ATAGGCATTTTCAGC
    
    $
  • Nのある特殊な例
    SRR094905.3584471_HWI-BRUNOP20X:0617:2:2:3835:82177_length=75/2 131     chr2  92326102        255     4M1I70M =       92325999        -177    GTTTTTGAGGATTTTGTTGGAAATGGGATTGTTTTTATATAAATTTTAGATAGAAGTATTTTTAGAAGTTTTATT     DBBCGHGGGH1144-HIGHHHHHHHGHGGGGBHHHHHHHHHGHHHBHHHHHGEHHGGHHHFHIIHD<<,@GHHH@     NM:i:18 XX:Z:1C2N9C8C8C2C7C1C4C5C3C1C5C2NNNN    XM:Z:.h............z........z........h..h.......h.h....x.....h...h.x.....h......        XR:Z:GA XG:Z:CT
    
    4M1I70M
    
    1C2N9C8C8C2C7C1C4C5C3C1C5C2NNNN
    
     C  N         C        C        C  C       C C    C     C   C C     C  NNNN
    
    GTTTTTGAGGATTTTGTTGGAAATGGGATTGTTTTTATATAAATTTTAGATAGAAGTATTTTTAGAAGTTTTATT
    
    .h............z........z........h..h.......h.h....x.....h...h.x.....h......
    
    
    
    GCTT TGAGGATTTCGTTGGAAACGGGATTGTCTTCATATAAACTCTAGACAGAAGCATTCTCAGAAGCTTNNNNNN
    
    <4=>I<==================================70=================================>
    
    SRR094905.5940808_HWI-BRUNOP20X:0617:2:3:1189:62849_length=75/2 131     chr14 80576475        255     5M2I68M =       80576385        -163    NNNGTNGAAAGTTTGATTTATTAGATATGGTAAATGTTTAGTTGTGGAGGAAATTAAATTGAAATAAGATGGGTT     !!!##!####GHGHGGFGGDG@GDGDFD@EHHFFHD@>A8'42+4<+0AGBCC=EADAFDGDGGH<CHGGHD=GB     NM:i:14 XX:Z:ACT1C1N23C6C3C11CC3C15C    XM:Z:....h.........................h......h...x...........hh...x...............z        XR:Z:GA XG:Z:CT
    
    5M2I68M
    
    ACT1C1N23C6C3C11CC3C15C
    
    ACT C N                       C      C   C           CC   C               C
    
    NNNGTNGAAAGTTTGATTTATTAGATATGGTAAATGTTTAGTTGTGGAGGAAATTAAATTGAAATAAGATGGGTT
    
    ....h.........................h......h...x...........hh...x...............z
    
    ACTGC  AAAGTTTGATTTATTAGATATGGCAAATGTCTAGCTGTGGAGGAAACCAAACTGAAATAAGATGGGTCGCA
    
    <=5=>II<==============================68==================================>
    
    12345  12345678901234567890123456789012345678901234567890123456789012345678

古い情報 v0.5.4


  • 以下v0.5.4の際に記載した過去仕様/バグです。
  • 結論から言えば、ゲノム上の後に来るリード(順方向のread2,逆方向のread1)の"XM:Z:"はreverseする必要があります。
    • versionはbismark_v0.5.4とbismark2SAM_v5_xm.plです。
  • 順鎖のヒット
    
    SRR039649.1 HWUSI-EAS291:6:1:1:584      +       chr2    27334672        27334775        GTTGGTAAGCGTCGCGCGTAGTTTCGGTNNNNNNNC    GCCGGCAAGCGCCGCGCGCAGCCCCGGTTCGCGCTCGC  .xz..h...Z.xZ.Z.Z.x..hhxZ..........Z    CCTTTCTATACTAACCAAAACCGAAACCAACTCCTT    CCCCTTTCTGTACTGGCCAAGGCCGGAGCCGGCTCCTT  .......x....xh....hh..Zx.h..zx......    CT      CT      BCCCCBCBBB<?<*@BBBABBA%%%%%%!!!!!!!%    BCCBCBBBBAABBCBBCBBCBCABCCBBBCBCCCBA
    
    当該領域のゲノム切り出し(hg19)
    
    GCCGGCAAGCGCCGCGCGCAGCCCCGGTTCGCGCTCGCGCCTCTCCCTACACACCTCCCTGCAAGCTGAAGGAGCCGGCTCCGGCCTTGGCCAGTACAGAAAGG
    
    read1 genome
    
    GCCGGCAAGCGCCGCGCGCAGCCCCGGTTCGCGCTCGC
    
    read1 orig read
    
    GTTGGTAAGCGTCGCGCGTAGTTTCGGTNNNNNNNC
    
    read1 methyl
    
    .xz..h...Z.xZ.Z.Z.x..hhxZ..........Z
    
    read2 genome complement
    
    ....................................................................AAGGAGCCGGCTCCGGCCTTGGCCAGTACAGAAAGGGG
    
    read2 orig read complement
    
    ....................................................................AAGGAGTTGGTTTCGGTTTTGGTTAGTATAGAAAGG
    
    read2 methyl complement
    
    ____________________________________________________________________......xz..h.xZ..hh....hx....x.......
    
    sam出力
    
    SRR039649.1 HWUSI-EAS291:6:1:1:584      67      chr2    27334672        255     36M     chr2    27334740        103     GTTGGTAAGCGTCGCGCGTAGTTTCGGTNNNNNNNC    BCCCCBCBBB<?<*@BBBABBA%%%%%%!!!!!!!%    NM:i:15 MD:Z:1TT2T5T6T2TTT4NNNNNNN1     XM:Z:.xz..h...Z.xZ.Z.Z.x..hhxZ..........Z
    
    SRR039649.1 HWUSI-EAS291:6:1:1:582      131     chr2    27334740        255     36M     chr2    27334672        -103    AAGGAGTTGGTTTCGGTTTTGGTTAGTATAGAAAGG    BCCBCBBBBAABBCBBCBBCBCABCCBBBCBCCCBA    NM:i:9  MD:Z:6TT2T1T3TT4TT4T7   XM:Z:.......x....xh....hh..Zx.h..zx......
    
    逆鎖のヒット
    
    SRR039649.8 HWUSI-EAS291:6:1:1:1965     -       chr20   13765483        13765598        ACCATTACAACAAAACATCACACAAAAANNNNNNAA    TTGCCATTACGGCGGAGCGTCACGCGGAGGATTGTGGG  h.......zx.zx.h.z....z.zx.hh......hh    TTTTTTATTTTGAGTTTTTTAATTTTAGTTTATTTT    CCCTTCATCTCGAGCTCTTCAACTCCAGTCCATCCCCA  hhh..h..h.z...h.h..h..h.hx...hh..hhh    GA      CT      BCCCBBCCBBBBBB?BCBCCCCCBA%%%!!!!!!%%    BCCBCBABBBBB@C@BBBBBB>CBBB9B@BBBCBBC
    
    当該領域のゲノム切り出し(hg19)
    
    CCCTTCATCTCGAGCTCTTCAACTCCAGTCCATCCCCACTCACCGTCCGCAGTCCTACCAAGCCTCACGTGGGGCTCACACCCACAATCCTCCGCGTGACGCTCCGCCGTAATGGC
    
    read2 genome
    
    CCCTTCATCTCGAGCTCTTCAACTCCAGTCCATCCCCA
    
    read2 orig read
    
    TTTTTTATTTTGAGTTTTTTAATTTTAGTTTATTTT
    
    read2 methyl
    
    hhh..h..h.z...h.h..h..h.hx...hh..hhh
    
    read1 genome complement
    
    ................................................................................CCCACAATCCTCCGCGTGACGCTCCGCCGTAATGGCAA
    
    read1 orig read complement
    
    ................................................................................TTNNNNNNTTTTTGTGTGATGTTTTGTTGTAATGGT
    
    read1 methyl complement
    
    ________________________________________________________________________________hh......hh.xz.z....z.h.xz.xz.......h
    
    sam出力
    
    SRR039649.8 HWUSI-EAS291:6:1:1:1965     115     chr20   13765563        255     36M     chr20   13765483        -115    TTNNNNNNTTTTTGTGTGATGTTTTGTTGTAATGGT    BCCCBBCCBBBBBB?BCBCCCCCBA%%%!!!!!!%%    NM:i:20 MD:Z:TTNNNNNNTT1TT1T4T1T1TT1TT7T        XM:Z:h.......zx.zx.h.z....z.zx.hh......hh
    
    SRR039649.8 HWUSI-EAS291:6:1:1:1962     179     chr20   13765483        255     36M     chr20   13765563        115     TTTTTTATTTTGAGTTTTTTAATTTTAGTTTATTTT    BCCBCBABBBBB@C@BBBBBB>CBBB9B@BBBCBBC    NM:i:17 MD:Z:TTT2T2T1T3T1T2T2T1TT3TT2TTT        XM:Z:hhh..h..h.z...h.h..h..h.hx...hh..hhh

参考


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