R中的DNA条件频率

DNA conditional frequency in R

我试图找出 R 中 2 个不同 DNA 序列中是否存在任何条件依赖性

这是我的代码,但是我遇到了错误;

Error in `[.data.frame`(data, i) : undefined columns selected

我不确定问题出在哪里,如果我将 data[i-1]==bases[b2] 括起来,我只会得到多个 unexpected},这是我唯一能想到的事情。

for (b1 in 1:length(bases))
{
   for (b2 in 1:length(bases))
   {
       count = 1
       for (i in 2:length(mydata1))
       {
           if ((mydata1[i]==bases[b1]) & mydata1[i-1]==bases[b2])
           {
               count = count+1
           }
       }
       b3 = c(bases[b1], bases[b2], count)
       print(b3)
   }
}

_我基本上期待某些 DNA 碱基的列表,例如,我认为 DNA 序列好像是以先前的碱基为条件的;

[1] "A" "C" "002"
[1] "A" "C" "005"
[1] "A" "C" "009"

等等,通过清楚地显示 A 的条件,可以告诉我某个碱基是否对后续碱基的身份有任何影响C.

Ok 所以本质上 mydata1(还有 mydata2)是 DNA 序列,也就是说 "A", "G", "C" and "T" 的列表,每个都是10,000 个碱基长。

如图所示;

  V1
1  T
2  C
3  G
4  G
5  T
6  G
7  G 
8  G
9  C
10 A

我的任务是尝试确定该序列是否具有相互依赖的碱基,因此如果 [1] T 影响 [2] C 的存在,等等。一个一个序列是相关的,另一个不是。

如果我没理解错的话,你想计算DNA序列中每对核苷酸i,i+1出现的次数。您可以使用 R 函数 table 实现此目的;下面提供了一个示例。

# input sequence
seq <- "ACGTACTGCACAAACTAC"

# length of input sequence
length_seq <- nchar(seq, type="chars")

# first substring: from 1 to second-last 
seq1 <- substr(seq, 1, (length_seq - 1))

# second substring: from 2 to last
seq2 <- substring(seq, 2, length_seq)

# split strings
seq1_split <- strsplit(seq1, "")[[1]]
seq2_split <- strsplit(seq2, "")[[1]]

# initialize vectors
first_nt <- vector(mode="character", length = (length_seq - 1))
second_nt <- vector(mode="character", length = (length_seq -1))

# fill vectors
count = 0
for (b in seq1_split)
{
    count = count + 1
    first_nt[count] <- b
}

count = 0
for (b in seq2_split)
{
    count = count + 1
    second_nt[count] <- b
}

# create matrix with character i and i+1 in each row
mat <- matrix(c(first_nt, second_nt), nrow=(length_seq - 1))

# collapse matrix
to_table <- apply(mat, 1, paste, collapse="")

# table
my_table <- table(to_table)

print(my_table)