Как я могу построить консенсусные последовательности в виде двоичной тепловой карты в R - PullRequest
Новая
       13

Как я могу построить консенсусные последовательности в виде двоичной тепловой карты в R

1 голос
/ 28 мая 2020

У меня есть много аминокислотных последовательностей, например, в формате fasta, в котором я сделал множественное выравнивание последовательностей. Я пытался построить что-то вроде двоичного кода в виде тепловой карты. Если бы в нем было изменение, оно было бы красным, если бы не изменилось, было бы, например, желтым.

Я наткнулся на msaplot из ggtree пакета. Я также проверил на это пакет ggmsa. Но пока не получил то, что хотел. Итак, в основном я хотел:

  1. изменить множественное выравнивание последовательностей на двоичную матрицу (если амино отличается от эталонной последовательности, отобразите x, если не y)

  2. график как тепловая карта

Множественное выравнивание последовательностей - это что-то вроде этого для тех, кто не знает. enter image description here

Я знаю, что должен предоставить какой-то пример данных, но я не уверен, как создать пример множественного выравнивания последовательностей

если вы установите ggmsa, вы можете получить пример данных и график в r используя: 

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

1 Ответ

1 голос
/ 29 мая 2020

Читаем в раскладе: 

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

enter image description here

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)

введите описание изображения здесь

Последняя строка - это Rattus, ссылка, и все, что похоже на это, читается, как вы можете видеть в приведенном выше выравнивании, последние 4 последовательности идентичны.

...