Pythonバイオ?　Pythonバイオ/ツール?
1855　　　2019-08-12 (月) 15:06:21

AP012030

AP012030用のsnpEffの設定 †

想定：　Breseqの出力にVCNファイルがある。それをsnpEffに食わせる。

snpEffのマニュアル

↑

snpEffのデータベースのエントリーを作成 †

（アノテーション）データベースは、手持ちは

java -jar /usr/local/snpEff.jar databases | grep -i Escherichia_coli

で探せる。

k-12株はある（Escherichia_coli_k_12）が、chromosome欄が違う。実際、snpEffで

java -jar /usr/local/snpEff.jar Escherichia_coli_k_12 output.vcf > 1_output.ann.vcf

とすると、うまくゆかない（ほとんど空状態）

データベースの追加方法は、入力データとしてGenbank、GTFやGFF3があるらしい。

genbankがもっともオリジナルデータなので、これを使ってみる。ただしFASTA情報も必要。

まず、データディレクトリ ~/src/snpEff/data/Escherichia_coli_k_12_AP012030 を作成する。この下に

AP012030.gb -> genes.gbk　として置く（名前は決め打ち？）
AP012030.fasta -> sequences.fa　として置く （同上）

次に、~/src/snpEff/ ディレクトリ下のファイル　snpEff.config の、それらしい場所（アルファベットでソートされているので）に、

Escherichia_coli_k_12_AP012030.genome : Escherichia_coli_k_12_AP012030

の1行を追加する。

その上で、

java -jar /usr/local/snpEff.jar build -genbank -v Escherichia_coli_k_12_AP012030

結果は、

00:00:00	SnpEff version SnpEff 4.3t (build 2017-11-24 10:18), by Pablo Cingolani
00:00:00	Command: 'build'
00:00:00	Building database for 'Escherichia_coli_k_12_AP012030'
00:00:00	Reading configuration file 'snpEff.config'. Genome: 'Escherichia_coli_k_12_AP012030'
00:00:00	Reading config file: /home/yamanouc/src/snpEff/snpEff.config
00:00:00	done
Chromosome: 'AP012030'	length: 4621430

	Create exons from CDS (if needed): 
...............（中略）.
	Exons created for 4260 transcripts.

	Deleting redundant exons (if needed): 
		Total transcripts with deleted exons: 0

	Collapsing zero length introns (if needed): 
		Total collapsed transcripts: 0
		Adding genomic sequences to exons: 	Done (4261 sequences added, 0 ignored).

	Adjusting transcripts: 
	Adjusting genes: ................
	Adjusting chromosomes lengths: 
	Ranking exons: 
	Create UTRs from CDS (if needed): 
	Remove empty chromosomes: 

	Marking as 'coding' from CDS information: 
	Done: 4260 transcripts marked
00:00:01	Caracterizing exons by splicing (stage 1) : 
	....
00:00:01	Caracterizing exons by splicing (stage 2) : 
	....00:00:01	done.
00:00:01	[Optional] Rare amino acid annotations
00:00:01	Warning: Cannot read optional protein sequence file '/home/yamanouc/src/snpEff/./data/Escherichia_coli_k_12_AP012030/protein.fa', nothing done.
00:00:01	Protein check file: '/home/yamanouc/src/snpEff/./data/Escherichia_coli_k_12_AP012030/genes.gbk'

00:00:01	Checking database using protein sequences
00:00:01	Comparing Proteins...
	Labels:
		'+' : OK
		'.' : Missing
		'*' : Error
	. 
 .......（中略）.......................................

	Protein check:	Escherichia_coli_k_12_AP012030	OK: 0	Not found: 4260	 Errors: 0	Error percentage: NaN%
00:00:01	Saving database
00:00:02	[Optional] Reading regulation elements: GFF
00:00:02	Warning: Cannot read optional regulation file '/home/yamanouc/src/snpEff/./data/Escherichia_coli_k_12_AP012030/regulation.gff', nothing done.
00:00:02	[Optional] Reading regulation elements: BED 
00:00:02	Cannot find optional regulation dir '/home/yamanouc/src/snpEff/./data/Escherichia_coli_k_12_AP012030/regulation.bed/', nothing done.
00:00:02	[Optional] Reading motifs: GFF
00:00:02	Warning: Cannot open PWMs file /home/yamanouc/src/snpEff/./data/Escherichia_coli_k_12_AP012030/pwms.bin. Nothing done
00:00:02	Done
00:00:02	Logging

ということで、データベース作成は完了。

src/snpEff/data/Escherichia_coli_k_12_AP012030 下に snpEffectPredictor.bin が生成された。

↑

VCFファイルの処理 †

それで、VCFファイルの処理を試みる。

たとえば、~/KishimotoRNA2-orig/breseq/01_43B_S1/data ディレクトリ下で

java -jar /usr/local/snpEff.jar Escherichia_coli_k_12_AP012030 output.vcf > 1_output.ann.vcf

を実行。出力として、��snpEff_genes.txt、��snpEff\summary.html、��1_output.ann.vcf を得た。

出力（stats）ファイル名を変えるには、-s <filename> で指定。Statsデータをcsvで出力するには、-csvStats <filename>

java -jar /usr/local/snpEff.jar -s <sample>.html -csvStats <sample>.csv Escherichia_coli_k_12_AP012030 <sample>.vcf > <sample>.ann.vcf

同様に、~/1710JNHX-0008/rawdata/2_5-1/AP012030/data ディレクトリ下で

java -jar /usr/local/snpEff.jar Escherichia_coli_k_12_AP012030 output.vcf > output.ann.vcf

を実行。

結果は、http://pepper.is.sci.toho-u.ac.jp/Kishimoto2017/SnpEff/　
　なお、breseqの出力は　http://pepper.is.sci.toho-u.ac.jp/Kishimoto2017/breseq/

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2019-08-12 (月) 15:06:21 (1325d)