Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-10-13T12:03:12.847Z
Answers.LastActivityDate: 2023-10-13T12:03:12.847Z
Answers.UserDisplayName: terdon
Body: <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>
  
Tags: <phylogenetics>
Title: Find protein in DNA sequences with arbitrary encoding and possible frameshifts
ThreadText: # 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](https://www.ncbi.nlm.nih.gov/books/NBK62051/) 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](https://www.ebi.ac.uk/about/vertebrate-genomics/software/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](https://www.ebi.ac.uk/Tools/psa/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.
  
  
  ---
  
  