macでインフォマティクス

macでインフォマティクス

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

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

 

インストール

依存

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

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

 

本体 (Github)

https://github.com/mossmatters/HybPiper

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と出る。

 

ラン

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

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に移動する。

f:id:kazumaxneo:20171115201348j:plain

左下の方の"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を1つにまとめる。

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万5千のアミノ酸配列のベンチマークで、8GBのメモリ使用量で2時間で処理できたと記載されている。clustal-omegaやMAFFETでは、128GBのメモリを積んだマシンでも24時間で終了しなかったと書かれており、いかに優秀かわかる。

 

公式サイト

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

 

インストール

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

Github

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

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ś

Scientific Reports 6, Article number: 33964 (2016)

 

 

高速で高効率な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

 

 

ラン

圧縮

dsrc c input.fastq output.dsrc
  • c compression

解凍

dsrc t input.dsrc output.fastq
  • t decompression

 

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

dsrc c -m0 -c -t4 input.fastq input.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により極めて不均一なカバレッジになってしまうことが知られている。1本の鎖の中のカバレッジが大きく変動するため、一般のde brujinグラフのアセンブルツールはもとより、鋳型量が異なるmRNAやメタゲノム向けのアセンブルも、それだけではRNAウィルスのアセンブルに完全には対応できないと考えられる(カバレッジが大きく変化する箇所で切れる)。IVAは極端にカバレッジが変動し、また高度にヘテロガスな状態のRNAウィルス向けに開発されたde novo assemblyの方法論である。ペアリード情報を使うため、入力はペアリードである必要がある。

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

 

 

マニュアル

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 M1, 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. doi: 10.1093/bioinformatics/btv120. Epub 2015 Feb 28.

 

GCbiasを考慮したイルミナのシミュレーター ArtificialFastqGenerator

 

 ArtificialFastqGeneratorはカバレッジGCバイアスを考慮可能なNGSリードのシミュレーター。イルミナのペアードエンドfastqに対応している。

  

 

比較表 Biostars

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

ダウンロード

javaの実行ファイルがダウンロードできる。

https://sourceforge.net/projects/artfastqgen/

 

ラン

ピークカバレッジ50の101-bpペアエンドfastqを発生させる。

java -jar ArtificialFastqGenerator.jar -R reference.fa -S -CMP 50 -RL  101 ">chr" -O seq
  • -R Reference genome sequence file, (must be specified). 
  • -S Prefix of the sequence identifier in the reference after which read generation should begin (must be specified).
  • -CMP The peak coverage mean for a region (default = 37.7).
  • -RL The length of each read, (default = 76).
  • -O Path for the artificial fastq and log files, including their base name (must be specified).

 パラメータの詳細は.logに出力される。

 

 

 

リアルデータのクオリティデータを使い、かつphread quality scoreの確率に基づいてエラーを発生させるなら、以下のフラグをつける。

 

インサートサイズ平均600-bp、インサートサイズのSD60、ピークカバレッジ100の301-bpペアエンドfastqを発生。

java -jar ArtificialFastqGenerator.jar -R reference.fa -S -CMP 100 -RL  301 -TLM 600 -TLSD 60 ">chr" -O seq -F1 r1.fq -F2 r2.fq
  • -F1 First fastq file to use for real quality scores, (must be specified if useRealQualityScores = true).
  • -F2 Second fastq file to use for real quality scores, (must be specified if useRealQualityScores = true).
  • -SE Whether to simulate error in the read based on the quality scores, (default = false).
  • -URQS Whether to use real quality scores from existing fastq files or set all to the maximum, (default = false). 
  • -TLM The mean DNA template length, (default = 210).
  • -TLSD The standard deviation of the DNA template length, (default = 60).

 

Trueseqのリアルシーケンスデータを使っている。アダプターに近い先頭数bpと3'末端1bpのクオリティが低くなっている。

>R1

f:id:kazumaxneo:20171106173641j:plain

 >R2

f:id:kazumaxneo:20171106174414j:plain

 

 

引用

Generation of Artificial FASTQ Files to Evaluate the Performance of Next-Generation Sequencing Pipelines

Matthew Frampton , Richard Houlston

PLoS One. 2012; 7(11): e49110.