GAM形式のファイルの見方

[hanzou666]

GAMファイルをJSONに変換して、自分でパースするときに参考になるかもしれないのでメモを残しておく。

% vg view -a aln.gam | jq
{
  "sequence": "GTCACGCCGTACGGCTCCAGCACCAGATCCCTGCCCAGTACCGTGCCGTCAAGCACCGACTCGCATCCGGCGAACGGCGTATCCTGCCACGCCGTGGCCGTATCCGCATGGTTGATCACGAACGTCAACGACTCCCCGCGTGACGGGTAC",
  "path": {
    "mapping": [
      {
        "position": {
          "node_id": "10761",
          "offset": "500",
          "is_reverse": true
        },
        "edit": [
          {
            "from_length": 113,
            "to_length": 113
          }
        ],
        "rank": "1"
      },
      {
        "position": {
          "node_id": "10760",
          "is_reverse": true
        },
        "edit": [
          {
            "from_length": 37,
            "to_length": 37
          }
        ],
        "rank": "2"
      }
    ]
  },
  "name": "246303_chr1_1_2715113_2715262_2707113_NC_011593.1 Bifidobacterium longum subsp. infantis ATCC 15697, complete genome",
  "quality": "HR0dHR4eHR0fHR4iICIhISAhIh0gISIhIh8eHSAjHiIdIR8kHyIfHyQfISAhIiAhHSAdHSIfIiAjIyIeHR0hICEgISMjICIgJR0gHSEeICAdHR0iHx8dJB8dICMdIh0dIB8gIx0dHyMdIh0jJR0fHSMkICEfIB0eIyEf
  "mapping_quality": 60,
  "score": 160,
  "fragment_next": {
    "name": "246303_chr1_0_2714989_2715138_2706989_NC_011593.1 Bifidobacterium longum subsp. infantis ATCC 15697, complete genome"
  },
  "identity": 1,
  "fragment_length_distribution": "5000:0:0:0:1"
}

1行が1クエリ配列を表す。

親フィールド

• sequence: クエリ配列
• path: vg(のJSONバージョン)とだいたい同じ。どこにマップされたかの詳細な情報が書いてある。
• name: クエリ配列の名前
• quality: fastqをクエリにした場合、配列のクオリティ情報がここに乗る
• mapping_quality: 文字通り。
• score: blastとかで出るやつ。
• fragment_next: paired-endでマップしたときの相方の名前がここに入る。もう片方のリードだと、fragment_prev になる。
• identity: クエリ配列がリファレンスにどれだけ類似していたか。 vg map の場合、グローバルアラインメントの結果になるので注意。
• fragment_length_distribution
• secondary_score: 文字通り。ない場合もある。

大雑把にsam形式との対応をとると、↓のような感じでしょうか？

SAM	QNAME	FLAG	RNAME	POS	MAPQ	CIGAR	RNEXT	PNEXT	TLEN	SEQ	QUAL	Other
GAM	name		pathのposition	pathのedit?	mapping_quality	pathのedit?	fragment_next,fragment_prev			sequence		secondary_score, fragment_length_distribution

マッピングのスコアとかクオリティとか類似度でフィルタリングするだけだったら、vg filter , vg sift , jq らへんで簡単に対応できる。
samでいう RNAME,POS,CIGAR に当たるpathは少し複雑なので以下で詳しく述べる。

pathの中には、mappingフィールドがあり、この中にノードごとのマッピング結果がリストで入っている。ここの部分を抜き出して、簡略化して書くと、

"path": {
    "mapping": [
      {"position": {"node_id": "10761"}, "edit": {"from_length": 113}, "rank": 1},
      {"position": {"node_id": "10760"}, "edit": {"from_length": 37}, "rank": 2},
      ....
      }
    ]
  },

となっていて、この部分を翻訳すると、

クエリ配列の先頭はノードID10761に113塩基だけマップされて、残りはノードID10760に37塩基だけマップされた

と読むことができる。mappingの中の1つ1つの要素 {"position", "edit", "rank"} は各ノード(=配列)の情報になっている。Samと対応をとると、positionがRNAME,FLAG(の一部)を表し、edit(のto_length)がCIGAR、を表している。rankはゲノムグラフに特有のもので、クエリから見て何番目に通るノードなのかを表している。

今回はidentityが100%だったのでeditのfrom_lengthとto_lengthは同じ数字だったが、違う部分があると同じになるとは限らない。from_lengthはリファレンスグラフのノードからみた塩基数、to_lengthはクエリ配列からみた塩基数である。(この辺はもっと具体例が必要？)

[hanzou666]

pathの例1

"path": {
    "mapping": [
      {
        "position": {
          "node_id": "45607",
          "offset": "498",
          "is_reverse": true
        },
        "edit": [
          {
            "from_length": 6,
            "to_length": 6
          },
          {
            "from_length": 1,
            "to_length": 1,
            "sequence": "A"
          },
          {
            "from_length": 127,
            "to_length": 127
          },
          {
            "from_length": 1,
            "to_length": 1,
            "sequence": "C"
          },
          {
            "from_length": 16,
            "to_length": 16
          }
        ],
        "rank": "1"
      }
    ]
  },

このクエリ配列は、グラフのノードID45607に逆向きにマップされた。 offsetはSAMでいうPOSに対応するスタート位置の座標(ただし0始まり)を表すので、今回の例ではノードの座標498からスタート。
先頭から6塩基は完全マッチ、次の1塩基はリファレンスと違ってAに置換される、次の127塩基は完全マッチ、その次の1塩基はCに置換、最後の16塩基は完全マッチ。
rankは1しかないので、このクエリはこのノードだけにマップされたということ。

[hanzou666]

pathの例2

  "path": {
    "mapping": [
      {
        "position": {
          "node_id": "16167",
          "offset": "707"
        },
        "edit": [
          {
            "from_length": 146,
            "to_length": 146
          },
          {
            "to_length": 3,
            "sequence": "CTG"
          }
        ],
        "rank": "1"
      }
    ]
  },

このクエリ配列は、グラフのノードID16167にマップされて、スタートは座標707。
先頭から146塩基は完全マッチ、残りの3塩基はリファレンスにはないCTGという配列。

[hanzou666]

pathの例3

  "path": {
    "mapping": [
      {
        "position": {
          "node_id": "30574",
          "offset": "511"
        },
        "edit": [
          {
            "from_length": 80,
            "to_length": 80
          },
          {
            "from_length": 1,
          },
          {
            "from_length": 70,
            "to_length": 70
          }
        ],
        "rank": "1"
      }
    ]
  },

欠失がある場合

[6br]

from_lengthとto_lengthについては、以下の関係にあるようです。

https://github.com/vgteam/vg/blob/master/src/vg.proto から、

message Edit {
	    int32 from_length = 1; // Length in the target/ref sequence that is removed.
	    int32 to_length = 2; // Length in read/alt of the sequence it is replaced with.
	    string sequence = 3; // The replacement sequence, if different from the original sequence.
	}

となっています。このとき、ノード内の塩基配列の部分が、

• matchのとき
  from_length==to_length
• mismatchのとき(置換)
  from_length==to_length、置換される配列がsequenceに入る
• gapのとき、referenceに対してdeletionのとき
  from_lengthの値が入る
• gapのとき、referenceに対してinsertionのとき
  to_lengthの値と挿入される配列がsequenceに入る

となっているようなので、to_lengthとfrom_lengthが一致するか、どちらかの情報だけを出力するように、CIGARに対応するような内容をmatch/mismatchとdeletion/insertionに分解して出力しているように見受けられます。