<p>If you need to model introns and splicing, then you should use a tool that is designed for that. One of the limitations of blast and its cousins is that they only look at sequence similarity and do not attempt to join <a href="https://www.ncbi.nlm.nih.gov/books/NBK62051/" rel="nofollow noreferrer">HSPs</a> into possible transcript models.</p>
<p>I haven't worked on this in a long time, but back when I did, the programs I used were:</p>
<ul>
<li><p><a href="https://www.ebi.ac.uk/about/vertebrate-genomics/software/exonerate" rel="nofollow noreferrer">exonerate</a></p>
<blockquote>
<p>Exonerate is a generic tool for pairwise sequence comparison. It allows you to align sequences using a many alignment models, either exhaustive dynamic programming or a variety of heuristics.</p>
</blockquote>
<p>In this case, you would use your protein sequence as the query and the set of  DNA sequences as the target and try to align using the protein2genome model. Something like this:</p>
<pre><code>  exonerate --model protein2genome protein.pep dna.fa &gt; out
</code></pre>
<p>This will result in (many) alignments of your query protein agains the various genomic sequences. Here's an example I had handy:</p>
<pre><code>$ exonerate --model protein2genome -n 1 frog.histone.pep c.elegans.chrI.fa 
Command line: [exonerate --model protein2genome -n 1 frog.histone.pep c.elegans.chrI.fa]
Hostname: [oregano]

C4 Alignment:
------------
         Query: frog
        Target: CM041006.1 Caenorhabditis elegans strain DL226 chromosome I, whole genome shotgun sequence:[revcomp]
         Model: protein2genome:local
     Raw score: 732
   Query range: 32 -&gt; 300
  Target range: 3039620 -&gt; 3034073

      33 : LeuGlnLysLysLeuGluGluLeuGluLeuAspGluGlnGlnArgLysArgLeuGl :      51
           |||:!!:!!:!!|||:!!:!!|||! !|||   |||! !|||! !:!!||||||..
           LeuGluGluGlnLeuLysLysLeuGlyLeuThrGluProGlnThrGlnArgLeuSe
 3039620 : CTAGAAGAACAGCTGAAAAAGCTCGGGCTCACGGAGCCACAAACTCAACGGCTCAG : 3039566

      52 : uAlaPheLeuThrGlnLysGlnLysValGlyGluLeuLysAspAspAspPheGluL :      70
           .! !||||||  !   |||:!!   :!!  !||||||!..!!:|||   .!.:!!!
           rGluPheLeuGlnValLysGluGlyIleLysGluLeuSerGluAspMetLeuGlnT
 3039565 : CGAATTCTTGCAGGTCAAGGAGGGCATCAAGGAGTTGAGTGAGGATATGCTACAGA : 3039509

      71 : ysValSerGluLeuGlyAlaGlyAsnGlyGlyValValPheLys  &gt;&gt;&gt;&gt; Targe :      85
             ! !.!!|||||||||  !||||||||||||||||||  !|||            
           hrGluGlyGluLeuGlyHisGlyAsnGlyGlyValValAsnLys++          
 3039508 : CTGAAGGAGAACTTGGCCATGGAAATGGAGGAGTTGTTAATAAGgt.......... : 3039462

      86 : t Intron 1 &gt;&gt;&gt;&gt;  ValSerHisLysProThrSerLeuIleMetAlaArgLys :      97
           55 bp              !  !|||::!  !|||.!!:!!|||||||||||||||
                          ++CysValHisArgLysThrGlyValIleMetAlaArgLys
 3039461 : ...............agTGTGTTCATCGAAAAACCGGTGTAATTATGGCTCGGAAG : 3039373

      98 : LeuIleHisLeuGluIleLysProAlaIleArgAsnGlnIleIleArgGluLeuGl :     116
           |||:!!||||||||||||||||||:!!:!!|||.!.||||||:!!!:!||||||  
           LeuValHisLeuGluIleLysProSerValArgGlnGlnIleValLysGluLeuAl
 3039372 : CTTGTACACCTCGAGATCAAGCCGTCAGTCCGACAACAGATCGTCAAAGAGCTCGC : 3039316

     117 : nValLeuHisGluCysAsnSerProTyrIleValGlyPheTyrGlyAlaPheTyrS :     135
            |||||||||:!!||||||||||||!:!||||||||||||||||||||||||   .
           aValLeuHisLysCysAsnSerProPheIleValGlyPheTyrGlyAlaPheValA
 3039315 : CGTACTTCACAAGTGCAATAGTCCATTCATTGTTGGATTCTACGGGGCTTTTGTGG : 3039259

     136 : erAspGlyGluIleSerIleCysMetGluHisMetAspGlyGlySerLeuAspGln :     153
           .!:!!..!!!:||||||||||||||||||:!!|||||||||  !|||||||||  !
           spAsnAsnAspIleSerIleCysMetGluTyrMetAspGlyLeuSerLeuAspIle
 3039258 : ACAATAACGATATATCAATTTGTATGGAGTATATGGATGGGCTTTCACTTGATATA : 3039205

     154 : ValLeuLysLys  &gt;&gt;&gt;&gt; Target Intron 2 &gt;&gt;&gt;&gt;  AlaGlyLysIlePro :     162
           ||||||||||||           3837 bp           !.!|||::::!!|||
           ValLeuLysLys++                         ++ValGlyArgLeuPro
 3039204 : GTTTTGAAAAAGgt.........................agGTTGGCCGCCTCCCG : 3035341

     163 : GluLysIleLeuGlyLysValSerIleAlaValIleLysGlyLeuThrTyrLeuAr :     181
           ||||||.!!:!!|||::!:!!|||:!!||||||:!!:::|||||||||||||||!:
           GluLysPheValGlyArgIleSerValAlaValValArgGlyLeuThrTyrLeuLy
 3035340 : GAAAAGTTTGTCGGACGAATCTCGGTTGCAGTTGTCCGTGGATTGACGTACTTGAA : 3035284

     182 : gGluLysHisLysIleMetHisArg{A}  &gt;&gt;&gt;&gt; Target Intron 3 &gt;&gt;&gt;&gt;  :     190
           !!!::!!  !||||||:!!||||||{|}            50 bp           
           sAspGluIleLysIleLeuHisArg{A}++                         +
 3035283 : AGATGAAATTAAGATTTTGCATAGA{G}gt.........................a : 3035208

     191 :  {sp}ValLysProSerAsnIleLeuValAsnSerArgGlyGluIle--LysLeuC :     207
            {||}|||||||||||||||!!:||||||||||||!.!|||||||||##   ||||
           +{sp}ValLysProSerAsnMetLeuValAsnSerAsnGlyGluIle##---LeuC
 3035207 : g{AT}GTTAAACCATCGAATATGCTTGTCAATAGTAACGGAGAAATTAA---TTAT : 3035157

     208 : ysAspPheGlyValSerGlyGlnLeuIleAspSerMetAlaAsnSerPheValGly :     225
           ||||||||||||||||||||..!|||||||||||||||||||||||||||||||||
           ysAspPheGlyValSerGlyMetLeuIleAspSerMetAlaAsnSerPheValGly
 3035156 : GCGATTTTGGAGTCTCTGGAATGTTGATTGATTCGATGGCCAACTCATTTGTCGGA : 3035103

     226 : ThrArgSerTyrMetSer  &gt;&gt;&gt;&gt; Target Intron 4 &gt;&gt;&gt;&gt;  ProGluArg :     234
           |||||||||||||||:!!            799 bp           |||||||||
           ThrArgSerTyrMetAla++                         ++ProGluArg
 3035102 : ACGAGGAGTTATATGGCGgt.........................agCCCGAACGA : 3034277

     235 : LeuGlnGlyThrHisTyrSerValGlnSerAspIleTrpSerMetGlyLeuSerLe :     253
           |||  !|||:!!||||||!!!:!!..!|||||||||||||||.!.|||||||||||
           LeuThrGlySerHisTyrThrIleSerSerAspIleTrpSerPheGlyLeuSerLe
 3034276 : CTCACAGGATCCCACTATACAATTTCATCGGATATCTGGTCATTTGGATTATCACT : 3034220

     254 : uValGluMetAlaIleGlyArgTyrProIleProProProAspAlaLysGluLeuG :     272
           |||||||:!:  !|||||||||||||||:!:||| !!|||..!  !  !|||! !!
           uValGluLeuLeuIleGlyArgTyrProValProAlaProSerGlnAlaGluTyrA
 3034219 : TGTCGAACTTCTTATTGGACGATATCCAGTGCCTGCACCGTCTCAAGCAGAATATG : 3034163

     273 : luLeuIlePheGlyCysSerVal&lt;-&gt;GluGlyAspProAlaSerSerGluLeuAla :     289
               !!!:|||..!   :!!! !   |||  !!!:! !|||..!|||     !! !
           laThrMetPheAsnValAlaGluAsnGluIleGluLeuAlaAspSerLeuGluGlu
 3034162 : CTACAATGTTTAATGTGGCTGAAAATGAGATTGAACTGGCTGATAGTCTTGAAGAG : 3034109

     290 : ProArgProArgProProGlyArgProIleSer :     300
           |||!.!  !!.!||||||.!!!.!|||   |||
           ProAsnTyrHisProProSerAsnProAlaSer
 3034108 : CCCAACTATCATCCGCCGAGCAATCCGGCGTCA : 3034074

vulgar: frog 32 300 . CM041006.1 3039620 3034073 - 732 M 52 156 5 0 2 I 0 51 3 0 2 M 73 219 5 0 2 I 0 3833 3 0 2 M 32 96 S 0 1 5 0 2 I 0 46 3 0 2 S 1 2 M 14 42 F 0 2 G 1 0 M 26 78 5 0 2 I 0 795 3 0 2 M 48 144 G 0 3 M 21 63
-- completed exonerate analysis
</code></pre>
<p>In addition, exonerate lets you control the output you want. Check out the <code>--ryo</code> (roll your own output) option in <code>man exonerate</code>.</p>
</li>
<li><p><a href="https://www.ebi.ac.uk/Tools/psa/genewise/" rel="nofollow noreferrer">genewise</a></p>
<blockquote>
<p>GeneWise compares a protein sequence to a genomic DNA sequence, allowing for introns and frameshifting errors.</p>
</blockquote>
<p>Unlike exonerate, genewise <em>only</em> does protein to DNA alignment, but like exonerate it models introns and splicing. 20 years ago, when I did this every day, my impression was that genewise was slower but slightly more sensitive than the default exonerate options (exonerate can be tweaked to be made more sensitive). So I would usually try exonerate first and fall back to genewise if it failed to find a good match.</p>
</li>
</ul>

CC BY-SA 4.0
2023-10-13T12:03:12.847Z
2023-10-13T12:03:12.847Z
terdon
<p>I have a dataset with lots of DNA sequences (~10 M seqs, each ~6 KBp long).
I want to detect presence of a fixed given protein (~600 Bp / 200 aa) within each sequence and, if detected, obtain the frame and position of the protein within the sequence. However, the protein may be coded arbitrarily and differently in each sequence. In addition, I need to account for frameshifts induced by sequencing errors.</p>
<p>I've tried to use <code>tfasty36</code> from the <a href="https://github.com/wrpearson/fasta36" rel="nofollow noreferrer">FASTA36</a> package. This is basically a <code>tblastn</code> equivalent (protein query search against DNA database in all frames). Due to the statistical filtering built into these tools, I've come to the conclusion that <code>tfasty36</code> and <code>tblastn</code> are not the right tools for this job.</p>
<p><code>tfasty36</code> didn't work because it requires fiddling with statistical tolerances to show all matches. Also, the matches are quite fragmented (I guess, it uses heuristics instead of dynamic programming). For <code>tblastn</code>, building the database seems to be unnecessary (given that it may have the same problems as <code>tfasty36</code> afterwards).</p>
<p>I've looked into alignment tools, but couldn't find any that can align a protein sequence to a DNA sequence while taking frameshifts into account.</p>
<p>How can I solve this problem?</p>

<phylogenetics>
Find protein in DNA sequences with arbitrary encoding and possible frameshifts
# Find protein in DNA sequences with arbitrary encoding and possible frameshifts

(Asked by: username_0 on 2023-10-12)

I have a dataset with lots of DNA sequences (~10 M seqs, each ~6 KBp long).
I want to detect presence of a fixed given protein (~600 Bp / 200 aa) within each sequence and, if detected, obtain the frame and position of the protein within the sequence. However, the protein may be coded arbitrarily and differently in each sequence. In addition, I need to account for frameshifts induced by sequencing errors.

I've tried to use `tfasty36` from the [FASTA36](https://github.com/wrpearson/fasta36) package. This is basically a `tblastn` equivalent (protein query search against DNA database in all frames). Due to the statistical filtering built into these tools, I've come to the conclusion that `tfasty36` and `tblastn` are not the right tools for this job.

`tfasty36` didn't work because it requires fiddling with statistical tolerances to show all matches. Also, the matches are quite fragmented (I guess, it uses heuristics instead of dynamic programming). For `tblastn`, building the database seems to be unnecessary (given that it may have the same problems as `tfasty36` afterwards).

I've looked into alignment tools, but couldn't find any that can align a protein sequence to a DNA sequence while taking frameshifts into account.

How can I solve this problem?

**Tags:** phylogenetics

---

## 1 answers

### Answer by username_1 on 2023-10-13. Score: 2 | Accepted answer

If you need to model introns and splicing, then you should use a tool that is designed for that. One of the limitations of blast and its cousins is that they only look at sequence similarity and do not attempt to join HSPs into possible transcript models.
I haven't worked on this in a long time, but back when I did, the programs I used were:

exonerate

Exonerate is a generic tool for pairwise sequence comparison. It allows you to align sequences using a many alignment models, either exhaustive dynamic programming or a variety of heuristics.

In this case, you would use your protein sequence as the query and the set of  DNA sequences as the target and try to align using the protein2genome model. Something like this:
  exonerate --model protein2genome protein.pep dna.fa > out

This will result in (many) alignments of your query protein agains the various genomic sequences. Here's an example I had handy:
$ exonerate --model protein2genome -n 1 frog.histone.pep c.elegans.chrI.fa 
Command line: [exonerate --model protein2genome -n 1 frog.histone.pep c.elegans.chrI.fa]
Hostname: [oregano]

C4 Alignment:
------------
         Query: frog
        Target: CM041006.1 Caenorhabditis elegans strain DL226 chromosome I, whole genome shotgun sequence:[revcomp]
         Model: protein2genome:local
     Raw score: 732
   Query range: 32 -> 300
  Target range: 3039620 -> 3034073

      33 : LeuGlnLysLysLeuGluGluLeuGluLeuAspGluGlnGlnArgLysArgLeuGl :      51
           |||:!!:!!:!!|||:!!:!!|||! !|||   |||! !|||! !:!!||||||..
           LeuGluGluGlnLeuLysLysLeuGlyLeuThrGluProGlnThrGlnArgLeuSe
 3039620 : CTAGAAGAACAGCTGAAAAAGCTCGGGCTCACGGAGCCACAAACTCAACGGCTCAG : 3039566

      52 : uAlaPheLeuThrGlnLysGlnLysValGlyGluLeuLysAspAspAspPheGluL :      70
           .! !||||||  !   |||:!!   :!!  !||||||!..!!:|||   .!.:!!!
           rGluPheLeuGlnValLysGluGlyIleLysGluLeuSerGluAspMetLeuGlnT
 3039565 : CGAATTCTTGCAGGTCAAGGAGGGCATCAAGGAGTTGAGTGAGGATATGCTACAGA : 3039509

      71 : ysValSerGluLeuGlyAlaGlyAsnGlyGlyValValPheLys  >>>> Targe :      85
             ! !.!!|||||||||  !||||||||||||||||||  !|||            
           hrGluGlyGluLeuGlyHisGlyAsnGlyGlyValValAsnLys++          
 3039508 : CTGAAGGAGAACTTGGCCATGGAAATGGAGGAGTTGTTAATAAGgt.......... : 3039462

      86 : t Intron 1 >>>>  ValSerHisLysProThrSerLeuIleMetAlaArgLys :      97
           55 bp              !  !|||::!  !|||.!!:!!|||||||||||||||
                          ++CysValHisArgLysThrGlyValIleMetAlaArgLys
 3039461 : ...............agTGTGTTCATCGAAAAACCGGTGTAATTATGGCTCGGAAG : 3039373

      98 : LeuIleHisLeuGluIleLysProAlaIleArgAsnGlnIleIleArgGluLeuGl :     116
           |||:!!||||||||||||||||||:!!:!!|||.!.||||||:!!!:!||||||  
           LeuValHisLeuGluIleLysProSerValArgGlnGlnIleValLysGluLeuAl
 3039372 : CTTGTACACCTCGAGATCAAGCCGTCAGTCCGACAACAGATCGTCAAAGAGCTCGC : 3039316

     117 : nValLeuHisGluCysAsnSerProTyrIleValGlyPheTyrGlyAlaPheTyrS :     135
            |||||||||:!!||||||||||||!:!||||||||||||||||||||||||   .
           aValLeuHisLysCysAsnSerProPheIleValGlyPheTyrGlyAlaPheValA
 3039315 : CGTACTTCACAAGTGCAATAGTCCATTCATTGTTGGATTCTACGGGGCTTTTGTGG : 3039259

     136 : erAspGlyGluIleSerIleCysMetGluHisMetAspGlyGlySerLeuAspGln :     153
           .!:!!..!!!:||||||||||||||||||:!!|||||||||  !|||||||||  !
           spAsnAsnAspIleSerIleCysMetGluTyrMetAspGlyLeuSerLeuAspIle
 3039258 : ACAATAACGATATATCAATTTGTATGGAGTATATGGATGGGCTTTCACTTGATATA : 3039205

     154 : ValLeuLysLys  >>>> Target Intron 2 >>>>  AlaGlyLysIlePro :     162
           ||||||||||||           3837 bp           !.!|||::::!!|||
           ValLeuLysLys++                         ++ValGlyArgLeuPro
 3039204 : GTTTTGAAAAAGgt.........................agGTTGGCCGCCTCCCG : 3035341

     163 : GluLysIleLeuGlyLysValSerIleAlaValIleLysGlyLeuThrTyrLeuAr :     181
           ||||||.!!:!!|||::!:!!|||:!!||||||:!!:::|||||||||||||||!:
           GluLysPheValGlyArgIleSerValAlaValValArgGlyLeuThrTyrLeuLy
 3035340 : GAAAAGTTTGTCGGACGAATCTCGGTTGCAGTTGTCCGTGGATTGACGTACTTGAA : 3035284

     182 : gGluLysHisLysIleMetHisArg{A}  >>>> Target Intron 3 >>>>  :     190
           !!!::!!  !||||||:!!||||||{|}            50 bp           
           sAspGluIleLysIleLeuHisArg{A}++                         +
 3035283 : AGATGAAATTAAGATTTTGCATAGA{G}gt.........................a : 3035208

     191 :  {sp}ValLysProSerAsnIleLeuValAsnSerArgGlyGluIle--LysLeuC :     207
            {||}|||||||||||||||!!:||||||||||||!.!|||||||||##   ||||
           +{sp}ValLysProSerAsnMetLeuValAsnSerAsnGlyGluIle##---LeuC
 3035207 : g{AT}GTTAAACCATCGAATATGCTTGTCAATAGTAACGGAGAAATTAA---TTAT : 3035157

     208 : ysAspPheGlyValSerGlyGlnLeuIleAspSerMetAlaAsnSerPheValGly :     225
           ||||||||||||||||||||..!|||||||||||||||||||||||||||||||||
           ysAspPheGlyValSerGlyMetLeuIleAspSerMetAlaAsnSerPheValGly
 3035156 : GCGATTTTGGAGTCTCTGGAATGTTGATTGATTCGATGGCCAACTCATTTGTCGGA : 3035103

     226 : ThrArgSerTyrMetSer  >>>> Target Intron 4 >>>>  ProGluArg :     234
           |||||||||||||||:!!            799 bp           |||||||||
           ThrArgSerTyrMetAla++                         ++ProGluArg
 3035102 : ACGAGGAGTTATATGGCGgt.........................agCCCGAACGA : 3034277

     235 : LeuGlnGlyThrHisTyrSerValGlnSerAspIleTrpSerMetGlyLeuSerLe :     253
           |||  !|||:!!||||||!!!:!!..!|||||||||||||||.!.|||||||||||
           LeuThrGlySerHisTyrThrIleSerSerAspIleTrpSerPheGlyLeuSerLe
 3034276 : CTCACAGGATCCCACTATACAATTTCATCGGATATCTGGTCATTTGGATTATCACT : 3034220

     254 : uValGluMetAlaIleGlyArgTyrProIleProProProAspAlaLysGluLeuG :     272
           |||||||:!:  !|||||||||||||||:!:||| !!|||..!  !  !|||! !!
           uValGluLeuLeuIleGlyArgTyrProValProAlaProSerGlnAlaGluTyrA
 3034219 : TGTCGAACTTCTTATTGGACGATATCCAGTGCCTGCACCGTCTCAAGCAGAATATG : 3034163

     273 : luLeuIlePheGlyCysSerVal<->GluGlyAspProAlaSerSerGluLeuAla :     289
               !!!:|||..!   :!!! !   |||  !!!:! !|||..!|||     !! !
           laThrMetPheAsnValAlaGluAsnGluIleGluLeuAlaAspSerLeuGluGlu
 3034162 : CTACAATGTTTAATGTGGCTGAAAATGAGATTGAACTGGCTGATAGTCTTGAAGAG : 3034109

     290 : ProArgProArgProProGlyArgProIleSer :     300
           |||!.!  !!.!||||||.!!!.!|||   |||
           ProAsnTyrHisProProSerAsnProAlaSer
 3034108 : CCCAACTATCATCCGCCGAGCAATCCGGCGTCA : 3034074

vulgar: frog 32 300 . CM041006.1 3039620 3034073 - 732 M 52 156 5 0 2 I 0 51 3 0 2 M 73 219 5 0 2 I 0 3833 3 0 2 M 32 96 S 0 1 5 0 2 I 0 46 3 0 2 S 1 2 M 14 42 F 0 2 G 1 0 M 26 78 5 0 2 I 0 795 3 0 2 M 48 144 G 0 3 M 21 63
-- completed exonerate analysis

In addition, exonerate lets you control the output you want. Check out the --ryo (roll your own output) option in man exonerate.

genewise

GeneWise compares a protein sequence to a genomic DNA sequence, allowing for introns and frameshifting errors.

Unlike exonerate, genewise only does protein to DNA alignment, but like exonerate it models introns and splicing. 20 years ago, when I did this every day, my impression was that genewise was slower but slightly more sensitive than the default exonerate options (exonerate can be tweaked to be made more sensitive). So I would usually try exonerate first and fall back to genewise if it failed to find a good match.


---

