text: user285: I have a set of 4 amino acids which i have aligned,but i want to compare then with respective  nucleotide sequence in terms of change codon level.
  I have the alignment file where i want to take only for the C-terminal domain region region only the CTD starts with amino acids “NITNLC” till it ends with “HAPATV” and rest of the sequence i dont want to take.
  In short i want to take the part sequence which starts from "NITNLC" to "HAPATV" then i would like to compare what are the changes in the nucleotide sequence at codon level.
  
  To start with I have used the seqinr library to generate nucleotide alignment from aligned protein sequence.
  reverse.align {seqinr}
  
  Then the next step that i want to compare each of my amino acid sequence to my reference sequence which in my case is P0DTC2to rest of the sequence for example I want to compare P0DTC2 to K9N5Q8to the above mentioned CTD region which starts from NITNLC and ends with HAPATVand find what are the codon level changes if there is any I would like to report both in amino acid level and codon level.
  
  
  First part I guess my approach is correct. Second part I'm not sure how to proceed i guess it's much more than simple parsing I guess!!.
  Any help or suggestion would be really appreciated ,and if there exist a R based solution that would be really welcome.
  Files
  Amino acid file nucleotide file Aligned amino acid file Reverse alignment nucleotide file using seqinr
  
  user9172: To set the ball rolling, I tried to implement it manually with your files. Overall, there are three steps:
  
  Locate NITNLC (or HAPATV) in protein P0DTC2.
  Locate mismatched amino acids between proteins P0DTC2 and K9N5Q8 within the range from step 1.
  Print the amino acid and the DNA codon of both proteins from step 2.
  
  It works, but only for the first 60 amino acids. I wonder if it's due to the amino sequences in AMINOOO_seq_removed.fasta repeat itself every 60 acids. But why?
  #For example, the first three lines of protein P0DTC2
  >P0DTC2
  MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDLFLPFFS MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDLFLPFFS
  NVTWFHAIHVSGTNGTKRFDNPVLPFNDGVYFASTEKSNIIRGWIFGTTLDSKTQSLLIV NVTWFHAIHVSGTNGTKRFDNPVLPFNDGVYFASTEKSNIIRGWIFGTTLDSKTQSLLIV
  NNATNVVIKVCEFQFCNDPFLGVYYHKNNKSWMESEFRVYSSANNCTFEYVSQPFLMDLE NNATNVVIKVCEFQFCNDPFLGVYYHKNNKSWMESEFRVYSSANNCTFEYVSQPFLMDLE
  ...
  
  Step 0: Read the files.
  library('seqinr')
  
  #align file containing protein sequences
  count_added <- read.alignment('count_added_.clustal_num', format='clustal')
  names(count_added$seq) <- count_added$nam
  
  #DNA sequences
  rev3.aln <- read.alignment('rev3.aln', format='fasta')
  names(rev3.aln$seq) <- rev3.aln$nam
  
  Step 1: Locate NITNLC (or HAPATV) in protein P0DTC2. Again, there are two repeats of NITNLC 60 acids apart (at 807 and 870).
  locate <- function(seq, find)
    {address <- gregexpr(paste(strsplit(find, '')[[1]], collapse='[^a-z]*'), seq)
    #substr(seq, address[[1]][1], address[[1]][1]+attr(address[[1]], 'match.length')[1]-1)
  
    return(list(start=as.numeric(address[[1]]), 
               end=as.numeric(address[[1]] + attr(address[[1]], 'match.length') - 1)))
    }
  
  locate(seq=count_added$seq[['P0DTC2']], find='nitnlc')
  #start
  #807 870
  #end
  #812 875
  locate(seq=count_added$seq[['P0DTC2']], find='hapatv')
  #start
  #1212 1279
  #end
  #1217 1285
  
  Step 2: Locate mismatched amino acids between proteins P0DTC2 and K9N5Q8 within the range from step 1. Instead of using the range 812 - 1279, I chose range 1 - 20 for demonstration purposes.
  compare <- function(seq1, seq2, after=0, before=100000)
    {seq1_ = strsplit(seq1, '')[[1]]
    seq2_ = strsplit(seq2, '')[[1]]
    ind = which(seq1_ != seq2_ & grepl('[a-z]',seq1_) & grepl('[a-z]',seq2_))
    ind = ind[ind>after & ind<before]
    #seq1_[ind[1]]
    #seq2_[ind[1]]
  
    return(ind)
    }
  compare(seq1=count_added$seq[['P0DTC2']], seq2=count_added$seq[['K9N5Q8']], after=1, before=20)
  # [1] 5  7  8  9 10 13 14
  
  #protein comparison
  #K9N5Q8      MIHSVFLLMFLLTPTESYVD
  #P0DTC2      ----MFVFLVLLPL------
  #<mismatch>      5 7890  34
  
  Step 3: Print the amino acid and the DNA codon of both proteins from step 2. Note that ind is based on the align file.
  print_amino_codon <- function(ind, seq, seq_gene)
    {locate_amino <- gregexpr('[a-z]', seq)[[1]]
    if (!ind %in% locate_amino) return(NA)
    ind2 = match(ind, locate_amino)
  
    return(c(amino=substr(seq, ind, ind), codon=substr(seq_gene, ind2*3-2, ind2*3)))
    }
   
  codon(ind=5, seq=count_added$seq[['K9N5Q8']], seq_gene=rev3.aln$seq[['K9N5Q8']])
  #amino codon 
  #  "v" "atg" 
  codon(ind=6, vseq=count_added$seq[['K9N5Q8']], seq_gene=rev3.aln$seq[['K9N5Q8']])
  #amino codon 
  #  "f" "ttc" 
  codon(ind=7, seq=count_added$seq[['K9N5Q8']], seq_gene=rev3.aln$seq[['K9N5Q8']])
  #amino codon 
  #  "l" "ttg"
  
  codon(ind=5, seq=count_added$seq[['P0DTC2']], seq_gene=rev3.aln$seq[['P0DTC2']])
  #amino codon 
  #  "m" "atg" 
  codon(ind=6, seq=count_added$seq[['P0DTC2']], seq_gene=rev3.aln$seq[['P0DTC2']])
  #amino codon 
  #  "f" "ttt" 
  codon(ind=7, seq=count_added$seq[['P0DTC2']], seq_gene=rev3.aln$seq[['P0DTC2']])
  #amino codon 
  #  "v" "gtt" 
  
  #protein comparison
  #K9N5Q8      MIHSVFLLMFLLTPTESYVD
  #P0DTC2      ----MFVFLVLLPL------
  #<mismatch>      5 7890  34
  #<print>         ^^^
  
  #K9N5Q8 gene codon
  #gtg ttt cta ctg atg ttc ttg tta aca
  #                ^^^ ^^^ ^^^
  #P0DTC2 gene codon
  #atg ttt gtt ttt ctt
  #^^^ ^^^ ^^^ 
  
  
original_text: 285: <p>I have a set of 4 amino acids which i have aligned,but i want to compare then with respective  nucleotide sequence in terms of change codon level.</p>
  <p>I have the alignment file where i want to take only for the C-terminal domain region region only the CTD starts with amino acids <code>“NITNLC”</code> till it ends with <code>“HAPATV”</code> and rest of the sequence i dont want to take.</p>
  <p>In short i want to take the part sequence which starts from &quot;NITNLC&quot; to &quot;HAPATV&quot; then i would like to compare what are the changes in the nucleotide sequence at codon level.</p>
  <ol>
  <li><p>To start with I have used the seqinr library to generate nucleotide alignment from aligned protein sequence.</p>
  <p>reverse.align {seqinr}</p>
  </li>
  <li><p>Then the next step that i want to compare each of my amino acid sequence to my reference sequence which in my case is <code>P0DTC2</code>to rest of the sequence for example I want to compare <code>P0DTC2</code> to <code>K9N5Q8</code>to the above mentioned CTD region which starts from <code>NITNLC</code> and ends with <code>HAPATV</code>and find what are the codon level changes if there is any I would like to report both in amino acid level and codon level.</p>
  </li>
  </ol>
  <p>First part I guess my approach is correct. Second part I'm not sure how to proceed i guess it's much more than simple parsing I guess!!.</p>
  <p>Any help or suggestion would be really appreciated ,and if there exist a R based solution that would be really welcome.</p>
  <p>Files
  <a href="https://drive.google.com/file/d/1J4-eE99I6MQgNw6m6sFE60Qrt_4Bg7KJ/view?usp=sharing" rel="nofollow noreferrer">Amino acid file</a> <a href="https://drive.google.com/file/d/14cyaxEq84fYlHKbLjZmIhvGf60gakBKB/view?usp=sharing" rel="nofollow noreferrer">nucleotide file</a> <a href="https://drive.google.com/file/d/1FX-VAQixASrmiYZc-mIdHHOnwGtGZJ_j/view?usp=sharing" rel="nofollow noreferrer">Aligned amino acid file</a> <a href="https://drive.google.com/file/d/1dyQHqYEj5TXt1CBGwieRG48odghBIYTq/view?usp=sharing" rel="nofollow noreferrer">Reverse alignment nucleotide file using seqinr</a></p>
   9172: <p>To set the ball rolling, I tried to implement it manually with your files. Overall, there are three steps:</p>
  <ol>
  <li>Locate <code>NITNLC</code> (or <code>HAPATV</code>) in protein <code>P0DTC2</code>.</li>
  <li>Locate mismatched amino acids between proteins <code>P0DTC2</code> and <code>K9N5Q8</code> within the range from step 1.</li>
  <li>Print the amino acid and the DNA codon of both proteins from step 2.</li>
  </ol>
  <p>It works, but only for the first 60 amino acids. I wonder if it's due to the amino sequences in <code>AMINOOO_seq_removed.fasta</code> repeat itself every 60 acids. But why?</p>
  <pre><code>#For example, the first three lines of protein P0DTC2
  &gt;P0DTC2
  MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDLFLPFFS MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDLFLPFFS
  NVTWFHAIHVSGTNGTKRFDNPVLPFNDGVYFASTEKSNIIRGWIFGTTLDSKTQSLLIV NVTWFHAIHVSGTNGTKRFDNPVLPFNDGVYFASTEKSNIIRGWIFGTTLDSKTQSLLIV
  NNATNVVIKVCEFQFCNDPFLGVYYHKNNKSWMESEFRVYSSANNCTFEYVSQPFLMDLE NNATNVVIKVCEFQFCNDPFLGVYYHKNNKSWMESEFRVYSSANNCTFEYVSQPFLMDLE
  ...
  </code></pre>
  <p>Step 0: Read the files.</p>
  <pre><code>library('seqinr')
  
  #align file containing protein sequences
  count_added &lt;- read.alignment('count_added_.clustal_num', format='clustal')
  names(count_added<span class="math-container">$seq) &lt;- count_added$</span>nam
  
  #DNA sequences
  rev3.aln &lt;- read.alignment('rev3.aln', format='fasta')
  names(rev3.aln<span class="math-container">$seq) &lt;- rev3.aln$</span>nam
  </code></pre>
  <p>Step 1: Locate <code>NITNLC</code> (or <code>HAPATV</code>) in protein <code>P0DTC2</code>. Again, there are two repeats of <code>NITNLC</code> 60 acids apart (at 807 and 870).</p>
  <pre><code>locate &lt;- function(seq, find)
    {address &lt;- gregexpr(paste(strsplit(find, '')[[1]], collapse='[^a-z]*'), seq)
    #substr(seq, address[[1]][1], address[[1]][1]+attr(address[[1]], 'match.length')[1]-1)
  
    return(list(start=as.numeric(address[[1]]), 
               end=as.numeric(address[[1]] + attr(address[[1]], 'match.length') - 1)))
    }
  
  locate(seq=count_added<span class="math-container">$seq[['P0DTC2']], find='nitnlc')
  #start
  #807 870
  #end
  #812 875
  locate(seq=count_added$</span>seq[['P0DTC2']], find='hapatv')
  #start
  #1212 1279
  #end
  #1217 1285
  </code></pre>
  <p>Step 2: Locate mismatched amino acids between proteins <code>P0DTC2</code> and <code>K9N5Q8</code> within the range from step 1. Instead of using the range 812 - 1279, I chose range 1 - 20 for demonstration purposes.</p>
  <pre><code>compare &lt;- function(seq1, seq2, after=0, before=100000)
    {seq1_ = strsplit(seq1, '')[[1]]
    seq2_ = strsplit(seq2, '')[[1]]
    ind = which(seq1_ != seq2_ &amp; grepl('[a-z]',seq1_) &amp; grepl('[a-z]',seq2_))
    ind = ind[ind&gt;after &amp; ind&lt;before]
    #seq1_[ind[1]]
    #seq2_[ind[1]]
  
    return(ind)
    }
  compare(seq1=count_added<span class="math-container">$seq[['P0DTC2']], seq2=count_added$</span>seq[['K9N5Q8']], after=1, before=20)
  # [1] 5  7  8  9 10 13 14
  
  #protein comparison
  #K9N5Q8      MIHSVFLLMFLLTPTESYVD
  #P0DTC2      ----MFVFLVLLPL------
  #&lt;mismatch&gt;      5 7890  34
  </code></pre>
  <p>Step 3: Print the amino acid and the DNA codon of both proteins from step 2. Note that <code>ind</code> is based on the align file.</p>
  <pre><code>print_amino_codon &lt;- function(ind, seq, seq_gene)
    {locate_amino &lt;- gregexpr('[a-z]', seq)[[1]]
    if (!ind %in% locate_amino) return(NA)
    ind2 = match(ind, locate_amino)
  
    return(c(amino=substr(seq, ind, ind), codon=substr(seq_gene, ind2*3-2, ind2*3)))
    }
   
  codon(ind=5, seq=count_added<span class="math-container">$seq[['K9N5Q8']], seq_gene=rev3.aln$</span>seq[['K9N5Q8']])
  #amino codon 
  #  &quot;v&quot; &quot;atg&quot; 
  codon(ind=6, vseq=count_added<span class="math-container">$seq[['K9N5Q8']], seq_gene=rev3.aln$</span>seq[['K9N5Q8']])
  #amino codon 
  #  &quot;f&quot; &quot;ttc&quot; 
  codon(ind=7, seq=count_added<span class="math-container">$seq[['K9N5Q8']], seq_gene=rev3.aln$</span>seq[['K9N5Q8']])
  #amino codon 
  #  &quot;l&quot; &quot;ttg&quot;
  
  codon(ind=5, seq=count_added<span class="math-container">$seq[['P0DTC2']], seq_gene=rev3.aln$</span>seq[['P0DTC2']])
  #amino codon 
  #  &quot;m&quot; &quot;atg&quot; 
  codon(ind=6, seq=count_added<span class="math-container">$seq[['P0DTC2']], seq_gene=rev3.aln$</span>seq[['P0DTC2']])
  #amino codon 
  #  &quot;f&quot; &quot;ttt&quot; 
  codon(ind=7, seq=count_added<span class="math-container">$seq[['P0DTC2']], seq_gene=rev3.aln$</span>seq[['P0DTC2']])
  #amino codon 
  #  &quot;v&quot; &quot;gtt&quot; 
  
  #protein comparison
  #K9N5Q8      MIHSVFLLMFLLTPTESYVD
  #P0DTC2      ----MFVFLVLLPL------
  #&lt;mismatch&gt;      5 7890  34
  #&lt;print&gt;         ^^^
  
  #K9N5Q8 gene codon
  #gtg ttt cta ctg atg ttc ttg tta aca
  #                ^^^ ^^^ ^^^
  #P0DTC2 gene codon
  #atg ttt gtt ttt ctt
  #^^^ ^^^ ^^^ 
  </code></pre>
  