2023/03/12 誤字修正
ロングリードを使ったゲノムプロジェクトが爆発的に増えており、現在ではほとんどのゲノム解読プロジェクトでロングリードのシークエンスが主要に使用されています。ロングリードのゲノムプロジェクト増加に伴って、想定しない汚染生物のリードが含まれるシークエンシングデータに出会う頻度も高まっている可能性があります。汚染が見つかってしまった時は、ひどい時はサンプル調整からやり直すか、もしくは宿主と汚染生物（共生菌含む）のロングリードを分離する方法が求められます。バイオインフォマティクスのジャーナルでは、既にシークエンシングされたデータからの除染法を提案している研究がいくつか報告されています。手法としては、配列の特徴量から配列を生物ごとに分ける仕分ける"Binning"がメジャーだと思われます。このBinningはメタゲノムのcontigを種や株ごとに分けるMetagenomic Binningのアルゴリズムに似ていますが、（ショートリードの）メタゲノムのBinningがアセンブル後のコンティグを仕分けるために使われるのに対し、こちらのBinningはアセンブル前の生のロングリードを直接分ける点で異なります。ロングリードは数kbと長い距離に渡ってDNA配列を読めるため、コンティグ同様、配列の特徴量の類似性からリードを直接Binningすることが可能ということです。ただしそのためには、コンティグよりも高いエラー率、膨大なリードを扱う計算効率の高さが要求されます。
ロングリードを分類する他の方法として、リードにtaxaを直接アサインする方法、GC含量やk-mer組成によってリードをラフにクラスタリングする方法、バルクのリードセットを正規化して種それぞれのカバレッジをならす方法、などの派生が存在します。これらの方法は複雑な微生物叢のデータセットでは有効な可能性があります。しかしゲノムプロジェクトで汚染が見つかる時は、あくまで宿主のDNAが主要に読まれており、同時に汚染生物種が少しだけ読まれていることが大半です。これは環境サンプルが一般的に裾の長いロングテールな菌組成をしている事とは対照的です。差があるのは、ゲノムプロジェクトで読まれるサンプルは人間がコントロールしているからです。即ちゲノムプロジェクトで読まれるDNAは、汚染がなくて体細胞ではない組織が選ばれたり、先行して汚染が除去されて無関係の生物が徹底的に除外されています。汚染生物を減らす人為的な選択圧があるため、汚染があっても多くはないというわけです。しかし少ないからといって汚染生物のリードがある時に無視するわけにはいきません。無視して進めると、公開されたアセンブリ配列の中に汚染生物の配列が混ざる危険性が高くなります。アセンブル後早いうちかアセンブル前に除外する必要があります。
最初に説明した通り、ロングリードを仕分ける高度なアプローチが既に存在します。汚染がある時、そのような実装を試すのも１つです。しかしここではより簡単な方法として、既に存在するホストか汚染生物のゲノムにマッピングし、マッピングの有無によって分離することを考えます。このような方法はあまりに簡単であるため、新しい論文には報告されません。本ブログのようなサイトで扱うに適したテーマです。
間違ったリファレンスへのマッピングテスト
マッピングによってホストと汚染生物の混ざったロングリードを生物ごとに分離できるのかそれとも難しいのかを体感するため、ヒトゲノム（CHM13 *1）のONTとPacBioのシークエンシングデータを敢えてイーストゲノムのリファンスにマッピングしてみます。公開データのSRR12618296とSRR11292120をダウンロードして使います。どちらもCHM13を読んだヒトのシークエンシングデータです。ONTリードは10000リードだけサンプリングして使っています。マッピングにはminimap2を使います。
#ONT
minimap2 -ax map-ont -t 32 Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa SRR12618296.fastq.gz\
|samtools sort -@ 8 - > SRR12618296.bam
#PacBio HiFi
minimap2 -ax map-hifi -t 32 Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa SRR11292120.fastq.gz\
|samtools sort -@ 8 - > SRR11292120.bam
（ONTのリードはmap-ontプリセット、PacBio HiFiリードはmap-hifiプリセットを使用）
bmファイルをIGVで視覚化してMtアセンブリ領域を見たのが下の図です。上側のトラックがONT、下側のトラックがPacBioです。ONTのリードは散発的にアラインメントされてますが、PacBioのリードは全くアラインされていません。
ヒトゲノムを読んだリードをイーストゲノムにマッピングしているため、アラインメントはいずれもFalse Positive（FP）であり、ランダムアラインメントと考えられます。つまり下のトラックの状態が正しいのですが、ONTのリードはなぜ一部のリードがイーストゲノムにアラインメントされているのでしょうか？原因はマッピングのパラメータにあると考えられます。ONTではエラー率の高いリードをアラインするために "-x map-ont"プリセットを使いました。これは ノイズの多いロングリードをアラインするための設定で、エラーが多いため、感度がPacBioのリードより上げられています（精度をできるだけ維持したまま）。また、Mtゲノムには呼吸鎖複合体の遺伝子があり、これらは生物間で保存性が高い遺伝子です。Mtにあるヒトとイースト間で保存性の高い遺伝子領域だけ、"-x map-ont"プリセットではFPのアラインメントが起きていると推測されます。
Mt領域をさらに拡大すると白色のリードが見つかります（下画像）。白くなったリードは、IGVではMAPQ（マッピングクオリティ）が0のリードです。アライナーによっても異なりますが、通常はゲノムの２か所以上の領域に同じ品質でマッピングされたリードにMAPQ==0がアサインされます。マッピングされている事は間違いないのですが、２か所以上に当たることから生物学的な情報を取り出すことは出来ないため、信頼性ゼロの０がアサインされています。IGVで脱色されて表示されるのはそのためです。
補足；MAPQの最大値は255（ただし255は通常マッピング品質が利用できないことを示す）。MAPQはphredスコアと同じようにエラー率を負の対数変換したものなので、MAPQ=30なら、マッピングのエラー率はsamフォーマットの定義から0.1%となる。ただし経験則的マッピングの一意さやリードの品質スコアから得られる予測値であり、0.1%のエラー率は概ね正しいと考えられるが真に正しいのかはわからない。fastq/samのシミュレーターは神の目線で実行できるので、提供したゲノムについての真のMAPQ値をアサインしたsamを作ることが可能になる。ショートリードのMAPQについてはこちらが参考になる。
ただしMAPQ==0の全てがリピート性のアラインメントであるとは限りません。ミスマッチが多くてトリミングもあるようなアラインメントでも、奇妙な話ですが０になる可能性があります（間違っている可能性あり、MAPQ==0なら信頼度ゼロなのでアラインメントは起きないとも考えられる）。また、マルチマッピングでもMAPQが1以上になるアライナーとリードの組み合わせも存在するようです（*2）。ユニークマッピングだけを取り出すには、samtools vierwの"-q"オプションを使うのではなく、samtools viewで"-F 256"を指定してsam/bamをフィルタリングします。
ここでONTのbamのMAPQ値がどのような分布をしているか調べておきましょう。イーストゲノムにアラインしているので低くなるはずですが。
#bam => sam
samtools view -h SRR12618296.bam > SRR12618296.sam
#retrieve MAPQ field
grep -v @ SRR12618296.sam |cut -f 5 - > MAPQfield
#add "MAPQ"
echo MAPQ |cat - MAPQfield > MAPQfield_add_header
#R
data <- read.table("MAPQfield_add_header", header = TRUE)
hist(data$MAPQ, xlab="MAPQ", ylab="log frequency", breaks = 60, xlim=c(0,20), ylim=c(0,10000)) #hist関数を使って品度をプロット
稀にMAPQが10以上のアラインメントもありましたが、ほとんどのリードのアラインメントにはMAPQ==0がアサインされています。ヒトのONTリードを酵母ゲノムに当てているので、アラインメント可能なリードがあってもindelやミスマッチが多数生じ、それによってMAPQが0まで下がるのだと考えられます。最大でもMAPQ==16でした。（繰り返しますがMAPQが０というのは奇妙に聞こえます。ランダムアラインメントなら排除されるべきなので。）
（実際のショートデータでどのように分布するのか興味がある方はこちらのブログを参照 *3）
さて、試しにMAPQが0のリードだけsamtools viewの-qオプションを使って除いてみます。
#MAPQ1以上
samtools view -@ 8 -bSq 1 SRR12618296.bam > above_MAPQ1.bam
samtools index -@ 8 bove_MAPQ1.bam
IGVに読み込みます。
一番上のトラックが元のbam、一番下のトラックがabove_MAPQ1.bamです。白いリードは無くなった一方で、MAPQが1以上のリードはアラインメントされたままになっています。このようにsamtools viewでMAPQによるフィルタリングができます。
今回のデータではMAPQを17以上にすれば酵母ゲノムにアラインされたヒトリードは完全に除かれるはずです。
#MAPQ4以上
samtools view -@ 8 -bSq 17 SRR12618296.bam > above_MAPQ4.bam
=> 出力されたbamにはheaderしか残っていませんでした。
次に”正しいリファレンスであるヒトのGRCh38アセンブリ”にマッピングした時のMAPQ値分布を調べます。最初と同じONTリードを使います。リファレンスには完全長のCHM13アセンブリを用意すべきですが（完全長だがまだエラーは少しあるらしい）、差は大きくないと予想してGRCh38を使っています。ONTリードは10,000リードだけサンプリングして使っています。
#ONT
minimap2 -ax map-ont -t 32 GRCh38_primary_assembly.fasta SRR12618296.fastq.gz > SRR12618296.sam
#retrieve MAPQ field
grep -v @ SRR12618296.sam |cut -f 5 - > MAPQfield
#add "MAPQ"
echo MAPQ |cat - MAPQfield > MAPQfield_add_header
#R
data <- read.table("MAPQfield_add_header", header = TRUE)
hist(log(data$MAPQ), xlab="MAPQ", ylab="log frequency", breaks = 60, xlim=c(0,60), ylim=c(0,5000)) #hist関数を使って品度の対数値をプロット
興味深い結果です。ヒトのCHM13を読んだONTのリードをminimap2の"map-oint"でGRCh38アセンブリにマッピングしても、一部のリードはMAPQは0になっていました。おそらくリピート性のアラインメントがあるためと考えられます。対数にしてないので見にくいのですが、MAPQが1以上の低い低いカウントも少量見つかります。
２つのヒストグラムから、ヒトと酵母がmixしたONTリードをどちらかのリファレンスにアラインし、MAPQ値でカットオフすると、どの数値で線を引いてもFPとFNが発生すると予想されます。
他にFPのアラインメントを除去する方法として、クリッピングがあるかどうかを指標する方法を思いつきます。ロングリードでFPのアラインメントが発生した時、アラインメントはあくまでランダムの一致であり、リードの部分長でしか有意な一致は発生しない可能性が高くなります。そこで、クリッピングされた領域を持つアラインメントのリードを排除するという方法です。
sam formatドキュメントのCIGAR文字列の説明を載せておきます。Sが"soft-clipping"、Hがhard clippingです。
(sam formatより)
たとえば下の画像のFPのアラインメントを見てみます。samのFLAGは "46385S48M6I1M6I6M1I4M1I6M1I23M1I45M4I8M7699S"でした。左端46,385塩基と右端7,699塩基がソフトクリップ（S）されていることを示しています。アラインメントがある領域にも複数の挿入（I）と欠失（D）が起こっています。
（IGVでsoft clip部分を表示するにはView => preference => alignment => show soft clipped basesにチェックをつける。）
クリッピングされたリードを除くために、クリッピングの長さまで考慮できるsamclipを使います。
#samclip 100bp以上
samtools view -@ 4 -h SRR12618296.bam\
|samclip --ref Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa --max 100 \
|samtools sort > non_clipped.bam
#もしくはawkを使って文字列SかHがあるアラインメント行を除く（~は部分一致）。
samtools view -@ 4 -h SRR12618296.bam | awk '$6 ~ /S|H/ {next} {print}' |samtools sort > softclip_filtered.bam
samclipの"--max 100"は100塩基以上のクリッピングまで許容してsamで返すという意味です。よって101-bp以上のclippingがあるリードだけフィルタリング対象になります。--invertをつけると、クリッピングがないリードだけが返されます。このコマンドで127リード中67リードが排除されました。この方法の欠点は、ゲノムが不完全でコンティグ末端からはみ出しているリードにもソフトクリッピングは適用される点にあります。また、ロングリード全長でFPのアラインメントが起こっているリードも除くことができません。オフターゲットのFPのロングリードだけを除ける完全方法ではないことになります。
まとめ
MAPQで線を引くとFPとFNが必ず発生します。クリッピングされたリードを除けば推定されるFPのアラインメントを減らせるが、完全ではありません。
今回紹介した方法は宿主か汚染生物のドラフトゲノムがすでに構築されていることが前提であり、新規生物のゲノムアセンブリで利用できる状況は限定的である点に注意が必要です。De novoアセンブルして得られたcontigに分類やGC、カバレッジ、k-merなどをアサインして総合的に判断するのが確実と考えられます。
引用
The Sequence Alignment/Map format and SAMtools.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R; 1000 Genome Project Data Processing Subgroup.
Bioinformatics. 2009 Aug 15;25(16):2078-9
参考
*1
https://www.jstage.jst.go.jp/article/jsbibr/3/1/3_jsbibr.2022.primer2/_html/-char/ja
*2
http://biofinysics.blogspot.com/2014/05/how-does-bowtie2-assign-mapq-scores.html
*3
Understanding MAPQ scores in SAM files: does 37 = 42? — ACGT
関連
