macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

環状ゲノムのシーケンシングデータのマッピングを改善する CircularMapper

 

 Graphフォーマットを使えば環状のリファレンスゲノムを正しく表現できるが、プレーンのFASTA形式には環状のリファレンスゲノムを正しく表現する方法が整備されていない。そのため、環状ゲノムのシーケンシングデータを"線状の"リファレンスゲノムFASTAマッピングすると、FASTAの末端付近に相当する領域を読んだリードのアライメントが難しくなる(*1)。この環状ゲノムを考慮しないマッピングは、特に古いマッパーで起こりやすい。

 末端領域のカバレッジが落ち込むと、アライメント後の分析ステップで複数の問題が生じ得る。1つ目に、リードがアライメントされないため末端付近に変異(SNV、indel、CNV)が起きていても検出されなくなる(false negative)。2つ目に、ペアエンドリードが5'末端と3'末端にジャンプしてマッピングされると、実際は連続した領域にも関わらず、大きな構造変化が起きているとSV callerが勘違いする(false positiveの問題)。3つ目に、アライメントされなかったリードがゲノムのセカンドベストな部位にアライメントされ、そこで偽のバリアントコールを引き起こす(false positiveの問題)。1つ目の問題は、一般に5'末端と3'末端に遺伝子のコード領域がまたがった状態でリファレンスが登録されることはないため、たとえ変異の見逃しがあっても、大きな問題にはなりにくい(目的による)。2つ目の問題も、例えば、"1-7191318  Translocation"というようなゲノム全域に渡るSVコールは目視チェックで除けるため問題にならない。3つ目はより大きな問題かもしれない。オルガネラゲノムはコピー数が変動し、クロモソームよりコピー数が非常に多い場合もある。アライメントされなかったリードがリファレンスの別の領域にアライメントされると、そこで高度にミスアライメントを引き起こす(FPコールのホットスポットになる)可能性がある(*3)。このようなアライメントの誤りに由来するリードは、オルガネラゲノムのカバレッジの深さから、バリアントのフィルタリングを行っても信頼性の高いコールとして残る可能性がある(strand biasが強かったり、カバレッジが異常ならフィルタリングできるかもしれない)。この問題が根深いと表現したのは、例えば第三者が解析したデータにだけそのようなFPコールがあると、それが特異的なコールとして抽出され、candidatesのリストに残ってしまうことがありえるためである。

  CircularMapperは環状ゲノムのシーケンシングデータのリファレンス末端のアライメントを改善できるツール。延長したリファレンスにマッピングし、その後リアライメントして補正することで環状ゲノムのマッピングを改善する。Documentsを読む限り、本ツールはbwa alnで検証されている。

 

マニュアル

https://circularmapper.readthedocs.io/en/latest/contents/userguide.html

 

インストール

mac os10.13 java8でテストした。

依存

  • Java Runtime Environment 8 or later

本体 Github

#Anaconda環境ならcondaで導入できる。それ以外はコンパイル済み実行jarファイルをダウンロードするか、ビルドする。
conda install -c bioconda circularmapper 

circulargenerator

$ circulargenerator 

usage: CircularGenerator

 -e,--elongation <ELONGATION>   the elongation factor [INT]

 -h,--help                      show this help page

 -i,--input <INPUT>             the input FastA File

 -s,--seq <SEQ>                 the names of the sequences that should to

                                be elongated

Missing required options: i, e, s

> realignsamfile 

$ realignsamfile 

usage: RealignSAMFile

 -e,--elongation <ELONGATION>                          the elongation

                                                       factor [INT]

 -f,--filterCircularReads <FILTER>                     filter the reads

                                                       that don't map to a

                                                       circular identifier

                                                       (true/false),

                                                       default false

 -h,--help                                             show this help page

 -i,--input <INPUT>                                    the input SAM/BAM

                                                       File

 -r,--reference <REFERENCE>                            the unmodified

                                                       reference genome

 -x,--filterNonCircularSequences <FILTERSEQUENCEIDS>   filter the sequence

                                                       identifiers that

                                                       are not circular

                                                       identifiers

                                                       (true/false),

                                                       default false

Missing required options: i, r, e

 

ラン

まず500-bp modifyしたリファレンスゲノムを作る。ターゲットはリファレンスのミトコンドリアゲノムとする(FASTAのヘッダーで指定する *2)。

circularGenerator -e 500 -i ref.fa -s "chrMT"

 modifyされたfasta (ref_500.fa) ができる。

 

 modifyされたfastaマッピングする。

#indexをつける。
bwa index ref_500.fa

#bwa alnでmapping
bwa aln -t 8 ref_500.fa merged.fastq -n 0.04 -l 32 > output.sai
bwa samse -r "specify_read_group_shere" ref_500.fa output.sai merged.fastq > output_circmapper.sam

samファイルができる。 

 

リアライメントを実行してオリジナルsamにマッピングした状態を再現する。オリジナルfastaを指定する。 

realignSAMFile -e 500 -i output_circmapper.sam -r ref.fa -f true -x false

 bamができる。 

 

bamをソートする。

samtools sort -@ 8 output_circmapper_realigned.BAM -o output_circmapper_realigned.sorted.bam
samtools index output_circmapper_realigned.sorted.bam

 

IGVで末端付近を表示した。上がCircularMapperでリアライメントを実行したbam、下がbwa aln=> samtoolsで特別な作業なしに作ったbam。

f:id:kazumaxneo:20180901200623p:plain

 

引用

GitHub - apeltzer/CircularMapper: A method to improve mappings on circular genomes, using the BWA mapper

 

参考

*1

結果は同じように見えるが、線状ゲノムの末端がシーケンスができなくてマッピングされない事とは全く話が異なる。

 

*2

headerが不明ならsamtools faidxでヘッダー情報を抜き出すか、grep -n ">" 等を実行する。

 

*3

NGSのマッパーはシーケンスエラーも考慮に入れた高感度なマッパーなので、アルゴリズムにもよるが、少々ミスアライメントがあってもアンマップにはなりにくい。そのため、ベストマッピング部位がリファレンスから意図的に取り除かれると、セカンドベストな部位にアライメントされることがままある。Exomeキャプチャのシーケンシングデータでも全ゲノムをリファレンスにしてマッピングを行うのは、キャプチャ時に少量のexome外領域が入り込み、そのような領域由来のリードがexomeに誤ってマッピングされfalse positiveな情報が発生するのを防ぐためである(他の理由も考えられる)。