2020 5/16 ダウンロードリンク更新
これまではk-merの値を増やしても、k-merのピークがノイズの中に埋もれてしまうので意味がないと思いこんでいたが、SKESA(紹介)は、ペアエンドリードをマージし、リード長以上の長いk-merも使ってde brujin gaphを構築し、より長いunitigを構築することで、アセンブリの連続性改善を行なっている。 このやり方がうまく行くのは、リード長が長めで、ハイクオリティなショートリードシーケンシングが得られたときに限定されると考えていたが、今年の秋のUnicyclerのアップデートで、SPAdesのロングk-merのサポートが入った(Unicylcer #108)。こういったリクエストが入る背景には、きっとアセンブリの連続性が改善できそう、というデータが出ているからだと思われる。
上の図は、ART シミュレータを使い。Miseq v3のエラープロファイルで250-bpのペアエンドリードをゲノムのx100 fold発生させ(インサートサイズ500、SD50)、BBtoolsでマージし、k-mer頻度を可視化した結果になる(jellyfishでk-merカウント)。マージしたリードの平均サイズは453-bpあったので、k-mer=251でカウントしたが、驚くべきことに、まだk-merのシングルの領域に由来すると思われるピークが見えている(*2)。
SPAdesは repeat resolutionのためにペアエンドリードを使ったscaffoldingプロセスを持つので、k-merを増やしても結果は大差ないかもしれないが、k-mer値を増やし、より長いunitigができるなら、そのほうがミスアセンブリのリスクは減る(ただし、k-mer値を増やすと計算リソースをかなり使う)。実際にアセンブリにどのような影響が出るのかテストするため、SPAdesとUnicyclerでロングk-merを使えるようにして、ソースからビルドし直す。
参考にしたisuues #44 (spadesレポジトリのisssuesより)
インストール
mac 10.14と10.12、ubuntu16.04でテストした。最初に、brewやapt-get、condaなどで導入したspades、unicyclerがあれば、アンインストールするかリンクを外しておく(debian系OSのdocker環境でテストして、うまくいけば本環境にコピーした方がいいかも)。
ビルド依存
- g++ (version 5.3.1 or higher)
- cmake (version 2.8.12 or higher)
- zlib
- libbz2
SPAdes Github release
2018年12月現在最新のSPAdes-3.13.0.tar.gzをダウンロードする。
2020年5月現在で最新の3.14.1をダウンロードする。
wget https://github.com/ablab/spades/archive/v3.14.1.tar.gz
tar -xvf v3.14.1.tar.gz
cd spades-3.14.1/assembler/
コンパイル前に以下の2つを修正
> emacs spades_compile.sh
"-DSPADES_MAX_K="をつけ、k-merを最大301にした(実用上は251程度で十分かも *1)。
こちらも修正
> emacs src/spades_pipeline/options_storage.py
MAX_Kをspades_compile.shの指定値と同じ値に修正する。メモリ上限"MEMORY = 250"も、毎回フラグを立てるのは面倒なので、増やせる環境にあるなら増やしておいて良いかもしれない。(SPAdesはスモールゲノム向けだが、リード数が多くてユニークK-mer数が多くなるとかなりメモリを使う、メモリについてはテストしてないので自己責任で)。他のパラメータも同様。
#buildする
./spades_compile.sh
/binに実行ファイルができる。パスを通す。
UnicyclerもGithubから引っ張ってきてビルドする(Unicycler v0.4.7以降)。
#2018年12月現在、最新は0.4.7
git clone https://github.com/rrwick/Unicycler.git
cd Unicycler
make
> unicycler-runner.py
unicyclerを実行すると、リード長が長い場合、SPAdesに合わせてラージk-merが自動選択される。
Choosing k-mer range for assembly (2018-12-10 22:51:18)
Unicycler chooses a k-mer range for SPAdes based on the length of the input reads. It uses a wide range of many k-mer sizes to maximise the chance of finding an ideal assembly.
SPAdes maximum k-mer: 301
(unusual value, probably indicates custom SPAdes compilation)
Median read length: 261
K-mer range: 53, 91, 123, 149, 173, 191, 209, 223, 235, 247
UnicyclerでMiseq 301-bpのリードを使うと、上のk-merが自動選択された。
テスト
これからk merを増やした時にどんな効果があるのか簡単に調べてみます。記事は2回に分け、下にリンクを貼ります。
参考
*1
Miseqなら301-bp読めるが、リード後半のクオリティがかなり悪くなる。2018年現在の実用上の上限は、どんだけハイクオリティなリードでも、200前後のはず(ONTのようにシーケンサーを民主化してbase callerの選択肢が増えれば、さらなる改善が可能かもしれませんが)(illumina basecaller)。
*2
ベースコールが99.9%正しいと仮定すると、0.999^251=77.7%だが、ブリッジPCRのロングシーケンスはどんどんクオリティが悪くなるので、実際はもっと低いはず。