根据 A、C、G 或 T 的比例自动拆分不明确的碱基

Question

在 Biostrings 中，我加载了一个包含 427,351 个长度为 11 个核苷酸的 DNA 序列的 fasta 文件。

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

您可以看到，我的一个序列中第 11 位（第 N 行，第 11 列）有一个 "N" 核苷酸。

下一步是制作位置矩阵频率，但是，这只有在行的列总和为 "A"、"C"、"G" 和 [=37 时才有可能=] 是相等的。在上面的示例中，由于 N 基数，第 11 列的总和将比所有其他列少一。

编写 consensusMatrix 函数以便将所有非 A、C、G 和 T 碱基适当分类为 A、C、G、T 或它们的组合的最佳方法是什么？由于 N 代表 4 个碱基中的任何一个，因此对于 N 的每个实例，第 11 列的 A、C、G 和 T 值将增加 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 以外的任何序列。

Answer 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

根据 A、C、G 或 T 的比例自动拆分不明确的碱基

Automatically split ambiguous bases proportionally to A, C, G or T

r

bioinformatics