deep sequenceされたウィルスのアセンブルツール sparNA

sparNAはウィルスゲノムのアセンブリツール。ウィルスゲノムはRNA ploplymeraseのエラー率の高さなどの要因でhetero genesityが非常に高いため、特別な仕分け方をしない限りpopulation genomeやmeta genomeのデータセットに近い状態でシーケンス解析が行われる。また、ゲノムサイズが小さいために１万カバレッジ~10万カバレッジオーバーでシーケンスされることが多く、一般のk-merアセンブリではエラーと真の多様性を区別することが難しい（エラーのk-mer配列ですら何百回も出てくる可能性がある）。

sparNAはこのような混じり物で、割合に大きな違いがあるウィルスゲノムを個別にアセンブルして、結果をビジュアルで分析できるツールである。内部ではspadesを動かしているので、やっていることはspadesの総当たりに近いが、複数coverageと複数k-merを全自動で総あたりでアセンブルし、アノテーションもかけてくれるのは大変便利である。

インストール

依存

• Java
• Python 3.5
• SPAdes
• Bowtie2
• Samtools
• VCFtools
• BCFtools
• SeqTK
• Argh
• Biopython
• Khmer
• Pandas
• Plotly

brewとpipで一括インストールできる。

brew install python3 spades bowtie2 samtools vcftools bcftools seqtk pip3 install argh numpy biopython khmer pandas plotly

本体 (Github)

ヘルプ

./sparna.py help

ラン

アダプタートリミングなしでアセンブル（あらかじめトリミングしておく）

./sparna.py --lca --fwd-fq paired1.fq --rev-fq paired2.fq --norm-c-list 1,5,20,100 --norm-k-list 21,25,31 --out-prefix output --threads 12

• --fwd-fq　FWD_FQ
• --rev-fq　 REV_FQ
• --norm-c-list　 NORM_C_LIST
• --norm-k-list　NORM_K_LIST
• --out-prefix　OUT_PREFIX
• -t　THREADS(4)

指定したk-merでアセンブルを実行する。またカバレッジを複数指定することで、それをターゲットとして独立にアセンブルを行うらしい。

アダプタートリミングしてからアセンブル。カバレッジのノーマライズは実行しない。

./sparna.py --lca --qual-trim --no-norm --fwd-fq paired1.fq --rev-fq paired2.fq --norm-c-list 1,5,20,100 --norm-k-list 21,25,31 --out-prefix tet --threads 12

•  --qual-trim　False
• --no-norm　False

デフォルトでは Khmerのnormalize-by-median.pyでカバレッジがノーマライズされてアセンブルされるが、--no-normをつけるとノーマライズはスキップされそのままアセンブルが実行される。

結果はhtmlで出力される。

ツールの挙動を確かめるため、6kbのウイルスのシミュレーションデータ（シーケンスエラーがあるだけでほぼ純粋なゲノム）をアセンブルしてみた。

指定したk-merとカバレッジの組み合わせ全てでアセンブリが行われる。上の図で色が違うのはカバレッジまたはk-merサイズの異なるcontigである。例えば濃い緑のcontigはカバレッジ100指定、k-mer=21でアセンブルされたcontigである。カバレッジ指定が異なるcontigもだいたい同じ位置にplotがあるが（緑、黄土色、灰色）、k-merサイズが異なるとGCとcontigサイズが違う位置にcontigがプロットされている（水色と濃い青色）。

領域を指定するとそこが拡大される。

右上のauto scaleボタンかダブルクリックで元に戻る。

blastでhitすればアノテーションも付く。

右下のexport をクリック。

plotlyを使ったcontigの分析リンクが開く(plotly)。

Xボタンを押して消すことで特定のcontigのみ表示できる。

deepシーケンスしたウィルスゲノムなら、k-merのヒストグラムから推定ゲノムサイズを予測し、カバレッジとk-merを振ってランすることになる。

推定カバレッジが10万くらいだとどのような結果が得られるのかテストしてみる。

引用

GitHub - sdwfrost/sparNA: Consensus assembly of deep sequenced RNA virus populations

fastqのクオリティスコアをASCII +64からASCII +33に変換する。

BBtoolsのreformat.shを使えば、ASCII+64でクオリティスコアを計算しているfastqをASCII+33に変換することができる。

シングルリード

reformat.sh in=input.fq out=output_phred33.fq qin=64 qout=33 

ペアリード

reformat.sh in1=input1.fq in2=input2.fq out1=output1.fq out2=output2.fq qin=64 qout=33

BBtoolsで確認する。 

testformat.sh in=output1.fq

>illumina (Q+64)

user$ testformat.sh in=input1.fastq 

illumina fastq raw single-ended 100bp

>sanger (Q+33)

user$ testformat.sh in=phred33_R1.fq 

sanger fastq raw single-ended 100bp

BBtoolsは以前紹介しています。

引用

Biostars

https://www.biostars.org/p/63225/

指定した遺伝子のターゲットエンリッチメントを行う HybPiper

HybPiperは系統解析などを行うために遺伝子領域のエンリッチメントを行うことができるツール。NGSのリードを出発点として、準備した遺伝子配列セット（bait）にリードをアライメントし（BWA, BLAST）、spadesで個別にアセンブルを実行する。出力はcDNA配列とその翻訳産物となる。これを系統解析に使用する。またパラログの出力もサポートしている。

マニュアル

https://github.com/mossmatters/HybPiper/wiki/Installation

インストール

セットアップwiki

https://github.com/mossmatters/HybPiper/wiki/Installation

依存

• Python 2.7 or later
• BIOPYTHON 1.59 or later
• EXONERATE
• BLAST command line tools
• SPAdes
• GNU Parallel
• BWA

• samtools

全てbrewとpipでインストールできる。

sudo pip install biopython 
brew install exonerate 
brew install bwa 
brew install samtools 
brew install spades 
brew install blast 
brew install parallel

GNU parallelは、使用前に以下のように打つ。

parallel --bibtex

 will citeと打ちEnter。

Type: 'will cite' and press enter.

> will cite

これで使用できる。

本体 (Github)

git clone https://github.com/mossmatters/HybPiper.git 
cd HybPiper/ 
./reads_first.py --check-depend #依存ツールのチェック 

usr $ ./reads_first.py --check-depend

HybPiper was called with these arguments:

./reads_first.py --check-depend

blastx found at /usr/local/bin/blastx

exonerate found at /usr/local/bin/exonerate

parallel found at /usr/local/bin/parallel

makeblastdb found at /usr/local/bin/makeblastdb

spades.py found at /usr/local/bin/spades.py

bwa found at /Users/kazumaxneo/bwa-0.7.5a/bwa

samtools found at /Users/kazumaxneo/Documents/samtools-1.3.1/samtools

Package Bio successfully loaded!

Everything looks good!

goodと出る。

実行方法

チュートリアル

https://github.com/mossmatters/HybPiper/wiki

 はじめにテストデータをランする。

cd test_dataset/ 
./run_tests.sh

時間があればまとめます。

引用

HybPiper: Extracting coding sequence and introns for phylogenetics from high-throughput sequencing reads using target enrichment

Matthew G. Johnson, Elliot M. Gardner, Yang Liu, Rafael Medina, Bernard Goffinet, A. Jonathan Shaw, Nyree J. C. Zerega and Norman J. Wickett

Appl Plant Sci. 2016 Jul; 4(7): apps.1600016.

NCBIからvirusゲノムをダウンロードする

Accession IDを使い、virusのゲノム配列(FASTA)をダウンロードする。

NCBIのvirus Genomesに移動する。

左下の方の"Accession list of all viral genomes"をクリックしてvirusのリストをダウンロードする。

このようなリストが入手できる。

user$ head taxid10239.nbr.txt 

## Neighbors data for complete genomes: Viruses (taxid 10239)

## Columns: "Representative" "Neighbor" "Host" "Selected lineage" "Taxonomy name" "Segment name"

NC_003663 KC813499 vertebrates,human Poxviridae,Orthopoxvirus,Cowpox virus Cowpox virus segment  

NC_003663 HQ420896 vertebrates,human Poxviridae,Orthopoxvirus,Cowpox virus Cowpox virus segment  

NC_003663 KC813500 vertebrates,human Poxviridae,Orthopoxvirus,Cowpox virus Cowpox virus segment  

NC_003663 KY463519 vertebrates,human Poxviridae,Orthopoxvirus,Cowpox virus Cowpox virus segment  

NC_003663 KY549148 vertebrates,human Poxviridae,Orthopoxvirus,Cowpox virus Cowpox virus segment  

NC_003663 HQ420897 vertebrates,human Poxviridae,Orthopoxvirus,Cowpox virus Cowpox virus segment  

NC_003663 LT896722 vertebrates,human Poxviridae,Orthopoxvirus,Cowpox virus Cowpox virus segment  

NC_003663 KY549149 vertebrates,human Poxviridae,Orthopoxvirus,Cowpox virus Cowpox virus segment  

Accesion IDだけ抽出する。また同じIDを１つにまとめる。

cut -f 1 taxid10239.nbr.txt |unique -f 1 accesssion_ID.txt

user$ head accesssion_ID.txt

NC_003663

NC_003310

NC_006998

NC_001611

NC_027213

NC_005336

NC_002188

NC_024447

NC_004002

NC_001266

editectを使い、ゲノムをダウンロードする。なければbrewでインストールしておく。

brew install edirect

efetchコマンドをループで回しダウンロードする。

cat accesssion_ID.txt | while read p; do echo $p; efetch -db nucleotide -id $p -format fasta > $p.fasta; done;

引用

Biostar

https://www.biostars.org/p/3297/

巨大なプロテインファミリーのマルチプルアライメントを行うFAMSA

FAMSAは大規模タンパク質ファミリーのマルチプルアライメントを可能にするアルゴリズムを持つ方法論。CPUの並列化に対応しており、数千-数十万のタンパク質ファミリーの高速なマルチプルアライメントが可能になっている。

論文中では、オーサーが定義したおよそ41万５千のアミノ酸配列のベンチマークで、clustal-omegaやMAFFETでは、128GBのメモリを積んだマシンでも24時間で終了しなかったのが、8GBのメモリ使用量で２時間で処理できたと記載されている。

公式サイト

http://sun.aei.polsl.pl/REFRESH/index.php?page=projects&project=famsa&subpage=about

インストール

version1.2でmacに対応している。

Github

git clone https://github.com/refresh-bio/FAMSA.git
cd FAMSA/
make

FAMSAをパスの通ったディレクトリに移動させる。

実行方法

farmsa input.faa output.faa

引用

FAMSA: Fast and accurate multiple sequence alignment of huge protein families

Sebastian Deorowicz, Agnieszka Debudaj-Grabysz & Adam Gudyś

Sci Rep. 2016 Sep 27;6:33964.

高速で高効率なfastqの圧縮ツール DSRC

DSRCはマルチスレッドに対応したfastq（ABI SOLiD, and 454/Ion Torrent）の圧縮ツール。gzipやbzipなどの汎用的な圧縮ツールと比較して15~60%高効率とされる。圧縮・解凍速度も極めて速く、8スレッドで500MB/s出るとされる。

インストール

binaryのダウンロード

http://sun.aei.polsl.pl/REFRESH/index.php?page=projects&project=dsrc&subpage=download

mac向けをダウンロードして、実行権をつけパスを通す。

ソースコード（Gitで管理されている）

https://github.com/refresh-bio/DSRC

user$ dsrc

DSRC - DNA Sequence Reads Compressor

version: 2.00 RC @ 28.03.2014

usage: dsrc <c|d> [options] <input filename> <output filename>

compression options:

-d<n> : DNA compression mode: 0-3, default: 0

-q<n> : Quality compression mode: 0-2, default: 0

-f<1,..>: keep only those fields no. in tag field string, default: keep all

-b<n> : FASTQ input buffer size in MB, default: 8

-o<n> : Quality offset, default: 0

-l : use Quality lossy mode (Illumina binning scheme), default: 0

-c : calculate and check CRC32 checksum calculation per block, default: 0

automated compression modes:

-m<n> : compression mode, where n:

* 0 - fast version with decent ratio (-d0 -q0 -b16)

* 1 - slower version with better ratio (-d2 -q2 -b64)

* 2 - slow version with best ratio (-d3 -q2 -b256)

both compression and decompression options:

-t<n> : processing threads number, default (available h/w threads): 24, max: 64

-s : use stdin/stdout for reading/writing FASTQ data

実行方法

圧縮

dsrc c input.fastq output.dsrc

• c　compression

解凍

dsrc d input.dsrc output.fastq

• d　decompression

4CPUを使い、謝り符号訂正も行いながら高速に圧縮する。

dsrc c -m0 -c -t4 input.fastq output.dsrc

• -t　processing threads number, default (available h/w threads): 24, max: 64
• -c　calculate and check CRC32 checksum calculation per block, default: 0
• -m　compression mode. where n: -m<n> : compression mode, where n: 

　　　m0 - fast version with decent ratio (-d0 -q0 -b16)

　　　m1 - slower version with better ratio (-d2 -q2 -b64)　

　　　m2 - slow version with best ratio (-d3 -q2 -b256)

-- 簡単なベンチマーク --

timeコマンドを使いfastqの圧縮にかかる時間を計測する。gzipはシングルコアしか使えないため、複数コアを使ってgz圧縮可能なpigzとも比べてみる（参考リンク）。

計算には3.46GHz x 12コアのmac pro(2012モデル)を使う。5.6GBのfastqを圧縮する。

DSRC default

time dsrc c -t 22 input.fastq input.dsrc

 5.6GB ==> 788.4MB。要した時間は0分10.6秒。

gzip

time gzip -c input.fastq > input.fastq.gz

5.6GB ==> 755.4MB。要した時間は5分9.0秒。

tar

time tar cvjf input.tar.bz2 input.fastq 

5.6GB ==> 459.4MB。要した時間は8分35秒。

pigz

time pigz -c -p 22 out.fq > out.fq.gz 

5.6GB ==> 756.2MB。要した時間は0分22.9秒。

DSRC 高効率モード

time dsrc c -t 22 -m2 input.fastq input.dsrc

 5.6GB ==> 314.2MB。要した時間は0分22.3秒。

-m2オプションをつければ、高速にかなりの圧縮率を達成できそうである。

クオリティについては" -l "フラグをつけることでlossyな圧縮も可能である。もう使わない可能性が高いデータについては、非可逆圧縮までしてしまってもいいかもしれない（慎重に検討してください）。

*同じラボからfastqのlossyな圧縮が可能なFaStoreも発表されているが、FasToreについては実行ファイルをダウンロードして使っても、ソースからビルドしても、fastqの圧縮途中でsegmentation errorを起こした（cent OS6環境）。

GitHub - refresh-bio/FaStore: FaStore - high-performance FASTQ files compressor

引用

DSRC 2—Industry-oriented compression of FASTQ files

Łukasz Roguski Sebastian Deorowicz

Bioinformatics, Volume 30, Issue 15, 1 August 2014, Pages 2213–2215

Compression of DNA sequences in FASTQ format

Deorowicz, S., Grabowski, Sz

Bioinformatics, 2011; 27(6):860–862

ウィルスゲノムのde novo assemblyツール IVA

RNAウィルスのシーケンスでは、逆転写やPCR増幅のbiasにより極めて不均一なカバレッジになってしまうことが知られている。１本の鎖の中のカバレッジが大きく変動するため、一般のde brujinグラフのアセンブルツールはもとより、鋳型量が異なるmRNAやメタゲノム向けのアセンブルも、それだけではRNAウィルスのアセンブルに完全には対応できないと考えられる（カバレッジが大きく変化する箇所で切れる）。IVAは極端にカバレッジが変動し、また高度にヘテロガスな状態のRNAウィルス向けに開発されたde novo assemblyの方法論である。ペアリード情報を使うため、入力はペアリードである必要がある。

de novo アセンブリのツールはアセンブリ以外のプロセスは行わないツールも多いが、IVAは、イルミナのアダプターとPCRプライマーのトリミングから、final contigの出力まで１つのフローで行うことができる。

マニュアル

https://github.com/sanger-pathogens/iva/wiki

ダウンロード

macにインストールする場合、virtual machineに導入することが推奨されている。試しにmac OS10.12に直接インストールすると、ラン途中のシェルスクリプトでエラーを起こした。最終的にcent OS6サーバーにインストールした。

Python3環境と以下のツールが必要である。

• kmc 
• smalt
• samtools
• MUMmer

kmcは上記リンクからmac向けのbinaryをダウンロードしてパスを通す。他はbrewで導入できる。

本体

pipで導入する。

pip3 install iva

iva -hでヘルプ。

テストラン

iva --test outdir 

user]$ iva --test outdir

Running iva in test mode...

Copied input test files into here: /home/disk1/uesaka/mizutani/outdir

Current working directory: /home/disk1/uesaka/mizutani/outdir

Running iva on the test data with the command:

/usr/local/python/bin/iva --threads 1 --pcr_primers hiv_pcr_primers.fa -f reads_1.fq.gz -r reads_2.fq.gz iva.out

Finished running iva

Looks OK. Final output contigs file is: /home/disk1/uesaka/mizutani/outdir/iva.out/contigs.fasta

/outdir/iva.out/にcontig.fastaができる。 

実行方法

ペアエンドfastqを指定してラン。

iva -f reads_1.fastq -r reads_2.fastq Ouptut_dir --threads 8 -v

• -v　Be verbose by printing messages to stdout. Use up to three times for increasing verbosity.
• -t　INT Number of threads to use [1] 

インターレースのペアエンドfastqを指定してラン。

iva --fr reads.fastq Ouptut_dir --threads 8 -v 

引用

IVA: accurate de novo assembly of RNA virus genomes

Hunt M, Gall A, Ong SH, Brener J, Ferns B, Goulder P, Nastouli E, Keane JA, Kellam P, Otto TD

Bioinformatics. 2015 Jul 15;31(14):2374-6