2023/02/13 誤字修正
ロングリードを使ったメタゲノムシークエンシングが徐々に増えてきています。一般に、ロングリードシークエンシングでは、メタゲノムアセンブリによってショートリードよりも連続性の高いMAGを得ることができます。連続性の高いMAGが得られれば、下流解析の精度を上げることができ、解析も進めやすくなります。
次の図は、SRAアクセッションSRR18495437で公開されているPacbio HiFiメタゲノムシークエンシングをダウンロードし、metaFlyeを使ってアセンブリした結果です。Bandageで可視化しています。いくつか環状MAGが得られていました。一部は分岐していますが、３つだけはそれぞれ１個の閉じた配列として出力されていました。
１個の閉じた配列は１つの細菌の完全なchromosome配列が表現されているはずです。つまり、エラーが出やすいbinningなどの追加操作を必要せずに下流解析を行えることを意味しています。このような結果が得られるのは、精度の高いロングリードを使ったメタゲノムシークエンシングの大きなアドバンテージの１つです（ただしVibrioなど２つchromosomeを持つ細菌の分類もあるので、分類を調べておく必要がある）。
さて、このような１本もしくは少ないcontigから構成される環状配列が得られた時、エラーチェックや株レベルのバリエーションを詳しく調べるために、環状配列の組み立てに使われたロングリードだけを別ファイルとして取り出したいことがあります。
例えば、上の画像の中央の左から２番目の一本の環状配列（黄緑色、3.9Mb）をfasta形式で書き出し、以下のコマンドを実行して、この配列にマッピングされたロングリードを書き出すことを考えます。
minimap2 -ax map-pb -t 20 circular_cointig.fasta SRR18495437.fastq \
|samtools view -@ 4 -F 4 -h - \
|samtools sort |samtools fastq >
MAG1_mapped.fastq
このコマンドでは、アセンブリに使われたロングリード（非圧縮で12GB）を3.9Mbの環状配列にminimap2のmap-pbプリセットでマッピングし、samtools にパイプしてマッピングされたリードだけ（-F 4）をsam形式（-h）で取り出し、再びsamtools にパイプしてソート・bam形式変換を行い（samtools sort）、最後にsamtools fastqを使ってfastqとして書き出しています。
このコマンドにより、3.9Mbの環状配列にマッピングされたリードだけが回収されます。試したところ、得られたfastqのファイルサイズは、元の12GBに対し、およそ1.5GBでした。
このリードを使ってmetaFlyeでアセンブルすると再び3.9Mbの閉じた環状配列が再構成されるはずです。試してみます。
flye --meta --pacbio-hifi MAG1_mapped.fastq --threads 20 -o MAG1_reassembly
結果をBandageで可視化しました。下の図を見ると、さっきより短い3.1Mbの線状配列（黄緑色）が得られ、さらにたくさんの短いcontigが出力されてしまっています。中には主要なコンティグと繋がりがない650-kbもの長さのcontigまで出来てしまっています。
何が起こっているのでしょうか。ここで使用したシークエンシングデータはpacbioのHiFiシークエンシングのデータです。エラーは非常に少なく、間違ったアラインメントは起こりにくいと考えられます。しかし、メタゲノムのシークエンシングであるため、例えこの細菌ゲノムがクローン性の非常いものであると仮定しても、共存している他の細菌のゲノム配列との間で高度に保存された配列を読んだリードがこの細菌ゲノム配列にマッピングされることは起こり得ます（rRNAやほかの高度に保存された遺伝子配列など）。このようなリードがたくさんあれば、上のような結果になることは説明出来ます。
ここでは異種生物由来の可能性が高いリードを除くことを考えます。上手くできれば、再アセンブルでも一本の閉じた配列が得られるはずです。これを高効率に行うにはどうしたら良いのでしょうか？１つのアイデアは、異種生物のリードは全長がリファレンスにアラインメントされず部分的なアラインメントと仮定し、クリッピングしてアラインされているリードをフィルタリングするというものです。
試してみましょう。
minimap2 -ax map-pb -t 20 circular_cointig.fasta SRR18495437.fastq \
|samtools view -@ 4 -F 4 -h - \
|samclip --ref circular_cointig.fasta \
|samtools sort \
|samtools fastq > mapped.nonclipped.fastq
- --max <NUM> Maximum clip length to allow (default=5)
今度のコマンドもマッピングされたリードだけを返すまでは最初のコマンドと同じです。しかしパイプしてsamclipに渡すところから変わっています。samclipはsamのクリッピングを行うツールで、デフォルトでは5塩基以上のクリッピングがあるリードは破棄してリードを返します（アンマップも返すので先にsamtools viewを実行している）。samclipを通すことで、どこかの領域がクリッピングされてアラインメントされているリードは除かれ、全長がアラインメントされたリードだけを回収できます。
＜補足＞上のコマンドは、マッピング => アンマップの除去 => クリッピングリードの除去 => fastqの書き出しを１コマンドで実行しています。この時のパイプの順番ですが、1)マッピング => 2)アンマップの順番で行って最初にアンマップを捨てることで、効率的に計算量を減らすことができます。トータルの実行時間は20秒ほどでした（3990x、20スレッド、PCIe SSD）。
このコマンドで書き出されたリードのファイルサイズは170MBほどでした。最初のコマンドの1.5GBより10倍近く減少していることになります。
このリードを使って再びmetaFlyeでアセンブルしてみます。
flye --meta --pacbio-hifi mapped.nonclipped.fastq --threads 20 -o MAG1_reassembly
結果をBandageで可視化しました。最初と同じ環状配列が再構成されていました。長さも数bpしか変化していませんでした。
このように、クリップされたリードを捨てることで、無関係のリードを効率的に減らすことが出来ました。ちなみに--metaを外して単離ゲノムの設定でアセンブルしても同じ配列が得られました。良かったら参考にしてください。
気になった事
- HiFiシークエンシングリードではうまく出来たが、よりエラー率の高いリードでも同じ事ができるのかは分からない。
- ターゲットの菌の多様性が高かったり、複雑な菌叢のデータではどうなるか分からない。もしかすると、このような方法をリード単位の分類アサインや配列情報量で仕分ける方法と組み合わせる必要があるのかもしれない。
- アダプターを含んでいたり、末端のリードクオリティが低くれば、samclipで過剰にリードが捨てられてしまう可能性がある。事前にトリミングしたリードを使うのが望ましい。
- ランダムサンプリングによってリード数を減らせばデプスの低い配列は消せると考えた。実際に試したところ、非常に断片化したcontigしか得られなかった。
- このような方法で集めたリードを使って再アセンブルすればより連続性の高いbin配列が得られるかもしれない（MetaWrapのre-assemblyの発想に近い。未検討）。
- よりペナルティに厳しいglobal alignmentできるアライナーでも同じ結果が達成できるはず。
引用
https://www.ncbi.nlm.nih.gov/sra/?term=SRR18495437
関連
