macでインフォマティクス

macでインフォマティクス

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

diginormによるシーケンスデータの軽量化

 

 "digital normalization"という名で発表されたこの手法は、k-merを指標にリードを間引いて、データサイズを軽量化する方法論。データサイズが大きすぎてアセンブルできないサンプルの軽量化に使えるとされる。トリミングターゲットは、低/高のk-merカバレッジのリードになる。k-merのカバレッジが中央値から極端に少ないのはシーケンスエラーと考えられ、k-merのカバレッジが極端に多いにはリピートと考えられる。それらを削って軽量化する。

 ただし低k-merカバレッジ配列はエラーでなくレアなバリアントの可能性もあるし、高カバレッジの配列も単純リピートでなく倍数体ゲノムの極端に似た配列の可能性がある。それらをエラーのk-merと一緒に消してしまう恐れがある。それについて論文の著者がHPで考えを述べている(Why you shouldn't use digital normalization)。 興味がある人はその辺りを理解した上で使ってください。

  本手法はde novo transcrptomeのデータにも使うことが可能とされ、de novo transcriptomeの論文にも使われた実績を持つ。

 

diginorm公式サイト

http://ged.msu.edu/papers/2012-diginorm/

解析にはkhmerを使う。

http://khmer.readthedocs.io/en/v2.0/user/scripts.html#scripts-diginorm 

チュートリアル

A tutorial in basic digital normalization — ANGUS 2.0 documentation

ランに必要なツールのインストール

khmerとscreedが必要。

git clone git://github.com/ged-lab/screed.git 
git clone git://github.com/ged-lab/khmer.git

cd screed
python setup.py install
cd ../khmer
make test
export PYTHONPATH=/root/khmer/python

 khmerとscreedはpipでも導入できます。

ベストなk-merについて

 

ラン  

トランスクリプトームデータ

 1-passノーマライズを行う。

はじめにインターレースのpaired.fqにmergeしておく。

seqtk mergepe R1.fq R2.fq > paired.fq

 

 分布でトリミング。

khmer/scripts/normalize-by-median.py -C 20 -k 20 -N 4 -x 5e8 paired.fq

CとKの 値について(上と同じリンク)。

 

 

ゲノムデータ

スリーパスノーマライズ (3-pass)を行う。

はじめにインターレースのpaired.fqにmergeしておく(seqtk(リンク))。

seqtk mergepe R1.fq R2.fq > paired.fq

k-merカウントグラフを書く。(countgraph.py)

load-into-counting.py -k 20 -x 5e8 -T 8 graph paired.fq 
  • -k k-mer size to use (default: 32)
  • -T Number of simultaneous threads to execute
  • -x upper bound on tablesize to use; overrides --max-memory-usage/-M.

 

1-pass 分布でリードをトリミング。

normalize-by-median.py ; Discard sequences based on whether or not their median k-mer abundance lies above a specified cutoff.

normalize-by-median.py -C 20 -k 20 -N 4 -x 5e8 paired.fq
  • -k k-mer size to use (default: 32)
  • -C when the median k-mer coverage level is above this number the read is not kept. (default: 20)
  • -N number of tables to use in k-mer countgraph
  • -x upper bound on tablesize to use; overrides --max-memory-usage/-M.

paired.fq.keepが出力される。試したデータでは、リードは287,956から273,788に減少 (5%カット)。 出力されるリード数は揃っていない。あとでペアに修正する。

 2-pass 今度はlow coverageのリードをトリミングする。

filter-abund.py ; Trim sequences at a minimum k-mer abundance.

filter-abund.py -T 8 -C 2 graph paired.fq.keep
  • -C CUTOFF, --cutoff CUTOFF Trim at k-mers below this abundance. (default: 2)
  • -T Number of simultaneous threads to execute (default: 1)

R1.fastq.keep.abundfiltとR2.fastq.keep.abundfiltが出力される。

paired.fq.keep.abundfiltが出力される。リードは60%ほどトリムされた。 

 3-pass 最後にもう一度-C 5でトリミング。

normalize-by-median.py -C 5 -k 20 -N 4 -x 1e8 paired.fq.keep.abundfilt

paired.fq.keep.abundfilt.keepが出力される。

最後に ペアリードとペアじゃないリードに分ける。

split-paired-reads.py -0 orphan -1 R1_trimming -2 R2_trimming paired.fq.keep.abundfilt.keep 
  • -0 output_orphaned Allow "orphaned" reads and extract them to this file (default: None)
  • -1 Output "left" reads to this file (default: None)
  • -2 Output "right" reads to this file (default: None)

 R1とR2の行数は同じになる。ペアじゃないオーファンなリードが20%近くを占めた。

 

 

KMCを使ってk-mer頻度を計算(高速なk-merカウントツール KMC)。

 ↓

20merでヒストグラムを書く。ややエラーの多いfastqを使っている。

diginorm処理前

f:id:kazumaxneo:20180221164944j:plain

頻度が少ない(エラー、low coverage、レアバリアント、コンタミ)20merが非常に多い。

 ↓

diginorm処理後

f:id:kazumaxneo:20180221163327j:plain

頻度が少ない20-mer配列が消えた。

 

 

diginormの前に、エラー訂正もやっておいた方がいいかもしれません。karatならrecall、precision評価ツールもついています。エラーコレクションツール karect

下は"before"をkaratでエラー訂正したfastqのヒストグラム

f:id:kazumaxneo:20180221163254j:plain 

まだ多いが、beforeと比べて頻度が少ないk-merの数が大きく減っている。捨てる前に、修復できるなら修復して、捨てられないようにしようという発想です。

 

 

引用

A Reference-Free Algorithm for Computational Normalization of Shotgun Sequencing Data

C. Titus Brown, Adina Howe, Qingpeng Zhang, Alexis B. Pyrkosz, Timothy H. Brom

Publication date: 2012