macでインフォマティクス

macでインフォマティクス

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

samblasterでduplicationリードにタグをつける

 

samblasterは、samファイルのduplicationのリードにタグをつけたり、構造変化の指標となるsplit-alingment readやdiscordant read pairを別ファイルに出力できるツール。samの時点でデータをより分けることで、discordant read pairやsplit-alingment readを使ったlarge indel検出などを劇的に軽量化することが可能になる。2014年に論文として発表された。

 

 

インストール

Github

GitHub - GregoryFaust/samblaster: samblaster: a tool to mark duplicates and extract discordant and split reads from sam files.

 

git clone git://github.com/GregoryFaust/samblaster.git 
cd samblaster
make
cp samblaster /usr/local/bin/.

 

 

duplication、discordant-read、split-readの判定基準(下の方にあります)。

GitHub - GregoryFaust/samblaster: samblaster: a tool to mark duplicates and extract discordant and split reads from sam files.

 

 

 

ラン

 bwa memのアライメントの過程でduplicationにタグをつける。

bwa index -a is input.fa
bwa mem -t 12 -R "@RG\tID:X\tLB:Y\tSM:Z\tPL:ILLUMINA" input.fa *.fastq | samblaster |samtools view -@ 12 -Sb - |samtools sort -@ 12 - > samp.sorted.bam

 上のコマンドはfastqからbwa mem => samblaster => samtools view => samtools sortの流れでbamを作っている。samblasterはbwa memからsamファイルを受け取り、 duplication readsにタグをつけてsamtools viewに渡していることになる。

 今回は以下のようなメッセージがプリントされた。

samblaster: Marked 1874 of 339039 (0.55%) read ids as duplicates using 13344768k memory in 1.056S CPU seconds and 45S wall time.

1874リードがduplicationと判定されている。

 

 

見にくいので、ここから下はsamファイル出力として記載。

 

 

discordant-readとsplit-readは別ファイルに出力する。

bwa mem -t 12 -R "@RG\tID:X\tLB:Y\tSM:Z\tPL:ILLUMINA" input.fa *.fastq | samblaster -e -d samp.disc.sam -s samp.split.sam > output.sam
  • -e Exclude reads marked as duplicates from discordant, splitter, and/or unmapped file.
  • -d FILE Output discordant read pairs to this file. [no discordant file output]
  • -s FILE Output split reads to this file abiding by paramaters below. [no splitter file output]

  

duplication-readは全出力から除く。

bwa mem -t 12 -R "@RG\tID:X\tLB:Y\tSM:Z\tPL:ILLUMINA" input.fa *.fastq | samblaster -r -e -d samp.disc.sam -s samp.split.sam > output.sam
  •  -r Remove duplicates reads from all output files.  

 

 

 

 

 

 

 

 注意;bwa memで-M(mark shorter split hits as secondary)をつけている時は、samblasterにも-Mをつけてランを行う

  • -M  Compatibility mode (details below); both FLAG 0x100 and 0x800 denote supplemental (chimeric). Similar to bwa mem -M option.

 

 

引用

SAMBLASTER: fast duplicate marking and structural variant read extraction

Gregory G. Faust1 and Ira M. Hall1,2,*

Bioinformatics. 2014 Sep 1; 30(17): 2503–2505.

FastQValidatorでfastqデータを検証する。

 

FastQValidatorは、fastqのフォーマットを検証しておかしなリードが含まれるのか調べることができるツール。具体的には、1つだけファイル名がおかしかったり(ヘッダーが@で始まっていないとか短すぎるとか)、数塩基しかないようなリードが混じっているかどうかなど調べ、エラー数をコールすることができる。使い道はないと思われるかもしれないが、大規模な解析を行っているとfastqのダウンロードが不完全に終わって、気づかないうちに変なfastqを使っているということが起きるかもしれない。FastQValidatorでフォーマットを検証しておけば、そういうトラブルを最小限の手間で防ぐことができる。

本ツールは自動解析ツールの入力ファイルのチェックに使われたりもしている。

 

 

 

インストール

Github

https://github.com/GregoryFaust/samblaster

 

コンパイルにはlibStatGenが必要。本体より先にビルドしておく。

git clone https://github.com/statgen/libStatGen.git 
cd libStatGen/
make all
make test

FastQValidatorのダウンロードとビルド。

git clone https://github.com/statgen/fastQValidator.git 
cd fastQValidator
make all

bin/にfastQValidatorができる。パスを通しておく。

 

ラン

fastqファイルを検証する。

fastQValidator --file input.fastq 

解析が終わるとメッセージが標準出力にプリントされる。

Finished processing input.fastq with 16000000 lines containing 4000000 sequences.

There were a total of 0 errors.

Returning: 0 : FASTQ_SUCCESS

このファイルにはエラーはない。

エラーがあれば以下のようなメッセージがプリントされる。

ERROR on Line 25: The sequence identifier line was too short.
ERROR on Line 29: First line of a sequence does not begin wtih @
ERROR on Line 33: No Sequence Identifier specified before the comment.

 

 

color space fastqかどうか調べる。

fastQValidator --colorSpace --file input.fastq 

ここで使ったfastqはATGCNなので、SOLiDのcolorspace fastq(0123)かどうかチェックすると全行でエラーが出る。

Finished processing input.fastq with 16000000 lines containing 4000000 sequences.

There were a total of 200000000 errors.

Returning: 1 : FASTQ_INVALID

  • --baseSpace ACTGN only
  • --colorSpace 0123. only

 

25-bp以下のリードがないかどうかも調べる。

fastQValidator --minReadLen 25 --file input.fastq
  •  --minReadLen Minimum allowed read length (Defaults to 10).

 

 

塩基の構成も各部位ごとに調べる。

 fastQValidator --baseComposition --file 
  •  --baseComposition Print the Base Composition Statistics.

塩基の割合がポジション別にプリントされる。

Read Index %A %C %G %T %N Total Reads At Index

         0    18.02   27.26   38.16   16.56    0.00 4000000

         1    29.03   23.09   25.11   22.77    0.00 4000000

         2    21.28   25.40   21.87   31.46    0.00 4000000

         3    21.20   32.89   25.26   20.65    0.00 4000000

         4    23.30   22.15   33.77   20.77    0.01 4000000

         5    23.19   23.99   31.71   21.12    0.00 4000000

         6    30.34   23.45   23.68   22.54    0.00 4000000

 

 

複数fastqを入力とし、Phread quality scoreの平均も計算させる。

fastQValidator --avgQual --params --file input*fastq
  • --params Print the parameter settings.
  • --avgQual Print the average phred quality per cycle & overall average quality.

各部位のPhread quality scoreがポジション別にプリントされる。

Average Phred Quality by Read Index (starts at 0):

Read Index Average Quality

0 32.79

1 32.69

2 32.87

3 32.96

4 33.23

5 36.57

6 36.77

 

 

50-bp以下のリードが100以上あればランを止める。

fastQValidator --minReadLen 50 --maxErrors 100 --file input.fastq
  • --maxErrors Number of errors to allow before quitting reading/validating the file. -1 (default) indicates to not quit until the entire file is read.

 

 

全配列をキャッシュするので、大きなfastqだとメモリが足りなくなる恐れがある。心配なら2-passモードでランする。

fastQValidator -2 --file input*fastq #フォルダの全fastqを調べる。

 

引用

https://genome.sph.umich.edu/wiki/FastQValidator

変異のフィルタリングを行うSnpSift

以前SnpEffという変異のアノテーションを行うことができるツールを紹介した(リンク)。このツールにはもう一つSnpSiftというツールが同梱されている。SnpSiftは変異コール結果のVCFファイルを扱うツールで、クオリティやp値など様々な指標に基づいて変異を仕分けることができる。大きなゲノムの解析では数万を超える変異がコールされてくるのが普通だが、そういった時、統計の理屈上偽陽性は一定数出てきてしまう。そこで、はじめに信頼度などの条件に従って変異コール結果をフィルタリングできれば便利である。

SnpSiftはそんな時に使えるツールです。使い方をまとめてみます。

 

公式ページ

SnpSift

インストール

SnpEff

SnpEffのパッケージに入っているので、SnpEffをインストールしている人はすでに持っているはず。SnpEffはbrewでも導入できます。

 

VCFフォーマットについて簡単にまとめています。


 

ラン

filter 条件に従ってフィルタリングする。

 

quality30以上のコールだけ取り出す。

cat variants.vcf | java -jar SnpSift.jar filter " ( QUAL >= 30 )" > filtered.vcf

 

 quality20以上のindelとquality30以上のコールを取り出す。(パイプ"|"でA or Bを表す)

cat variants.vcf | java -jar SnpSift.jar filter "(( exists INDEL ) & (QUAL >= 20)) | (QUAL >= 30 )" > filtered.vcf

 

カバレッジ (DP) 25以上のコールだけ取り出す。

cat variants.vcf | java -jar SnpSift.jar filter "(DP >= 25)" > filtered.vcf

3つ以上のサンプルでhomozygoteのコールだけ取り出す。

cat variants.vcf | java -jar SnpSift.jar filter "(countHom() > 3)" > filtered.vcf

heterozygoteでカバレッジ25以上のコールだけ取り出す。

cat variants.vcf | java -jar SnpSift.jar filter "((countHet() > 0) && (DP >= 25))" > filtered.vcf

 

VCFの標準フォーマットに従っていれば、色々な取り出し方が可能である。例えば下のように

f:id:kazumaxneo:20170819204540j:plain

カラムの説明で#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Zがコールの1行前につけられている。上のVCFからCHROMフィールドが"chr"だけ取り出すなら

cat variants.vcf | java -jar SnpSift.jar filter "( CHROM = 'chr' )" > filtered.vcf

フィルターがPASSしているものだけ取り出す。

cat variants.vcf | java -jar SnpSift.jar filter "(FILTER = 'PASS')"  > filtered.vcf

Positionが 123456-654321のコールだけ取り出す。

cat variants.vcf | java -jar SnpSift.jar filter "( POS > 123456 ) & ( POS < 654321 )" > filtered.vcf

-n オプションを使い、chrとマッチしないものだけ取り出す。

cat variants.vcf | java -jar SnpSift.jar filter -n "( CHROM = 'chr' )" > filtered.vcf
  • -n Inverse. Show lines that do not match filter expression

 

genotypeフィールドでフィルタリングする。

2サンプルのコール結果なら、例えばgenotypeフィールドは

"GT:PL:GQ 1/1:255,66,0:63 0/1:245,0,255:99"

のようになっている。1番サンプルがGQ>60以上の部位だけ取り出すなら

cat variants.vcf | java -jar SnpSift.jar filter "( GEN[0].GQ > 60 )"  > filtered.vcf

サンプルはindexで管理されているので、2番サンプルならGEN[1]。全サンプルならGEN[*]で表せる

 

PLのように複数フィールドを持っているなら(上の例の1番サンプルならPLは255,66,0)、PLの後にも[index]をつける。

cat variants.vcf | java -jar SnpSift.jar filter "( GEN[0].PL[2] = 0 )"  > filtered.vcf

他にも多種多様なフィルタリングが可能になっている。公式に膨大な例があるので、そちらを確認してみてください。

 

 

 

Annotate 別のVCFファイルのアノテーション情報を移す。

 

dbSNPや1000ゲノムなど大規模データのアノテーション情報を移したりするのに使う。

java -jar SnpSift.jar annotate dbSnp132.vcf variants.vcf > variants_annotated.vcf 

 例えばアノテーション情報を移す前のVCFが以下なら、

22      16346045    .            T    C    0.0    FAIL    NS=244

変異が参照先のVCFと合致すると

22      16346045    rs56234788   T    C    0.0    FAIL    NS=244;RSPOS=16346045;GMAF=0.162248628884826;dbSNPBuildID=129;SSR=0;SAO=0;VP=050100000000000100000100;WGT=0;VC=SNV;SLO;GNO

 のようにIDフィールドとINFOフィールドに追記される。-idオプションをつけると、IDフィールドの書き込みだけ行われる。

 

 

他にも様々なコマンドが用意されている。時間を見つけて順番にまとめていきます。

 

CaseControl

 

RmRefGen

 

TsTv

 

 ExtractFields

 

varType

 

gwasCat

 

Split

 

Concordance

 

Private

 

Vcf2Tped

 

Intersect

 

RmInfo

 

GeneSets

 

GT

 

VcfCheck 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

引用

Using Drosophila melanogaster as a model for genotoxic chemical mutational studies with a new program, SnpSift

Cingolani, P., et. al., Frontiers in Genetics, 3, 2012.

 

ユーザー定義の変異を再現可能なfastqのシミュレーター NEAT-genReads

 

NEAT-genReadsは2016年に発表されたfastqをシミュレートできるツール。変異のVCFファイルなどの情報も与えて現実に近いfastqを発生させることができる。fastq以外にポジコンとして使えるbamやVCFファイルも生成されるため、indel検出ツールの妥当性をポジコンと比較して考えることが可能になっている。

NEAT-genReadsにはすでに学習済みのプロファイルが用意されているが、ユーザーがbamやVCFファイルを指定することでユーザー定義のプロファイルでシーケンスbiasやエラー発生率を再現することができる。また領域をbedファイルで明示することで、特定の領域だけシミュレートすることができる(例えばexome解析のシミュレーションに使う)。

 

 

インストール

Github

GitHub - zstephens/neat-genreads: NEAT read simulation tools

本体はpythonで記述されている。ダウンロードしてパスを通す。

 

依存

 

 

ラン

1、デフォルト条件でランすると平均カバレッジ10のfastqが作られる。

./genReads.py -r ref.fa -R 101 -o simulated_data
  • -R read length. Required.
  • -r ref.fa. Required.

ランが終わると"simulated_data_read1.fq" ができる。デフォルトのプロファイルに従い変異は導入されている。

 

 

 

2、インサートサイズが平均600で、カバレッジ x50、301bp x 2の ペアードエンドfastqを作る。アライメントのbamと変異が起きた部位を示したvcfファイルも出力させる。

./genReads.py --pe 600 50 -c 50 -R 301 -r ref.fa -o simulated_data
  • --vcf Output golden VCF file
  • --bam Output golden BAM file
  • --pe Paired-end fragment length mean and standard deviation. To produce paired end data, one of --pe or --pe-model must be specified.
  • -c Average coverage across the entire dataset. Default: 10

"simulated_data_read1.fq" と "simulated_data4_read2.fq" ができる。 

 

bamとvcfをIGVに読み込ませて確認する。(igvがない人はこちらを参照)

samtools index simulated_data_golden.bam #index作成
igv -g ref.fa simulated_data_golden.bam,simulated_data_golden.vcf #IGV起動と読み込み

f:id:kazumaxneo:20170817191004j:plain

一番上が変異部位。ホモとヘテロの部位がある。defaultの条件では かなりの頻度でindel変異が入る。

 

 

 

3、ユーザーが指定したVCFと同じ変異を導入する。また、それ以外の変異の発生率をゼロにしてfastqを作る。

./genReads.py --pe 600 50 -c 50 -R 301 -v mutation_call.vcf -M 0 -r ref.fa -o simulated_data
  • -v Variants from this VCF will be inserted into the simulated sequence with 100% certainty.
  • -M The mutation rate model is rescaled to make this the average value. Must be between 0 and 0.3. These random mutations are inserted in addition to the once specified in the -v option.

 

動作をする場合、使ったvcfと作られたvcfをIGVに読み込んで確認するとわかりやすい。

samtools index simulated_data_golden.bam #index作成
igv -g ref.fa simulated_data_golden.bam,simulated_data_golden.vcf,mutation_call.vcf #IGV起動と読み込み

 

 

ーtでbedを指定すれば、特定の領域だけシミュレートすることができます。exome解析やchip-seqのデータのシミュレートのほか、計算リソースを削減するために特定のクロモソームだけシミュレートしたい時に使えるオプションです。

 

 

 

4、既存データのGCバイアスに基づいてGC biasッモデルを作り、それを使いfastqを作る。

bedtools genomecov -d -ibam normal.bam -g ref.fa > coverage
python2.7 computeGC.py -r ref.fa -i coverage -w 100 -o model
./genReads.py --gc-model GCcoverage -r ref.fa -R 101 -o simulated_data

はじめにbedtoolsでcoverageを計算させている(リンク)。

注意;リード長もモデルのbamと同じと仮定してシミュレートされるため、ラストのコマンドの-R の指定値は元のbamのリード長と揃える必要がある。

 

 

4と似た流れで変異、リードサイズ、シーケンスエラーのbiasを既存データからモデル化して利用することもできます。詳細はGitのHPで確認してください。

GitHub - zstephens/neat-genreads: NEAT read simulation tools

 

 

引用

Simulating Next-Generation Sequencing Datasets from Empirical Mutation and Sequencing Models

Zachary D. Stephens , Matthew E. Hudson, Liudmila S. Mainzer, Morgan Taschuk, Matthew R. Weber, Ravishankar K. Iyer

Published: November 28, 2016https://doi.org/10.1371/journal.pone.0167047

 

 

 

ナノポアのONTリードのシミュレーター NanoSim

 

NanoSImは2017年に発表されたOxford nanoporeのロングリードのシミュレーター。ユーザーが指定したONTリードからプロファイルを作成し、それに基づいてロングリードを発生させることができる。

 

 

インストール

Github

https://github.com/bcgsc/NanoSim

 

依存

  • LAST (Tested with version 581)
  • R (Tested with version 3.2.3)
  • Python (2.6 or above)
  • Numpy (Tested with version 1.10.1 or above)

LASTはbrewで導入できる。

Githubからダウンロードして/srcにパスを通す。

 

 

ラン

ランは二段階で行う。第一ステップはONTのシーケンスデータを指定してのモデルの構築となる。

./read_analysis.py -i ONT.fasta -r reference.fa
  • -r eference genome of the training reads
  • -m User can provide their own alignment file, with maf extension
  • -b number of bins (for development), default = 20
  • -o The prefix of output file, default = 'training'

LASTを使いリファレンスゲノムにONTリードをアライメントしている。エラーを評価するため、ONTリード自身からアセンブルしたcontigをリファレンスに使ったりしてはならない。

カレントディレクトリにref_genome~とtrainning~というファイルがいくつかできる(-o 指定がない時)。

 

 

オーサーらにより、yeastと、E.coliの1dと2dで読んだONTリードのプロファイルやシーケンスデータが用意されている(R7とR9両方あり)。指定のONTリードがないならそれを使う。FTPサーバーリンク

wget ftp://ftp.bcgsc.ca/supplementary/NanoSim/yeast* #例えばyeastのデータをダウンロード

yeast_2D.fasta

yeast_S288C_ref.fa

yeast_profile.zip

がダウンロードされる。yeast_2D.fastaがONTリード、yeast_S288C_ref.faがリファレンスゲノムになる。

 

 

第二段階- 配列のシミュレーション。先ほど作ったtraining~を指定して走らせる。

./simulator.py linear -r referenceg.fa -c training
  •  -r reference genome in fasta file, specify path and file name  
  •  --max_len Maximum read length, default = Inf
  • --min_len Minimum read length, default = 50
  • --perfect Output perfect reads, no mutations, default = False  
  • --KmerBias prohibits homopolymers with length >= 6 bases in output reads, can be omitted  
  • -o The prefix of output file, default = 'simulated' 
  • -n Number of generated reads, default = 20,000 reads  
  • -c the prefix of training set profiles, same as the output prefix in read_analysis.py, default = training  
  • circular | linear Do not choose 'circular' when there is more than one sequence in the reference <options>:

 -oで指定しなければ、simulated_reads.fastaが出力される。環状ゲノムならcircularにすること。

出力された配列をseqkitで簡単に分析する。

user$ seqkit stats simulated_reads.fasta 

file                   format  type  num_seqs      sum_len  min_len  avg_len  max_len

simulated_reads.fasta  FASTA   DNA     20,000  139,782,497      168  6,989.1   25,253

最長25253-bp、最短168-bp、平均6989-bpのリードが出力された。

 

 seqkitは以下で紹介しています。


 

 

引用

NanoSim: nanopore sequence read simulator based on statistical characterization.

Chen Yang, Justin Chu, René L Warren, Inanç Birol

Gigascience 2017 gix010. doi: 10.1093/gigascience/gix010

SNVやindel変異を再現できるfastqのシミュレーターwgsim

 

wgsimはfastqをシミュレートできるツールである。シーケンスエラーを再現したり、diploidゲノムの多型を想定して、一定の確率で変異を入れることができる(indelシーケンスエラーは再現されない)。

wgsimはARTなどのツールでは不可能な300bp以上の配列を発生させることでも可能である。ただしシーケンスプロファイル等はないので、シーケンスのbias等を再現するには向いていない点には注意が必要。論文にはまだなっていない(リンク)。

 

 

インストール

Github

https://github.com/lh3/wgsim

コンパイル

gcc -g -O2 -Wall -o wgsim wgsim.c -lz -lm

 

 

ラン

re.faの配列を元にデフォルトの条件でfastqを発生させる。

wgsim ref.fasta R1.fastq R2.fastq
  • -e base error rate [0.020]
  • -d outer distance between the two ends [500]
  • -s standard deviation [50]
  • -N number of read pairs [1000000]
  • -1 length of the first read [70]
  • -2 length of the second read [70]
  • -r rate of mutations [0.0010]
  • -R fraction of indels [0.15]
  • -X probability an indel is extended [0.30]
  • -S seed for random generator [-1]
  • -A disgard if the fraction of ambiguous bases higher than FLOAT [0.05]
  • -h haplotype mode

カレントディレクトリにR1.fastqとR2.fastqができる。デフォルトだと70-bpでインサートサイズが500、 SDが50、リード数が1000000などとなっている。また、上記のデフォルトに設定された確率でindelやSNV、シーケンスエラーが導入されている。カバレッジはリード数xリード長をゲノムサイズで割って計算できる(ペアエンドはx2)。

f:id:kazumaxneo:20170816201633j:plain

カバレッジにはある程度ばらつきができる。 

 

 

リード長150-bpのfastqを30000リード数発生させる。

wgsim -1 150 -2 150 -N30000 ref.fasta R1.fastq R2.fastq

1-kbp以上のfastqを発生させることも可能であるが、シミュレーションには相当な時間がかかる。ロングリードを作るなら、pacbioかnanoporeのシミュレーターを使うことをお勧めします。

 

 indelやSNVについては以下のパラメータで調整する。

  • -e0 => シーケンスエラーがゼロになる。indelやSNVは起きる。
  • -r0 => indelがゼロになる。シーケンスエラーやSNVは起きる。

 よってindelもSNVもないfastqをシミュレートするならこのようになる。

wgsim -e0 -r0 ref.fasta R1.fastq R2.fastq

 

 

 

引用

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

 

https://popmodels.cancercontrol.cancer.gov/gsr/packages/wgsim/

 

  

NCBIで全データを一度にblast解析し、得られたリストをEntrez Directでアノテーションに変換する。

 

複数の配列のblast解析を行う場合、ローカルでデータベースなどを構築して進めるのが一つの手である。しかしローカルだとデータベースの更新や、データサイズが問題になる(例えばnrのデータも2015年にダウンロードすると200GBを超えていた)。

ネットワーク越しのランで十分な速度が得られる限りは、専門のアノテーターが管理している最新のデータベースを使うのがベストである。そのためにコマンドラインからネットワーク越しにblast解析を行う方法も用意されているが、特定のポートを開ける必要があり、環境によってはやや敷居が高い。今回はNCBIプラウザ越しにアクセスして、数百のcontigのblast解析を行う方法を書いていく。

後半では、全blast結果をテキストでダウンロードし、それからbest hitのIDを取り出し生物名を得る流れをまとめている。ダウンロードされたファイルにはヒットした生物のEntrez_IDしか書かれていないので、生物名を得るためにUNIXのコマンドとNCBIのツールを組み合わせている。

 

生物系研究者の方はSafariなどのプラウザを使ってNCBIにアクセスされている人が多いと思います。その時、全配列を同時にblastにかけて結果をダウンロードする方法を知っていると便利です。一例ですが参考にしてみてください。

 

 

8/20 追記

 

NCBIのNucleotide blastに移動する。

Nucleotide blastを選択。

f:id:kazumaxneo:20170810175809j:plain

 

 

ウィンドウ下のファイルを選択、をクリック。

f:id:kazumaxneo:20170810175844j:plain

 生物を限定するならorganismの項目に学名を書く。

 

de novoアセンブルして得たcontigファイルを選択。

f:id:kazumaxneo:20170810180112j:plain

 BLAST解析を実行する。

 

途中でページ更新に失敗した場合、1ページ戻り右上のRecent resultsを選択(recent resultsリンク)。

f:id:kazumaxneo:20170820135533j:plain

IDをクリックすれば再開できる。

f:id:kazumaxneo:20170820135722j:plain

ジョブが終わっている場合、右端の downloadでダウンロード可能(xmlファイルのみ)。

 

 

 

数分で結果が得られた。

f:id:kazumaxneo:20170810180231j:plain

 

上の方のダウンロードをクリック。

f:id:kazumaxneo:20170810180301j:plain

 Hit Table(text) を選択してダウンロード。

 

ダウンロードされたファイルにはhitした全ての部位の情報がタブ形式で保存されている。

usr $ head alignment.txt

# blastn

# Iteration: 0

# Query: NODE_1_length_40202_cov_39.6477

# RID: SS3HY4K7014

# Database: nr

# Fields: query id, subject ids, query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score

# 264 hits found

NODE_1 gi|545742605|emb|FR873486.1|    NODE_1 FR873486.1      83.535  2976    407     41      9041    11995   12939   15852   0.0     2704

NODE_1 gi|545742605|emb|FR873486.1|    NODE_1 FR873486.1      84.876  2215    295     22      34599   36796   838     3029    0.0     2198

NODE_1 gi|545742605|emb|FR873486.1|    NODE_1 FR873486.1      85.036  1517    197     19      12592   14100   17670   19164   0.0     1517

hitしたIDのところはEntrezのIDとembのIDだが、このままでは何がヒットしたのか分かりにくい。IDをコンバートしてわかりやすく修正する。

 

query_IDとsubject_IDだけ取り出す。

cut -f 1-2 Alignment.txt > subject_and_query_ID

head -15 subject_and_query_ID

# blastn

# Iteration: 0

# Query: NODE_1_length_40202_cov_39.6477

# RID: SS3HY4K7014

# Database: nr

# Fields: query id, subject ids, query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score

# 264 hits found

NODE_1 gi|545742605|emb|FR873486.1|

NODE_1 gi|545742605|emb|FR873486.1|

NODE_1 gi|545742605|emb|FR873486.1|

NODE_1 gi|545742605|emb|FR873486.1|

NODE_1 gi|545742605|emb|FR873486.1|

NODE_1 gi|545742605|emb|FR873486.1|

NODE_1 gi|545742605|emb|FR873486.1|

NODE_1 gi|545742605|emb|FR873486.1|

 

 

1つのcontigに複数当たっている。ここではbest hitだけで十分なので、uniqコマンドを使い1列目のNODE名が重複する行を消す。awkにパイプしてコメント行の削除もワンライナーで行う。

uniq -f 2 subject_and_query_ID | awk '! /^#/' > subject_and_query_ID_unique

 user$ head -5 subject_and_query_ID_unique

NODE_1 gi|545742605|emb|FR873486.1|

NODE_2 gi|332139268|gb|HQ009547.1|

NODE_3 gi|531033897|emb|HF586473.1|

NODE_4 gi|190702157|gb|EF710641.1|

NODE_5 gi|54109776|emb|AJ632317.1|

 

重複が除けたので、ここからsubject_IDだけ保存する。

cut -f 2 subject_and_query_ID_unique > subject_ID_unique

user$ head subject_ID_unique

gi|545742605|emb|FR873486.1|

gi|332139268|gb|HQ009547.1|

gi|531033897|emb|HF586473.1|

gi|190702157|gb|EF710641.1|

gi|54109776|emb|AJ632317.1|

gi|332139193|gb|HQ009535.1|

gi|531033897|emb|HF586473.1|

gi|531033897|emb|HF586473.1|

gi|531034048|emb|HF586474.1|

gi|531033897|emb|HF586473.1|

 

 Entrez IDだけ保存する。

cut -f 2 -d "|" subject_ID_unique > Entrez_ID

user$ head Entrez_ID 

545742605

332139268

531033897

190702157

54109776

332139193

531033897

531033897

531034048

531033897

IDリストが入手できた。

 

ここから先はDAVIDの変換ツールを使っても良いが、今回はEntrez Directのコマンドを使う。ターミナルで以下のコマンドをペーストしてツールをダウンロード。

cd ~
perl -MNet::FTP -e
'$ftp = new Net::FTP("ftp.ncbi.nlm.nih.gov", Passive => 1);
$ftp->login; $ftp->binary;
$ftp->get("/entrez/entrezdirect/edirect.tar.gz");'
gunzip -c edirect.tar.gz | tar xf -
rm edirect.tar.gz
export PATH=$PATH:$HOME/edirect
./edirect/setup.sh

 

ホームにedirect/ができる。パスを通しておく。

echo 'export PATH=$PATH:/Users/user/edirect/' >> ~/.bash_profile
source ~/.bash_profile

 

efetchコマンドでxmlファイルを出力し、xmlからxtractコマンドで<tile>を抽出する。IDを1つだけ探すなら以下のようなコマンドになる。

efetch -id ID -db nucleotide -format docsum | xtract -pattern DocumentSummary -element 'Title' #IDのところはEntrez_IDの番号を入力。

実行例

 

user$ efetch -id 545742605 -db nucleotide -format docsum | xtract -pattern DocumentSummary -element 'Title'

Cotesia congregata bracovirus segment 28, complete sequence

前半の出力は下のようなxmlファイルである。

f:id:kazumaxneo:20170820132227j:plain

Titleの項目が生物名に使えそうである。

f:id:kazumaxneo:20170820132237j:plain

この行のタグの間をEntrez Directのxtractコマンドで抽出している。

 

 

今回はEntrez_IDのリストに対して実行するので、ループで順番に処理させる。

while read line; do title=$(efetch -id $line -db nucleotide -format docsum |  xtract -pattern DocumentSummary -element 'Title'); echo "$line      $title"; done<Entrez_ID > title_list.txt

user$ head -5 output 

545742605      Cotesia congregata bracovirus segment 28, complete sequence

332139268      Cotesia vestalis bracovirus segment c24, complete sequence

531033897      Cotesia congregata sequence containing Cotesia congregata bracovirus proviral locus 2 (PL2)

190702157      Cotesia sesamiae Mombasa bracovirus clone BAC 6L7, complete sequence

54109776      Cotesia congregata virus complete genome, segment Circle14

 IDから<Title>部分の情報に変換できた。

このような流れで、全contigをblastにかけてベストヒットのアノテーションリストを入手できる。

 

集計する。1列目は数値順、2列めは辞書順でソート。

sort -k 1,1n -k 2 title_list.txt |uniq -c |sort -r > ranking.txt

 

591 1228866639 Papilio machaon Ycf2-like protein

7 1228866639 PREDICTED: Helianthus annuus

7 1228808943Medicago truncatula cytochrome b559

2 1044533238 PREDICTED: Vigna angularis LINE-1

.

.

.

省略

 

上のような集計結果を出力できた。使ったcontigはPapilio machaonにbest hitしているcontigが大多数を占める結果となった。

 

 

 

 

今回はリストを入手しbest hitだけ残してblast解析を行ったが、NCBIのblast結果からxmlファイルもダウンロードできる。ケースによってはそちらを使った方が早い。また行程がやや煩雑なので、頻繁に行うならもう少し工夫してスクリプト化した方が良いと思います。

 

blastp、blastx等も流れは同じです。 

 

 

引用

Entrez Direct: E-utilities on the UNIX Command Line

https://www.ncbi.nlm.nih.gov/books/NBK179288/

 

Biostar

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