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

CC BY-SA 4.0
2020-07-23T08:03:22.817Z
Perfect @RyanSYKwan its alot of code
M__
CC BY-SA 4.0
2020-07-23T12:31:05.270Z
thanks a lot this will be really really helpful for my future work
kcm
CC BY-SA 4.0
2020-07-23T13:21:02.640Z
Glad it helps! I learn a lot from you too.
Ryan SY Kwan
CC BY-SA 4.0
2020-07-23T06:05:01.893Z
2020-07-23T13:16:06.813Z
2020-07-23T13:16:06.813Z
Ryan SY Kwan
<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>

<r><phylogenetics>
Comparing aligned amino acids to codon
# 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:

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$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------
#      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------
#      5 7890  34
#         ^^^

#K9N5Q8 gene codon
#gtg ttt cta ctg atg ttc ttg tta aca
#                ^^^ ^^^ ^^^
#P0DTC2 gene codon
#atg ttt gtt ttt ctt
#^^^ ^^^ ^^^

> **Comments:**
> **username_2**: Perfect @RyanSYKwan 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.  

---

