32  99  chrA    1758    60  100M    =   1995    337 CGACCTATTACCCAGCCATAAGTTGACAACCAAATCCCTTAAACGCACTGTCTATAAGCGCATATAGTAGACGGACCGTTATCCTCATTAAGGACCATTA    5555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555    NM:i:2  MD:Z:67A21C10
18  147 chrA    1762    60  65M3D35M    =   1519    -343    CTATTACCCAGCCATAAGTTGACAACCAAATCCCTTAAACGCACTGTCTATAAGCGCATATAGAAGGACCGTTATCCTCATTCAGGACCATTAATACTGC    5555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555    NM:i:3  MD:Z:65^GAC35

sam 파일은 mapping 후 읽기를 유지합니다 spec
참고 1

한 줄에 하나의 리드가 쓰여 있다.
각각의 리드가 레퍼런스와 비교해, 어떤 차이를 가지고 있는지 알고 싶다.
예를 들어, 맨 위의 리드에는 두 개의 돌연변이가 있습니다.
1824-1825, A=>T
1846-1847, C=>A

각 열은 순서대로 seq name,flag(2 bit로 표현됨),chr,pos,map quality,CIGAR,seq name of mate read,pos of mate read, insert size, sequence, basequality, tags. 가 되는 spec 참조.

CIGAR

CIGAR은 insertion, deletion, softclip 등의 위치를 ​​나타낸다.
예를 들어, indel이 없는 경우 100M(리드 길이=100bp의 경우), 도중에 del이 있는 경우 30M3D70M또한,
N : skipped bases on reference
S : soft clipping
H : hard clipping
P : padding
도 사용된다. Soft clip은 sequence의 항에 포함되어 있지만 align하지 않은 것, hard clip은 align하지 않고 sequence의 항에도 포함되어 있지 않은 부분을 가리킨다. soft clipped reads는 IGV로 보면 미스매치가 이어져 보인다. (img) https://www.google.co.jp/search?q=soft+clipping+reads&client=safari&rls=en&source=lnms&tbm=isch&sa=X&ei=ml_0UqKrCofYkAWL9oGoDg&ved=0CAkQ_AUoAQ&biw=soft&twd isch

주의점은, pos는 position of left-most aligned base이므로, soft clip 된 부분은 포함하지 않는다. (물론 하드 클립도 포함하지 않음)
그러므로

TCGATCGA
1  2

25M5I70M의 경우는 pos는 위의 1을 가리키고, 9M의 경우는 2의 부분을 가리킨다.

MD 태그

MD tag 는 3S6M 와 61G22C15 와 같은 태그로, 참조 배열을 보지 않고, 미스매치나 Deletion 의 배열을 확인할 수 있다. (참고)[ ㅡㅡㅜㅜㅜㅜㅜㅜㅜㅜㅜㅜ bgs포 t. jp / 2012/07 / 데이 ply - 응 rs 단지 g-m-gs. HTML ]
65^GAC35 그렇다면, 61 베이스 레퍼런스와 일치(혹은 Insertion) 후, 레퍼런스에서는 G이지만 read에서는 "무언가"로 바뀌고 있다. 그 후, 22 베이스 일치(혹은 Insertion) 하고, 레퍼런스가 C 의 미스매치가 있고, 그 후 15 베이스 일치(혹은 Insertion) 하고 있다.
61G22C15 의 ^GAC는 참조 서열의 GAC 부분이 deletion이 되어 있음을 나타낸다.

Insertion의 정보는 MD tag에는 포함되지 않는다(빠져 있다). 따라서, 리드의 각 베이스의 레퍼런스 배열에 있어서의 위치를 ​​확인하는 경우는, CIGAR와 조합할 필요가 있다.

덧붙여서, 만약 MD tag가 빠져 있는 경우는 65^GAC35 로 보완해 준다.

1-based 및 0-based

기준 서열에 대한 게놈상의 위치를 ​​나타내는 경우, 1-based 및 0-based position이 있다.

따옴표:

대단히 쉬운 자료

sam file은 1-based position을 사용하고 있다.

이번에는 돌연변이를 나타내는데 0-based position을 사용한다.

MD와 CIGAR를 결합하여 돌연변이 추출

조금 까다롭고, 이런 느낌의 코드가 되었다.
htps : // 기주 b. 코 m/우스야마/ゔぁ리안 t_파르세 r/bぉb/마s테 r/사 m_파르세 r. rb

별로 자신이 없다.

메인 부분은 아래와 같은 느낌.

    @variants = []
    while cigar_idx < cigar.length and md_idx < md.length
      match([current_cigar.type, current_md.type]) do
        with _[:soft, _] {
          @variants << current_cigar.clone
          get_next_cigar.call
        }
        with _[_, :del] {
          current_cigar.ref = current_md.ref
          @variants << current_cigar.clone
          get_next_cigar.call
          get_next_md.call
        }
        with _[_, :mismatch] {
          @variants << Variant.new({:type => :mismatch,
                                    :ref => current_md.ref,
                                    :len => 1})
          current_cigar.len -= 1
          get_next_md.call
        }
        with _[:ins, :match] {
          @variants << current_cigar.clone
          get_next_cigar.call
        }
        with _[:match, :match] {
          if current_cigar.len > current_md.len
            @variants << Variant.new({:type => :match, :len => current_md.len })
            current_cigar.len -= current_md.len
            get_next_md.call
            get_next_cigar.call if current_cigar.len == 0
          else
            @variants << current_cigar.clone
            current_md.len -= current_cigar.len
            get_next_cigar.call
            get_next_md.call if current_md.len == 0
          end
        }
      end
    end

C++에서도 구현했습니다. 여기