如何在R中绘制共识序列的二元热图？

Question

我有许多氨基酸序列，如在fasta格式，我做了多个序列对齐。我试图绘制一些像二进制代码作为热图。如果它有一个变化，它将是红色的，如果它没有改变将是黄色的，例如。

我遇到了 msaplot 从 ggtree包。我还检查了 ggmsa 包的。但到目前为止，我并没有得到我想要的东西。所以基本上我想

1. 将多序列对齐改为二进制矩阵（如果氨基与参考序列不同，则绘制x，如果不是y）。

2. 作热图

多序列对齐是这样的，对于那些不知道的人来说。

我知道我应该提供某种类型的数据例子，但我不知道如何创建一个多序列对齐的例子。

如果您安装 ggmsa 你可以有一个例子的数据和绘图在R使用。

protein_sequences <- system.file("extdata", "sample.fasta", package = "ggmsa")
ggmsa(protein_sequences, start = 265, end = 300) 

Answer 1

我们读的对齐。

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

设置参考作为第一序列， 一些Rattus，你也可以使用共识与 consensusString() :

aln = unmasked(aln)
names(aln)[1]
[1] "PH4H_Rattus_norvegicus"

ref = aln[1]

这里我们对序列进行迭代，并对其中序列与参考相同的地方进行二进制。

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

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

你看上面的图有序列反转， 所以为了使同样的事情可视化，我们把它反转， 并在265 - 300上子集。

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

最后一行，是Rattus，参考， 和任何类似的是读取， 你可以看到在上面的对齐， 最后4个序列是相同的。