k-mersからゲノムの類似性を高速計算する kWIP

　DNAシークエンシングの主な用途は、試料の遺伝的構成を互いに比較して共通性を同定し、したがって関連性を検出するか、またはその差を利用して機能を解明することである。最初に、仮定された遺伝的系統および複製を確認するか、またはサンプルを家族、集団および種にグループ化することを試みる。幅広いサンプルのコレクション間の遺伝的関連性を推定するには、バイアスを避け、サンプルあたりのコストを最小限に抑える必要がある。

　今日、集団ゲノミクスの研究の大部分は、NGSを用いて行われている[論文より　ref.1]。全ゲノムシーケンスデータを分析するために一般的に用いられる方法は２つの補間する方法に基づいている。すなわちリファレンスとするゲノムをアセンブリし、リシーケンスしたサンプルをこのゲノムにマッピングしてバリアントをコールする。このアプローチは、モデル生物においては機能的であるが、理想的な方法ではない。リファレンスとなる個体の選択は大部分がランダムであるため、リファレンスゲノムのアセンブリには時間とコストがかかるし[ref.2,3]、不適切なリファレンスゲノム配列に基づく分析は、リファレンスと異なっているか欠けている場合などバイアスを受けやすい[ref.4,5]。遺伝的関連性を測定するアラインメントフリーの方法は、このリファレンスゲノムバイアスを克服するのに役立つ。 

　懸念事項のもう一つは、サンプルの識別となる。最近のレビュー[ref.6]は、誤った標本化が驚くべき速度で起こっていることを発見した（pubmed）。集団遺伝子プロジェクトにおけるサンプル数の増加に伴い、正確かつ一貫したメタデータの問題は、技術的（mix-up）と生物学的（misidentification）といういくつかのレベルで発生する。大規模なフィールド、および遺伝子バンクコレクション全体がDNAシーケンス解析されている。ラボから実験室までのサンプル処理では、シーケンスリードファイルや最終的なデータ・リポジトリへのアップロードまで、サンプルとファイルが混在して誤って表示される可能性が十分にある。この問題は、そのような取り組みのしばしば非常に協力的な性質によってしばしば悪化する。しかしながら、いくつかの誤った同定は、（compilo）種複合体における倍数性、潜在性の種、またはサブのゲノムレベルの変化など、分子遺伝学的分析がなければ事実上検出されない可能性がある[ref.7]。残念なことに、この隠された変異の多くは、ショートリードシーケンシングデータからゲノム全体の遺伝的関連性を計算するための前述の現在のベストプラクティスに従うことによって容易に見過ごされる。誤ったサンプル同定および発散レベルの過小評価は、GWASに使用するサンプルおよび集団、欠落している遺伝性は実際にはメタデータ内に存在する可能性がある。

　アライメントフリー配列比較は、配列アライメントのプロセスを回避することにより、これらの困難に対処することを目的とする（一部略）。いくつかのアライメントフリー配列比較ツールは、依然として基礎となるゲノム配列の事前知識を必要とするため、デノボツールとしてのそれらの使用が妨げられる。最近、デノボ比較を可能にするいくつかのアルゴリズムが公開されている。これらのエクステンションは全てシーケンシング・リードからの系統発生的関係を再構築しようと試みる。 Spaced [ref.13,16]は、系統学的再構成の性能を向上させるために、間隔の取られたシード（小さなk-mersを互いから、または点在する無視された塩基と短い距離）でJensen-Shannon distanceを用いる。 Cnidaria [ref.17]とAAF [ref.18]はJaccard距離を用いて系統を再構築するが、mash [ref.19]はJaccard距離のMinHash近似を用いて同じ効果を得る。

　最も確立され、研究されたアライメントフリーの配列比較測定基準の1つは、D2統計である[ref.8,10]。これは、k-merマッチの数によって2つの配列間の差異を測定する。まず、全てのk量体を各配列中で計数し、計数ベクターに記録する。次に、これらのベクトルの差を測定する。元のD2統計量の場合、これは単にベクトル積を構築することによって達成される。観測されたk-merスペクトラウムをモデル化し相関させることによって精度を改善することを目的としたD2統計量のいくつかの派生、例えばD*2、Ds2が長年にわたって開発されている[ref.8,20-23]。これらの統計は次世代シークエンスデータ[ref.24]に拡張され、メタゲノム比較[ref.25]にうまく適用されたが、D2やDS2のようなD2統計的導関数は計算速度が遅く、定義が難しい。

　ここでは、k-merベースの配列比較に2つの概念を導入して組み合わせる遺伝的関連性を推定するための新しい指標であるk-mer Weighted Inner Productを提示する。 D2統計量と同様に、類似性測度はk-mer数の内積であるが、まず、各k-merを比較するのではなく、むしろすべてのk-mersを確率データ構造にハッシュする：Sketch[ref.26]。得られたSketchは、事実上、k-merカウントのベクトルである。重要なことに、すべてのサンプルのSketchは一定のサイズを持つ。第2に、情報理論上の重みづけを導入して、関連する遺伝子信号をノイズよりも高くする。次に、集団全体の頻度から導出された情報量によって重み付けされた、k-merカウント間の内積によって対の類似性が計算される。この手順は、シーケンシングリードから直接、メトリックであるk-mer加重内積を計算するソフトウェアツール（kWIP）で実装されている。著者らはシミュレーションと、公開されたデータセットの再分析によって、kWIPがサンプル間の遺伝的関連性を素早く正確に検出できることを示す。

マニュアル

https://kwip.readthedocs.io/en/latest/

インストール

ubuntu14.04に導入した。

Github

##依存
sudo apt-get install zlib1g-dev cmake build-essential git 

git clone https://github.com/kdmurray91/kWIP.git
cd kWIP

mkdir build && cd build
cmake .. -DCMAKE_INSTALL_PREFIX=${HOME}
make 
make test 
make install
cd bin/
./kwip

> ./kwip 

$ ./kwip

kwip version 0.2.0-13-g3cf8a9e

USAGE: ./kwip [options] hashes

OPTIONS:

-t, --threads Number of threads to utilise. [default N_CPUS]

-k, --kernel Output file for the kernel matrix. [default None]

-d, --distance Output file for the distance matrix. [default stdout]

-U, --unweighted Use the unweighted inner proudct kernel. [default off]

-w, --weights Bin weight vector file (input, or output w/ -C).

-C, --calc-weights Calculate only the bin weight vector, not kernel matrix.

-h, --help Print this help message.

-V, --version Print the version string.

-v, --verbose Increase verbosity. May or may not acutally do anything.

-q, --quiet Execute silently but for errors.

Each sample's oxli Countgraph should be specified after arguments:

./kwip [options] sample1.ct sample2.ct ... sampleN.ct

ビルド済みのバイナリもダウンロードできる。

https://kwip.readthedocs.io/en/latest/installation.html 

ラン

rice 3000 genome project (pubmed) のシーケンスを使用して、流れを確認する（リンク）。 

ワーキングディレクトリの作成

mkdir ~/rice_kwip 
cd ~/rice_kwip

公式ではvirtualenvwrapperでpythonの固有環境を作って作業しているが、ここでは使わない。

pip install srapy 
pip install khmer

ヒアドキュメントでダウンロードリストファイルを作成。いずれもシングルエンドのシーケンスデータ。

cat >sra_run_ids.txt <<EOF 
ERR626208 
ERR626209 
ERR626210 
ERR626211 
ERR619511 
ERR619512 
ERR619513 
ERR619514 
EOF

 シーケンスデータダウンロードディレクトリを作成し、SRApy（リンク）でダウンロードを実行。合計4.3GBある。

mkdir sra_files
get-run.py -s -f sra_run_ids.txt -d sra_files/

（追記）sra-toolikitを使ってもダウンロードできる。asperaを使うなら"-a"オプションでascpの実行ファイルパスと作成した秘密鍵の場所を指定する。

prefetch --option-file sra_run_ids.txt --max-size 1000000000

#終わったらsra_files/に移動
mv ~/ncbi/public/sra/ERR626* sra_files/

 fastqへの変換。fastq-dumpコマンドがないならbrew  install sratoolkitでSRA toolkitを導入する。出力はgz圧縮する。forで回して一括変換する。

mkdir fastqs
for srr in $(cat sra_run_ids.txt) 
do
  fastq-dump \
  --split-spot \
  --skip-technical \
  --stdout \
  --readids \
  --defline-seq '@$sn/$ri' \
  --defline-qual '+' \
  sra_files/${srr}.sra | \
  gzip -1 > fastqs/${srr}.fastq.gz 
done

fastq/に抽出されたfq.gzが保存される。

hash計算ディレクトリを作成し、khmerのload-into-counting.pyで

mkdir hashes
for srr in $(cat sra_run_ids.txt) 
do 
 load-into-counting.py \
  -N 1 \
  -x 1e9 \
  -k 20 \
  -b \
  -f \
  -s tsv \
  hashes/${srr}.ct.gz \
  fastqs/${srr}.fastq.gz 
done

全hashファイルを指定してkwip実行。

kwip -v -t 2 -k rice.kern -d rice.dist hashes/*.ct.gz 

  kernel matrixファイルとdistance matrixファイルが出力される。

> cat kernel matrix

$ cat rice.kern

ERR619511 ERR619512 ERR619513 ERR619514 ERR626208 ERR626209 ERR626210 ERR626211

ERR619511 6.11354e+08 3.61427e+08 3.69794e+08 3.54549e+08 1.79463e+08 1.78418e+08 1.86159e+08 1.7704e+08

ERR619512 3.61427e+08 5.58085e+08 3.51861e+08 3.3726e+08 1.67822e+08 1.67005e+08 1.74017e+08 1.65437e+08

ERR619513 3.69794e+08 3.51861e+08 5.79062e+08 3.45128e+08 1.72347e+08 1.7169e+08 1.78807e+08 1.70043e+08

ERR619514 3.54549e+08 3.3726e+08 3.45128e+08 5.41048e+08 1.63905e+08 1.63253e+08 1.70158e+08 1.61586e+08

ERR626208 1.79463e+08 1.67822e+08 1.72347e+08 1.63905e+08 4.81631e+08 2.83639e+08 2.94275e+08 2.81749e+08

ERR626209 1.78418e+08 1.67005e+08 1.7169e+08 1.63253e+08 2.83639e+08 4.78399e+08 2.94243e+08 2.80658e+08

ERR626210 1.86159e+08 1.74017e+08 1.78807e+08 1.70158e+08 2.94275e+08 2.94243e+08 5.08735e+08 2.91097e+08

ERR626211 1.7704e+08 1.65437e+08 1.70043e+08 1.61586e+08 2.81749e+08 2.80658e+08 2.91097e+08 4.73382e+08

>  cat rice.dist 

$ cat rice.dist 

ERR619511 ERR619512 ERR619513 ERR619514 ERR626208 ERR626209 ERR626210 ERR626211

ERR619511 0 0.873197 0.870041 0.87582 1.15695 1.15766 1.15429 1.15837

ERR619512 0.873197 0 0.872979 0.87891 1.16301 1.16344 1.16053 1.16459

ERR619513 0.870041 0.872979 0 0.875678 1.16073 1.16086 1.15807 1.16208

ERR619514 0.87582 0.87891 0.875678 0 1.16526 1.16543 1.16247 1.1668

ERR626208 1.15695 1.16301 1.16073 1.16526 0 0.904545 0.900558 0.905468

ERR626209 1.15766 1.16344 1.16086 1.16543 0.904545 0 0.8984 0.905801

ERR626210 1.15429 1.16053 1.15807 1.16247 0.900558 0.8984 0 0.902021

ERR626211 1.15837 1.16459 1.16208 1.1668 0.905468 0.905801 0.902021 0

distance matrixが得られたので、最後に図で出力する。

Rのfieldsパッケージも使うので、なければCRANからインストールする。

（R内で "install.packages("fields", dependencies = TRUE") 

img.Rをカレントにコピーする
cp kWIP/utils/img.R .

#描画してPDF保存
Rscript img.R rice

MDS

２つのグループに分けられた。 

引用

kWIP: The k-mer weighted inner product, a de novo estimator of genetic similarity

Kevin D. Murray, Christfried Webers, Cheng Soon Ong, Justin Borevitz, Norman Warthmann.

PLoS Comput Biol. 2017 Sep; 13(9)