macでインフォマティクス

macでインフォマティクス

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

TransDecoderでタンパク質のコード領域を推定する

 

TransDecoderはアセンブリなどで作ったcDNA配列からコード領域を見つけるツール。 RNA seq実験でdo novo assemblyした配列や、cuflinksなどのgenome guide assemblyなツールで作った配列からコード領域を探す時などに使われる。trinityや Trinotateにも取り込まれている。

 

インストー

brewでインストールできる。

brew install TransDecoder

動作に必要なcd-hitもインストールされる。

test.faやユーティリティコマンド群を使うならソースコードもダウンロードする。

https://github.com/TransDecoder/TransDecoder/releases

 

 

マニュアル

https://transdecoder.github.io

 

サポートされているgenetic code (v 3.01)

universal (default)

Euplotes

Tetrahymena

Candida

Acetabularia

Mitochondrial-Canonical

Mitochondrial-Vertebrates

Mitochondrial-Arthropods

Mitochondrial-Echinoderms

Mitochondrial-Molluscs

Mitochondrial-Ascidians

Mitochondrial-Nematodes

Mitochondrial-Platyhelminths

Mitochondrial-Yeasts

Mitochondrial-Euascomycetes

Mitochondrial-Protozoans

デフォルトではuniversalになっており、-Gで変更できる。

 

 

ラン

 

1、de novo assemblyで作ったfastaを使う。

sample-dataのTrinityのアセンブルデータからコード領域を探す。(TransDecoder-3.0.1/sample_data/simple_transcriptome_target/Trinity.fasta)。

step1

コード領域を探索する。

TransDecoder.LongOrfs -G universal -m 100 -t Trinity.fasta
  • -t transcripts.fasta
  • -m minimum protein length (default: 100)
  • -S strand-specific (only analyzes top strand)
  • -G genetic code (default: universal)

デフォルトでは100アミノ酸以上のORFが抽出される。strand specific sequenceなどでリード鎖が決まっている場合-Sをつける。

解析が終わると、Trinity.fasta.transdecoder_dir/ができる。中身は以下のようになった。

f:id:kazumaxneo:20170730101335j:plain

100アミノ酸以上になるcode領域を記載したgff3ファイルができる。またそのアミノ酸配列ファイル(.pep)もできている。

 

step2 (option)

見つかったコード領域をクエリにしてデータベースを検索。

 

blastxで相同なタンパク質を検索。(公式マニュアルの例)

blastp -query transdecoder_dir/longest_orfs.pep  -db uniprot_sprot.fasta  -max_target_seqs 1 -outfmt 6 -evalue 1e-5 -num_threads 10 > blastp.outfmt6

Pfamで相同なタンパク質を検索。(公式マニュアルの例)

hmmscan --cpu 8 --domtblout pfam.domtblout /path/to/Pfam-A.hmm transdecoder_dir/longest_orfs.pep

 

step3

blastやpfamの結果をもとに妥当なcode領域を選抜する。

TransDecoder.Predict -t target_transcripts.fasta --retain_blastp_hits blastp.outfmt6
  • -t   transcripts.fasta
  • -G    genetic code (default: universal)
  • --retain_pfam_hits domain table output file from running hmmscan to search Pfam.
  • --retain_blastp_hits blastp output in '-outfmt 6' format.
  • --single_best_orf Retain only the single best ORF per transcript. (Best is defined as having (optionally pfam and/or blast support) and longest orf)

 

Trinity.fasta.transdecoder_dir/にさらに中間ファイルが作成され、最終的にカレントディレクトリに上記の条件を元に選ばれたcode領域と除去されたcode領域のスコアを記載したファイルが出力される。

f:id:kazumaxneo:20170730101923j:plain

一番上に.bed、.gff3、.pep(fasta)、.cds(fasta)ができている。

gff3とbedには最終的に選ばれた領域を記載されている。また最終的に選ばれたコード領域のアミノ酸配列は.pepに記載され、コード領域の塩基配列は.cdsに記載され出力される。

この最終段階では選ばれなかったコード領域のアミノ酸配列は完全に除去されている。内部でcd-hitを動かしていると推測される。

igvで結果を確認する。

igv -g Trinity.fasta Trinity.fasta.transdecoder.bed

f:id:kazumaxneo:20170730102445j:plain

 太い部分が予想されたコード領域になる。

 

 

2、genome guide assemblyで作った gtfファイルと参照したgenome配列を使う。

sample-dataのcufflinksのゲノムガイドアセンブルで作られたgtfからコード領域を探す。(TransDecoder-3.0.1/sample_data/cufflinks_example/)。

step1

gtfで指定した領域をfastaにして抽出する。

util/cufflinks_gtf_genome_to_cdna_fasta.pl transcripts.gtf test.genome.fasta > transcripts.fasta 

transcripts.fastaができる。このコマンドはgtfからfasta配列を作るツールとしても使えそうである。

 

step2

gff3も作る。最後に必要なためで、それ以外にgff3を作る意味はない。

util/cufflinks_gtf_to_alignment_gff3.pl transcripts.gtf > transcripts.gff3

このスクリプトもgtfからgff3フォーマット変換するのに使えそうである。

 

step3

コード領域を探索する。de novo で始めた場合のstep1にあたる。

TransDecoder.LongOrfs -t transcripts.fasta

解析が終わると、transcripts.fasta.transdecoder_dir/ができる。

f:id:kazumaxneo:20170730103822j:plain

 

step4

妥当なcode領域を選抜する。de novo で始めた場合のstep3にあたる。

TransDecoder.Predict -t transcripts.fasta

カレントディレクトリに.gff3、.bed、.fasta、.cdsファイルが出力される。

f:id:kazumaxneo:20170730103927j:plain

 

step5

最後にゲノムベースのコード領域のアノテーションファイル(gff3)を作る。

util/cdna_alignment_orf_to_genome_orf.pl transcripts.fasta.transdecoder.gff3 transcripts.gff3 transcripts.fasta > transcripts.fasta.transdecoder.genome.gff3

f:id:kazumaxneo:20170730104117j:plain

 

igvで結果を確認する。

igv -g test.genome.fasta transcripts.gtf,transcripts.fasta.transdecoder.genome.gff3

f:id:kazumaxneo:20170730104757j:plain

 上がcufflinksのmergeで作ったgtfファイル(transcripts.gtf)、下がコード領域を予想して最終的に得られたtranscripts.fasta.transdecoder.genome.gff3。太い部分が予想されたコード領域である。

 

 

 

igvについては以前紹介しています。

 

cufflinksも以前何度か取り上げています。例えば以下にcufflinksを使ったワークフローを書いています。

 

引用

De novo transcript sequence reconstruction from RNA-Seq: reference generation and analysis with Trinity

Brian J. Haas,#1,# Alexie Papanicolaou,#2 Moran Yassour,1,3 Manfred Grabherr,4 Philip D. Blood,5 Joshua Bowden,6 Matthew Brian Couger,7 David Eccles,8 Bo Li,9 Matthias Lieber,10 Matthew D. MacManes,11 Michael Ott,2 Joshua Orvis,12 Nathalie Pochet,1,13 Francesco Strozzi,14 Nathan Weeks,15 Rick Westerman,16 Thomas William,17 Colin N. Dewey,9,18 Robert Henschel,19 Richard D. LeDuc,19 Nir Friedman,3 and Aviv Regev1,20,#

Nat Protoc. 2013 Aug; 8(8): 10.1038/nprot.2013.084.