Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-23T08:03:22.817Z
Answers.Text: Perfect @RyanSYKwan its alot of code
Answers.UserDisplayName: M__
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-23T12:31:05.270Z
Answers.Text: thanks a lot this will be really really helpful for my future work
Answers.UserDisplayName: kcm
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-23T13:21:02.640Z
Answers.Text: Glad it helps! I learn a lot from you too.
Answers.UserDisplayName: Ryan SY Kwan
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-23T06:05:01.893Z
Answers.LastActivityDate: 2020-07-23T13:16:06.813Z
Answers.LastEditDate: 2020-07-23T13:16:06.813Z
Answers.UserDisplayName: Ryan SY Kwan
Body: <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>
  
Tags: <r><phylogenetics>
Title: Comparing aligned amino acids to codon
ThreadText: # Comparing aligned amino acids to codon
  
  (Asked by: username_0 on 2020-07-06)
  
  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.
  
  1. To start with I have used the seqinr library to generate nucleotide alignment from aligned protein sequence.
  
     reverse.align {seqinr}
  2. Then the next step that i want to compare each of my amino acid sequence to my reference sequence which in my case is `P0DTC2`to rest of the sequence for example I want to compare `P0DTC2` to `K9N5Q8`to the above mentioned CTD region which starts from `NITNLC` and ends with `HAPATV`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.
  
  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](https://drive.google.com/file/d/1J4-eE99I6MQgNw6m6sFE60Qrt_4Bg7KJ/view?usp=sharing) [nucleotide file](https://drive.google.com/file/d/14cyaxEq84fYlHKbLjZmIhvGf60gakBKB/view?usp=sharing) [Aligned amino acid file](https://drive.google.com/file/d/1FX-VAQixASrmiYZc-mIdHHOnwGtGZJ_j/view?usp=sharing) [Reverse alignment nucleotide file using seqinr](https://drive.google.com/file/d/1dyQHqYEj5TXt1CBGwieRG48odghBIYTq/view?usp=sharing)
  
  **Tags:** r, phylogenetics
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2020-07-23. Score: 2 | Accepted answer
  
  To set the ball rolling, I tried to implement it manually with your files. Overall, there are three steps:
  
  1. Locate `NITNLC` (or `HAPATV`) in protein `P0DTC2`.
  2. Locate mismatched amino acids between proteins `P0DTC2` and `K9N5Q8` within the range from step 1.
  3. 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
  #^^^ ^^^ ^^^ 
  
  ```
  
  > **Comments:**
  > **username_2**: Perfect @username_1 its alot of code  
  > **username_0**: thanks a lot this will be really really helpful for my future work  
  > **username_1**: Glad it helps! I learn a lot from you too.  
  
  ---
  
  