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

 

 

ベストな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 mergepe R1.fq R2.fq > paired.fq

seqtkは以前紹介しました(リンク)。

 

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

khmer/scripts/load-into-counting.py -k 20 -x 5e8 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.

khmer/scripts/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が出力される。リードは287,956から273,788に減少 (5%カット)。 出力されるリード数は揃っていない。あとでペアに修正する。

 

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

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

khmer/scripts/filter-abund.py graph paired.fq.keep

-C CUTOFF, --cutoff CUTOFF Trim at k-mers below this abundance. (default: 2)

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

リードは60%ほどトリムされた。

 

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

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

 

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

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%近くを占めた。

 

 

 アセンブルがどう変わるかテストして、追記します。

 

 

 

引用

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