Как объединить (объединить) AAStringSets по имени?

Question

В литературе по биоинформатике / микробной экологии довольно распространенной практикой является объединение множественных выравниваний последовательностей нескольких генов до построения филогенетических деревьев.В терминологии R может быть яснее сказать «объединить» эти последовательности организмом, из которого они произошли, но я уверен, что примеры лучше.

Скажем, это два множественных выравнивания последовательности.

library(Biostrings)

set1<-AAStringSet(c("IVR", "RDG", "LKS"))
names(set1)<-paste("org", 1:3, sep="_")

set2<-AAStringSet(c("VRT", "RKG", "AST"))
names(set2)<-paste("org", 2:4, sep="_")

set1

A AAStringSet instance of length 3
    width seq    names               
[1]     3 IVR    org_1
[2]     3 RDG    org_2
[3]     3 LKS    org_3

set2

A AAStringSet instance of length 3
    width seq    names               
[1]     3 VRT    org_2
[2]     3 RKG    org_3
[3]     3 AST    org_4

Правильная конкатенация этих последовательностей будет

A AAStringSet instance of length 4
    width seq    names               
[1]     6 IVR--- org_1
[2]     6 RDGVRT org_2
[3]     6 LKSRKG org_3
[4]     6 ---AST org_4

"-"отмечает «разрыв» (недостаток аминокислоты) в этом положении или в этом случае отсутствие гена для конкатенации.

Я думал, что будет функция для этого в BioStrings, MSA, DECIPHER или других связанных пакетах, но я не смог найти ее.

Я обнаружил следующие вопросы и ответы, каждый из которых не обеспечивает требуемый выход, как описано.

1: https://support.bioconductor.org/p/38955/

выход

  A AAStringSet instance of length 6
    width seq names               
[1]     3 IVR org_1
[2]     3 RDG org_2
[3]     3 LKS org_3
[4]     3 VRT org_2
[5]     3 RKG org_3
[6]     3 AST org_4

Можно лучше описать как «добавление» последовательностей (объединяет два набора по вертикали).

2: https://support.bioconductor.org/p/39878/

output

  A AAStringSet instance of length 2

        width seq
    [1]     9 IVRRDGLKS
    [2]     9 VRTRKGAST

Объединяет последовательности в каждом наборе, полную химеру каждого набора (определенно не желательно).

3: Как объединить две последовательности DNAStringSet на выборку в выводе R?

  A AAStringSet instance of length 3
    width seq
[1]     6 IVRVRT
[2]     6 RDGRKG
[3]     6 LKSAST

Создает химеры последовательностей в том порядке, в котором они находятся. Еще хуже с различным числом последовательностей (более короткое множество циклов и объединений ...)

4: https://www.biostars.org/p/115192/

Вывод

  A AAStringSet instance of length 2
    width seq
[1]     3 IVR
[2]     3 VRT

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

Я бы обычно думал, что такого рода процессы будутбыть сделано с некоторой комбинацией bash и Python, но я использую DECIPHER множественный выравниватель последовательности яn R, поэтому имеет смысл выполнить остальную часть обработки в R.В процессе написания этого вопроса я нашел ответ, который я опубликую, но я ожидаю, что кто-то укажет мне на пропущенное мной руководство, описывающее функцию, которая делает это.Спасибо!

Answer 1

Так что я несколько фанатичный пользователь data.table в R, среди многих вещей здорово объединять наборы данных по именам.Я обнаружил, что Biostrings::AAStringSet s могут быть преобразованы в матрицы с использованием as.matrix, и они могут быть преобразованы в data.table и объединены.

set1.dt<-data.table(as.matrix(set1), keep.rownames = TRUE)
set2.dt<-data.table(as.matrix(set2), keep.rownames = TRUE)
set12.dt<-merge(set1.dt, set2.dt, by="rn", all=TRUE)
    set12.dt
      rn V1.x V2.x V3.x V1.y V2.y V3.y
1: org_1    I    V    R <NA> <NA> <NA>
2: org_2    R    D    G    V    R    T
3: org_3    L    K    S    R    K    G
4: org_4 <NA> <NA> <NA>    A    S    T

Это правильное слияние, но требуется больше работы для получения окончательного вариантарезультат.

Необходимо заменить «NA» на «-».Мне всегда нужно искать этот вопрос, чтобы запомнить лучший способ сделать это с data.table.

Самый быстрый способ замены NA в больших данных. Таблица

#slightly modified from original, added arg "x"
f_dowle = function(dt, x) {     # see EDIT later for more elegant solution
      na.replace = function(v,value=x) { v[is.na(v)] = value; v }
      for (i in names(dt))
        eval(parse(text=paste("dt[,",i,":=na.replace(",i,")]")))
    }

f_dowle(set12.dt, "-")

Объединение последовательностей (не включая имена с !"rn")

set12<-apply(set12.dt[ ,!"rn"], 1, paste, collapse="")

Преобразование обратно в AAStringSet и добавление обратных имен

set12<-AAStringSet(set12)
names(set12)<-set12.dt$rn

Желаемый вывод

set12
 A AAStringSet instance of length 4
    width seq names               
[1]     6 IVR--- org_1
[2]     6 RDGVRT org_2
[3]     6 LKSRKG org_3
[4]     6 ---AST org_4

Это работает, но кажется довольно громоздким, особенно при преобразовании между различными форматами данных.Очевидно, что это может обернуть его в функцию для более легкого использования, но опять же кажется, что это уже должна быть функция в некотором пакете Bioconductor ...