BWAの出力するSAMのミスマッチ表記

BWAの出力するSAMのミスマッチ表記について

• ミスマッチやInDelはSAM形式の出力の場合以下のようになります。
  • 厳密に言えば、SAM形式ではInDelに関係する列=フィールドは存在するものの、ミスマッチ情報は表記対象に含まれておらず、BWAの拡張(オプショナルなフィールド)にミスマッチ情報が含まれます。
  • 6列目のCIGARフィールドにInDel情報が含まれます。これはSAMの書式に含まれます。
  • 19列目前後に含まれるオプショナルなフィールドにミスマッチ情報を格納するMDタグが存在します。
  • SAM形式では、CIGARとMDタグの両方を考慮しないとゲノムとのテキストアライメントを作れないので、自前でテキストアライメントするSAMパーサを作成する場合には骨が折れます。
  • 経験的には、以下の3パターンのような特異な出力行を適切にパースできれば、パーサは他の入力でもうまく動きます。

BWAの出力例とそれに対応するテキストアライメント

1. 出力例1
   • 　InsertとDeleteを両方含む場合。
   1. BWAの出力文字列

      • HWUSI-EAS733_0001:1:114:1538:876#0      0       chr2    56022580        37      18M2D6M1I11M    *       0       0       AAAAAAAAAAAAAAAAAAATTTTTAAAATTTTTAAA    ffacc^bbbbggfggcgg]feb`e^RIIITeea[eB    XT:A:U  NM:i:3  X0:i:1  X1:i:0  XM:i:1  XO:i:1  XG:i:1  MD:Z:18^AG17

   2. 出力例中のミスマッチ等を表す項目1
      • 6列目のCIGAR列の文字列 "18M2D6M1I11M"
   3. 出力例中のミスマッチ等を表す項目2
      • 19列目のMDタグの文字列 "MD:Z:18^AG17"
   4. ゲノムの配列の取得
      • ゲノムにどのように当たっているかを確認するために、ゲノムの特定部位を切り出します。ゲノムのFASTA(1行50塩基の場合)を使用して簡易的に染色体位置を表示させるコマンド

      • cat -n /panasas/vol01/CIPdata/genomes/hg19/fasta/chr2.fa | awk '{print ($1-2)*50+1 "\t" $2}' | less
        
        56022501        gcagtgagccgagattgcgcccctgcagtccgcagtccggcctgggcgac
        56022551        agagcgagactccgtctcaaaaagaaaaaaaaaaaaaaaaaaaaaaaaGA
        56022601        TTTTTAaatttttaaatttttgcagagatgaaggttttactatgttaccc

   5. テキストアライメント結果
      • ゲノムとタグの情報を手作業でアライメントしたもの

      • An image of meaning of CIGAR                 <-------18------->DD<-6-->I<---11---->
        An image of meaning of MD tag                <-------18------->AG<-------17------->
        read                                         AAAAAAAAAAAAAAAAAA--ATTTTTAAAATTTTTAAA
        alignment                                    ||||||||||||||||||  |||||| |||||||||||
        56022551        agagcgagactccgtctcaaaaagaaaaaaaaaaaaaaaaaaaaaaaaGATTTTT-Aaatttttaaatttttgcagagatgaaggttttactatgttaccc
        last 1 digit    1234567890123456789012345678901234567890123456789012345-678901234567890123456789012345678901234567890
        last 2 digit    5555555556666666666777777777788888888889999999999000000-000011111111112222222222333333333344444444445
        last 3 digit    5555555555555555555555555555555555555555555555555666666-666666666666666666666666666666666666666666666

2. 出力例2
   • 配列の頭と尻にミスマッチがある場合。
   1. BWAの出力文字列

      • HWUSI-EAS733_0001:1:1:382:685#0 16      chr1    569052  25      36M     *       0       0       ATATAACAAACCCTGAGAACCAAAATGAACGAAAAC    PbYM_`N]Nd]]T]\ObXdaOXLKKRW\STQHHKM[    XT:A:U  NM:i:2  X0:i:1  X1:i:0  XM:i:2  XO:i:0  XG:i:0  MD:Z:0C34T0

   2. 出力例中のミスマッチ等を表す項目1
      • 6列目のCIGAR列の文字列 "36M"
   3. 出力例中のミスマッチ等を表す項目2
      • 19列目のMDタグの文字列 "MD:Z:0C34T0"
   4. ゲノムの配列の取得
      • ゲノムにどのように当たっているかを確認するために、ゲノムの特定部位を切り出します。ゲノムのFASTA(1行50塩基の場合)を使用して簡易的に染色体位置を表示させるコマンド

      • cat -n /panasas/vol01/CIPdata/genomes/hg19/fasta/chr1.fa | awk '{print ($1-2)*50+1 "\t" $2}' | less
        
        568951  ATTACCCCCATACTCCTTACACTATTCCTCATCACCCAACTAAAAATATT
        569001  AAATACAAATTACCACCTACCTCCCTCACCAAAGCCCATAAAAATAAAAA
        569051  ACTATAACAAACCCTGAGAACCAAAATGAACGAAAATCTGTTCACTTCAT
        569101  TCATTGCCCCCACAATCCTAGGCCTACCCGCCGCAGTACTGATCATTCTA

   5. テキストアライメント結果
      • ゲノムとタグの情報を手作業でアライメントしたもの

      • An image of meaning of CIGAR                                    <----------------36---------------->
        An image of meaning of MD tag                                  0C<---------------34--------------->T0
        read                                                            ATATAACAAACCCTGAGAACCAAAATGAACGAAAAC
        aligment                                                         ||||||||||||||||||||||||||||||||||
        569001       AAATACAAATTACCACCTACCTCCCTCACCAAAGCCCATAAAAATAAAAAACTATAACAAACCCTGAGAACCAAAATGAACGAAAATCTGTTCACTTCAT
        last 1 digit 1234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890
        last 2 digit 0000000001111111111222222222233333333334444444444555555555566666666667777777777888888888899999999990
        last 3 digit 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001

3. 出力例3
   • ソフトクリッピングを含む例。
   1. BWAの出力文字列

      • HWUSI-EAS733_0001:1:2:135:1029#0        16      chr8    26267991        37      33M2D2M1S       *       0       0       TATTTGAACTTGAGACCATTGTAAGCATGACCCCCA    e`bVbZdebOdadb]dbd^^^_Y[bbZ`XPZ]b[bb    XT:A:U  NM:i:2  X0:i:1  X1:i:0  XM:i:1  XO:i:1  XG:i:1  MD:Z:33^AA2

   2. 出力例中のミスマッチ等を表す項目1
      • 6列目のCIGAR列の文字列 "33M2D2M1S"
   3. 出力例中のミスマッチ等を表す項目2
      • 19列目のMDタグの文字列 "MD:Z:33^AA2"
   4. ゲノムの配列の取得
      • ゲノムにどのように当たっているかを確認するために、ゲノムの特定部位を切り出します。ゲノムのFASTA(1行50塩基の場合)を使用して簡易的に染色体位置を表示させるコマンド

      • $ cat -n ../CIPdata/genomes/hg19/fasta/chr8.fa | awk '{print ($1-2)*50+1 "\t" $2}' | less
        
        26267951        ATGCGTGTGACCTGTGAAGTGGTGTATTGTCACAGTAGCTTATTTGAACT
        26268001        TGAGACCATTGTAAGCATGACCCAACCTACCACCCTGTTTTTACATATCC
        26268051        AATTCCAGTAACTCTCAAATTCAATATTTTATTCAAACTCTGTTGAGGCA

   5. テキストアライメント結果
      • ゲノムとタグの情報を手作業でアライメントしたもの。

      • An image of meaning of CIGAR                            <--------------33--------------->DD<>S
        An image of meaning of MD tag                           <--------------33--------------->AA<> 
        read                                                    TATTTGAACTTGAGACCATTGTAAGCATGACCC--CCA
        alignmet                                                |||||||||||||||||||||||||||||||||  || 
        26267951        ATGCGTGTGACCTGTGAAGTGGTGTATTGTCACAGTAGCTTATTTGAACTTGAGACCATTGTAAGCATGACCCAACCTACCACCCTGTTTTTACATATCC
        last 1 digit    1234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890
        last 2 digit    5555555556666666666777777777788888888889999999999000000000011111111112222222222333333333344444444445
        last 3 digit    9999999999999999999999999999999999999999999999999000000000000000000000000000000000000000000000000000