<p>@StupidWolf answered me on stackoverflow and he's answer solved my question. I will share it here as well.</p>

<p>@StupidWolf:</p>

<p>We read in the alignment:</p>

<pre><code>library(Biostrings)
library(ggmsa)
protein_sequences &lt;- system.file("extdata", "sample.fasta", package = "ggmsa")
aln = readAAMultipleAlignment(protein_sequences)
ggmsa(protein_sequences, start = 265, end = 300) 

</code></pre>

<p><a href="https://i.sstatic.net/tU7Nd.png" rel="noreferrer"><img src="https://i.sstatic.net/tU7Nd.png" alt="enter image description here"></a>
Set the reference as the 1st sequence, some Rattus, you can also use the consensus with consensusString() :</p>

<pre><code>aln = unmasked(aln)
names(aln)[1]
[1] "PH4H_Rattus_norvegicus"

ref = aln[1]
</code></pre>

<p>Here we iterate through the sequence and make the binary for where the sequences are the same as the reference:</p>

<pre><code>bm = sapply(1:length(aln),function(i){
as.numeric(as.matrix(aln[i])==as.matrix(ref))
})

bm = t(bm)
rownames(bm) = names(aln)
</code></pre>

<p>The plot you see above has sequences reversed, so to visualize the same thing we reverse it, and also subset on 265 - 300:</p>

<pre><code>library(pheatmap)
pheatmap(bm[nrow(bm):1,265:300],cluster_rows=FALSE,cluster_cols=FALSE)

</code></pre>

<p><a href="https://i.sstatic.net/BA0vk.png" rel="noreferrer"><img src="https://i.sstatic.net/BA0vk.png" alt="enter image description here"></a>
The last row, is Rattus, the reference, and anything similar to that is read, as you can see in the alignment above, last 4 sequences are identical.</p>

CC BY-SA 4.0
2020-05-28T22:33:19.393Z
2020-05-28T22:33:19.393Z
rodrigarco
<p>I have many amino acids sequences in a fasta format that I did multiple sequence alignment. I was trying to plot something like binary a code as heatmap. If it had a change in the AA compared to the reference sequence it would be red, if it did not change would be yellow, for example.</p>

<p>I came across to <code>msaplot</code> from <code>ggtree</code> package. I also checked <code>ggmsa</code> package for that. But so far, I did not get what I wanted. So basically I wanted to:</p>

<ol>
<li><p>change the multiple sequence alignment to a binary matrix (if the amino differs from the reference sequence, plot x, if not y)</p></li>
<li><p>plot as a heatmap with two colors based if the amino acid changed or not</p></li>
</ol>

<p>here's some data example of what I doing, just instal <code>ggmsa</code> and you can try this:</p>

<pre><code>protein_sequences &lt;- system.file("extdata", "sample.fasta", package = "ggmsa")
ggmsa(protein_sequences, start = 265, end = 300, color = "Chemistry_AA",font = NULL)
</code></pre>

CC BY-SA 4.0
2020-05-28T20:44:33.277Z
Could you post your code?
M__
CC BY-SA 4.0
2020-05-28T20:55:20.037Z
I added a very simple example code present on ```ggmsa``` package
rodrigarco
<r><phylogenetics><heatmap>
How can I plot multiple sequence alignment in R as a heatmap
# How can I plot multiple sequence alignment in R as a heatmap

(Asked by: username_0 on 2020-05-28)

I have many amino acids sequences in a fasta format that I did multiple sequence alignment. I was trying to plot something like binary a code as heatmap. If it had a change in the AA compared to the reference sequence it would be red, if it did not change would be yellow, for example.
I came across to msaplot from ggtree package. I also checked ggmsa package for that. But so far, I did not get what I wanted. So basically I wanted to:

change the multiple sequence alignment to a binary matrix (if the amino differs from the reference sequence, plot x, if not y)
plot as a heatmap with two colors based if the amino acid changed or not

here's some data example of what I doing, just instal ggmsa and you can try this:
protein_sequences <- system.file("extdata", "sample.fasta", package = "ggmsa")
ggmsa(protein_sequences, start = 265, end = 300, color = "Chemistry_AA",font = NULL)

**Tags:** r, phylogenetics, heatmap
> **Comments:**
> **username_1**: Could you post your code?  
> **username_0**: I added a very simple example code present on ```ggmsa``` package  

---

## 1 answers

### Answer by username_0 on 2020-05-28. Score: 5 | Accepted answer

@StupidWolf answered me on stackoverflow and he's answer solved my question. I will share it here as well.
@StupidWolf:
We read in the alignment:
library(Biostrings)
library(ggmsa)
protein_sequences <- system.file("extdata", "sample.fasta", package = "ggmsa")
aln = readAAMultipleAlignment(protein_sequences)
ggmsa(protein_sequences, start = 265, end = 300) 



Set the reference as the 1st sequence, some Rattus, you can also use the consensus with consensusString() :
aln = unmasked(aln)
names(aln)[1]
[1] "PH4H_Rattus_norvegicus"

ref = aln[1]

Here we iterate through the sequence and make the binary for where the sequences are the same as the reference:
bm = sapply(1:length(aln),function(i){
as.numeric(as.matrix(aln[i])==as.matrix(ref))
})

bm = t(bm)
rownames(bm) = names(aln)

The plot you see above has sequences reversed, so to visualize the same thing we reverse it, and also subset on 265 - 300:
library(pheatmap)
pheatmap(bm[nrow(bm):1,265:300],cluster_rows=FALSE,cluster_cols=FALSE)



The last row, is Rattus, the reference, and anything similar to that is read, as you can see in the alignment above, last 4 sequences are identical.


---

