Transdecoder doesn't integrate Blast and Pfam results

[Alexander Rivero]

Hello, im new to TD and have the following problem running:

TransDecoder.Predict -t target_transcripts.fasta --retain_pfam_hits pfam.domtblout --retain_blastp_hits blastp.outfmt6 

My files for each one have this:

For target_transcripts.fasta (I write the "...." here for make this thread readable)

>contig00001  gene=isogroup00001  length=543  

TCCGCC......tttg

>contig00007  gene=isogroup00001  length=1238  

AGAAGAGCG....AaT

>contig00014  gene=isogroup00001  length=735  

CTTAAA.....TACtCG

For pfam.domtblout

#                                                                                     --- full sequence --- -------------- this domain -------------   hmm coord   ali coord   env coord

# target name        accession   tlen query name                    accession   qlen   E-value  score  bias   #  of  c-Evalue  i-Evalue  score  bias  from    to  from    to  from    to  acc description of target

#------------------- ---------- -----          -------------------- ---------- ----- --------- ------ ----- --- --- --------- --------- ------ ----- ----- ----- ----- ----- ----- ----- ---- ---------------------

PP2C                 PF00481.22   258 Gene.2::contig00007::g.2::m.2 -            272     3e-41  141.8   0.0   1   1   6.1e-45   3.7e-41  141.5   0.0    74   251    60   241    34   247 0.86 Protein phosphatase 2C

PP2C_2               PF13672.7    210 Gene.2::contig00007::g.2::m.2 -            272   0.00047   19.9   0.1   1   1   5.7e-07    0.0035   17.0   0.1    81   185    70   219    27   239 0.66 Protein phosphatase 2C

DUF4298              PF14131.7     87 Gene.2::contig00007::g.2::m.2 -            272     0.034   14.2   0.4   1   2     0.021   1.3e+02    2.7   0.0    24    42    68    86    60    97 0.82 Domain of unknown function (DUF4298)

DUF4298              PF14131.7     87 Gene.2::contig00007::g.2::m.2 -            272     0.034   14.2   0.4   2   2   0.00021       1.3    9.1   0.0    52    79   191   218   186   226 0.82 Domain of unknown function (DUF4298)

SYCP2_SLD            PF18584.2    112 Gene.4::contig00007::g.4::m.4 -            101     0.068   13.4   0.0   1   1   4.7e-06     0.085   13.1   0.0    40    98    14    74     7    81 0.87 Synaptonemal complex 2 Spt16M-like domain

and for blastp.outfmt6

Gene.1::contig00001::g.1::m.1 JI23_HORVU 46.012 163 84 3 16 177 1 160 1.05e-33 121

Gene.2::contig00007::g.2::m.2 P2C34_ARATH 69.027 226 69 1 39 263 125 350 6.26e-115 337

Gene.3::contig00007::g.3::m.3 P2C73_ARATH 68.932 103 32 0 1 103 36 138 9.03e-48 159

Gene.5::contig00014::g.5::m.5 JI23_HORVU 36.548 197 117 5 3 193 5 199 3.24e-26 103

But, in my transdecoder.pep i have this: (I write the "...." here for make this thread readable)

>Gene.1::contig00001::g.1::m.1 Gene.1::contig00001::g.1  ORF type:internal len:181 (-),score=23.07,JI23_HORVU|46.012|1.05e-33 contig00001:3-542(-)

NKN.........GHYAKA

>Gene.2::contig00007::g.2::m.2 Gene.2::contig00007::g.2  ORF type:complete len:272 (+),score=38.67,P2C34_ARATH|69.027|6.26e-115,PP2C|PF00481.22|3.7e-41,PP2C_2|PF13672.7|0.0035,DUF4298|PF14131.7|1.3e+02,DUF4298|PF14131.7|1.3 contig00007:303-1118(+)

MVRGVI........PNN*

>Gene.5::contig00014::g.5::m.5 Gene.5::contig00014::g.5  ORF type:complete len:207 (-),score=15.54,JI23_HORVU|36.548|3.24e-26 contig00014:113-646(-)

MAQVA........GQVN*

The second hit for contig00007 in blastp isn't in transdecoded.pep even i used the retain_blastp_hits as stated at the begining of this thread. What i can do?. Thanks.

[Brian Haas]

Hi,

For the blastp runs, you need to use the long orfs.pep file generated by the first step of transdecoder.  The accessions in the blastp output file should match up with the protein accessions in the long orfs fasta file.

The protein accessions usually have .p number extensions.  

hope this helps,

~b

--
You received this message because you are subscribed to the Google Groups "TransDecoder-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to transdecoder-us...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/transdecoder-users/388192fb-442f-459c-a5e9-fec55460dfe9n%40googlegroups.com.

--

--
Brian J. Haas
The Broad Institute
http://broadinstitute.org/~bhaas

 

[Alexander Rivero]

Hi,

I deleted all files related and made a new run, but still the same problem and the blastp and longorf.pep don't match exactly.

My transcripts file (isotig.fna) is the same from last post:

>contig00001  gene=isogroup00001  length=543  

TCCGCC......tttg

>contig00007  gene=isogroup00001  length=1238  

AGAAGAGCG....AaT

>contig00014  gene=isogroup00001  length=735  

CTTAAA.....TACtCG

My longest_orf.pep look like:

>Gene.1::contig00001::g.1::m.1 type:internal len:181 gc:universal contig00001:542-3(-)

NKNNI......NAGHYAKA

>Gene.2::contig00007::g.2::m.2 type:complete len:272 gc:universal contig00007:303-1118(+)

MVRGVIL..........NLSTNPNN*

>Gene.3::contig00007::g.3::m.3 type:complete len:116 gc:universal contig00007:113-460(+)

MAKEAEKN.......KNKKS*

>Gene.4::contig00007::g.4::m.4 type:complete len:101 gc:universal contig00007:482-180(-)

MLPYIKSFMIFCFYYCNHDQI.....PRQCNNLDLSLYCLFWRKQRQNY*

>Gene.5::contig00014::g.5::m.5 type:5prime_partial len:207 gc:universal contig00014:733-113(-)

SIVIGIPIDRAYLQ........HKKPLQSWEIHESDNMPTKSVDVLDGFSSVATISMGSSALFTAIMSTGQVN*

My blastp.outfmt6 look:

Gene.1::contig00001::g.1::m.1 JI23_HORVU 46.012 163 84 3 16 177 1 160 1.05e-33 121

Gene.2::contig00007::g.2::m.2 P2C34_ARATH 69.027 226 69 1 39 263 125 350 6.26e-115 337

Gene.3::contig00007::g.3::m.3 P2C73_ARATH 68.932 103 32 0 1 103 36 138 9.03e-48 159

Gene.5::contig00014::g.5::m.5 JI23_HORVU 36.548 197 117 5 3 193 5 199 3.24e-26 103

And transdecoder.pep (isotig.fna.transdecoder.pep) show:

>Gene.1::contig00001::g.1::m.1 Gene.1::contig00001::g.1  ORF type:internal len:181 (-),score=23.07,JI23_HORVU|46.012|1.05e-33 contig00001:3-542(-)

NKNNISKY......GHVYEKPYP.VYRGKSSTGVPKDWLQSWEIQEMLGDNMVHTEVRNAGHYAKA

>Gene.2::contig00007::g.2::m.2 Gene.2::contig00007::g.2  ORF type:complete len:272 (+),score=38.67,P2C34_ARATH|69.027|6.26e-115 contig00007:303-1118(+)

MVR.......TNPNN*

>Gene.5::contig00014::g.5::m.5 Gene.5::contig00014::g.5  ORF type:complete len:207 (-),score=15.54,JI23_HORVU|36.548|3.24e-26 contig00014:113-646(-)

MAQVARERQNVEE.......ATISMGSSALFTAIMSTGQVN*

Why it does not show my m.3 from contig00007? Here i attach the run of TransDecoder.Predict if it can help to solve my problem.

When later, i run Trinotate it does not show the hits i want. I think the problem is in transdecoder because i have a old report of trinotate (idk how they did it) where there are the hits i talk about and it have the "prot_coords" like in transdecoder.pep, and when i look all the files i load for the trinotate final report, the only one where are "prot_coords" is ttransdecoder.pep. 

Thanks!

[Brian Haas]

Interesting... . It looks like it's running correctly, but it's not clear to me why that entry would be missing.

Can you try:  

  grep Gene.3::contig00007::g.3::m.3  isotigs.fna.transdecoder_dir/longest_orfs.cds.best_candidates.gff3

and see if it's found there?


If you need to, you can send me (privately) the files  blastp.outfmt6 , isotigs.fna.transdecoder_dir/longest_orfs.cds.scores, and  isotigs.fna.transdecoder_dir/longest_orfs.gff3 and I'll take a look directly.

best,

~b

> To view this discussion on the web visit https://groups.google.com/d/msgid/transdecoder-users/d39e6cd8-7732-44a1-b356-008170be53bbn%40googlegroups.com.

[Alexander Rivero]

I made the grep Gene.3::contig00007::g.3::m.3  isotigs.fna.transdecoder_dir/longest_orfs.cds.best_candidates.gff3 but there isn't results, i cheked the file manually too and it is missing.

I sended a email with the filles attached to your broadinstitute email because i can't write privately from here. Thank you!

[Brian Haas]

Thanks for sending the files.

The issue in this case is that there's a much longer ORF that overlaps this one and also has blastp results.  Because it overlaps by more than 10% of the shorter length of the other orf, the shorter .m3 one is excluded.

best,

~b

To view this discussion on the web visit https://groups.google.com/d/msgid/transdecoder-users/0c5be787-bb93-466d-9d48-b6e0885ac954n%40googlegroups.com.

[Alexander Rivero]

Hi Brian! Thank you so much for your time and help! 

There is a way for overcome that issue and have that ORF included?

Some context: I have a 5 years old Trinotate report where there are the ORFs missing in my report, with the "prot_coord" from the transdecoder.pep. Im trying to mimic it from scratch and following the guides (im learning too, because is my first time doing fuctional annotation). Im starting to work with non-model medicinal and aromatic plants (or herbs) and im highly interested in that ORFs because, as stated in TransDecoder wiki, "To further maximize sensitivity for capturing ORFs that may have functional significance".

[Brian Haas]

Hi Alexander,

I can incorporate a parameter that should allow for it to not exclude hits.  I'll try to do it tonight and get back to you.

best,

~b

To view this discussion on the web visit https://groups.google.com/d/msgid/transdecoder-users/85360ef5-5266-4755-a254-9853e555474bn%40googlegroups.com.

[Brian Haas]

Here's an updated script that you can drop in as a replacement:

https://github.com/TransDecoder/TransDecoder/blob/devel/util/select_best_ORFs_per_transcript.pl

This should truly retain all blastp matching entries, even if they overlap other selected orfs on a given contig.

best,

~brian

[Alexander Rivero]

Thank you! It worked!

Best wishes, Alexander.