Graphフォーマットを使えば環状のリファレンスゲノムを正しく表現できるが、プレーンのFASTA形式には環状のリファレンスゲノムを正しく表現する方法が整備されていない。そのため、環状ゲノムのシーケンシングデータを"線状の"リファレンスゲノムFASTAにマッピングすると、FASTAの末端付近に相当する領域を読んだリードのアライメントが難しくなる（*1）。この環状ゲノムを考慮しないマッピングは、特に古いマッパーで起こりやすい。
末端領域のカバレッジが落ち込むと、アライメント後の分析ステップで複数の問題が生じ得る。１つ目に、リードがアライメントされないため末端付近に変異（SNV、indel、CNV）が起きていても検出されなくなる（false negative）。２つ目に、ペアエンドリードが5'末端と3'末端にジャンプしてマッピングされると、実際は連続した領域にも関わらず、大きな構造変化が起きているとSV callerが勘違いする（false positiveの問題）。３つ目に、アライメントされなかったリードがゲノムのセカンドベストな部位にアライメントされ、そこで偽のバリアントコールを引き起こす（false positiveの問題）。１つ目の問題は、一般に5'末端と3'末端に遺伝子のコード領域がまたがった状態でリファレンスが登録されることはないため、たとえ変異の見逃しがあっても、大きな問題にはなりにくい（目的による）。２つ目の問題も、例えば、"1-7191318 Translocation"というようなゲノム全域に渡るSVコールは目視チェックで除けるため問題にならない。３つ目はより大きな問題かもしれない。オルガネラゲノムはコピー数が変動し、クロモソームよりコピー数が非常に多い場合もある。アライメントされなかったリードがリファレンスの別の領域にアライメントされると、そこで高度にミスアライメントを引き起こす（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) ができる。
#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。
引用
参考
*1
結果は同じように見えるが、線状ゲノムの末端がシーケンスができなくてマッピングされない事とは全く話が異なる。
*2
headerが不明ならsamtools faidxでヘッダー情報を抜き出すか、grep -n ">" 等を実行する。
*3
NGSのマッパーはシーケンスエラーも考慮に入れた高感度なマッパーなので、アルゴリズムにもよるが、少々ミスアライメントがあってもアンマップにはなりにくい。そのため、ベストマッピング部位がリファレンスから意図的に取り除かれると、セカンドベストな部位にアライメントされることがままある。Exomeキャプチャのシーケンシングデータでも全ゲノムをリファレンスにしてマッピングを行うのは、キャプチャ時に少量のexome外領域が入り込み、そのような領域由来のリードがexomeに誤ってマッピングされfalse positiveな情報が発生するのを防ぐためである（他の理由も考えられる）。
