2
$\begingroup$

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 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?

$\endgroup$
0

1 Answer 1

2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.