Автоматически разбивает неоднозначные основания пропорционально A, C, G или T
В Biostrings я загрузил файл fasta из 427 351 последовательности ДНК длиной 11 нуклеотидов.
my.seq<-readDNAStringSet("my.fasta", "fasta")
Затем я сгенерировал матрицу, которая подсчитывает общее количество определенного нуклеотида в каждой из 11 позиций:
my.pfm<-consensusMatrix(my.seq)
>my.pfm
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
A 113370 120216 109984 40729 150681 11 340936 41684 75946 150648 84290
C 98927 107171 99251 110222 76286 427265 25668 256664 191010 103889 139625
G 118545 93632 95588 74975 138899 9 95 91414 64966 66896 113694
T 96509 106332 122528 201425 61485 66 60652 37589 95429 105918 89741
M 0 0 0 0 0 0 0 0 0 0 0
R 0 0 0 0 0 0 0 0 0 0 0
W 0 0 0 0 0 0 0 0 0 0 0
S 0 0 0 0 0 0 0 0 0 0 0
Y 0 0 0 0 0 0 0 0 0 0 0
K 0 0 0 0 0 0 0 0 0 0 0
V 0 0 0 0 0 0 0 0 0 0 0
H 0 0 0 0 0 0 0 0 0 0 0
D 0 0 0 0 0 0 0 0 0 0 0
B 0 0 0 0 0 0 0 0 0 0 0
N 0 0 0 0 0 0 0 0 0 0 1
- 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0
. 0 0 0 0 0 0 0 0 0 0 0
Вы можете видеть, что у меня есть один "N" нуклеотид, присутствующий в одной из моих последовательностей в 11-й позиции (строка N, столбец 11).
Следующим шагом является создание матрицы положения частоты, однако это возможно только тогда, когда суммы столбцов строк "A", "C", "G" и "T" равны. В приведенном выше примере сумма столбца 11 будет на единицу меньше, чем у всех остальных из-за основания N.
Каков наилучший способ написать функцию consensusMatrix, чтобы все не-A,C,G и T основания были соответствующим образом классифицированы как A, C, G, T или их комбинация? Поскольку N представляет собой любое из 4 оснований, то для каждого экземпляра N к значениям A, C, G и T столбца 11 добавляется 0,25. Однако следует написать функцию для всех остальных не-A,C,G и T нуклеотидов, поэтому они соответственно назначены в правильной пропорции к A,C,G,T?
Например, Y= C или T, поэтому для каждого экземпляра Y 0,5 будет добавлено к C, а 0,5 будет добавлено к значению T этого столбца. Я вижу проблему, если у нас есть что-то вроде кода V, так как это может быть G, A или C, в этом случае 0.33333 будет добавлено к каждому экземпляру V для этого столбца.
Что я пробовал:
my.pfm<-consensusMatrix(my.seq,ambiguityMap=IUPAC_CODE_MAP)
Error in .local(x, as.prob, shift, width, ...) :
unused argument (ambiguityMap = c("A", "C", "G", "T", "AC", "AG", "AT", "CG", "CT", "GT", "ACG", "ACT", "AGT", "CGT", "ACGT"))
Как я понимаю, должен быть какой-то символьный вектор, который сообщает функции, что делать, когда подсчитывается что-либо, кроме A, C, G, T, но я не могу понять это.
Идея здесь в том, чтобы иметь какую-то определяющую функцию, которая будет работать независимо от того, сколько и каких кодов неоднозначности встретятся в будущем.
Примечание. Я не хочу удалять из набора данных целые последовательности, в которых есть что-то кроме A, C, G или T.
1 ответ
Примерно так, но из комментариев кажется, что вы задаете неправильный вопрос для такого типа данных.
#get sum of non ACGT and divide by 4
props <- colSums(my.pfm[ !rownames(my.pfm) %in% c("A","C","G","T"),]) / 4
#add it back to ACGT rows
t(
apply(
my.pfm[ rownames(my.pfm) %in% c("A","C","G","T"),], 1, function(i)
props + i))
#output
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
# A 113370 120216 109984 40729 150681 11 340936 41684 75946 150648 84290.25
# C 98927 107171 99251 110222 76286 427265 25668 256664 191010 103889 139625.25
# G 118545 93632 95588 74975 138899 9 95 91414 64966 66896 113694.25
# T 96509 106332 122528 201425 61485 66 60652 37589 95429 105918 89741.25