0

我有一组遗传 SNP 数据，如下所示：

Founder1 Founder2 Founder3 Founder4 Founder5 Founder6 Founder7 Founder8 Sample1 Sample2 Sample3 Sample...
A A A T T T T T A T A T
A A A T T T T T A T A T
A A A T T T T T A T A T
A A A T T T T T A T A T
A A A T T T T T A T A T
A A A T T T T T A T A T
A A A T T T T T A T A T
A A A T T T T T A T A T
A A A T T T T T A T A T
A A A T T T T T A T A T
A A A T T T T T A T A T
A A A T T T T T A T A T

矩阵的大小为 56 列 x 46482 行。我需要首先将矩阵按每 20 行进行一次分类，然后将前 8 列（创始人）中的每一列（创始人）与每列 9-56 进行比较，并将匹配字母/等位基因的总数除以总行数（20）。最终我需要 48 个 8 列 x 2342 行矩阵，它们本质上是相似矩阵。我试图通过以下方式分别提取每一对：

"length(cbind(odd[,9],odd[,1])[cbind(odd[,9],cbind(odd[,9],odd[,1])[,1])[,1]=="T" & cbind(odd[,9],odd[,1])[,2]=="T",])/nrow(cbind(odd[,9],odd[,1]))"

但这远非高效，而且我不知道将函数应用于每 20 行和多对的更快方法。

在上面给出的示例中，如果行全部相同，如跨 20 行所示，则 Sample1 矩阵的第一行将是：

1 1 1 0 0 0 0
4

1 回答 1

0

我想这就是你想要的？它有助于将问题分解成更小的部分，然后对这些部分重复应用函数。我的解决方案需要几分钟才能在我的笔记本电脑上运行，但我认为它应该让您或其他人开始。如果您正在寻找更快的速度，我建议您查看data.table包装。我敢肯定还有其他方法可以使下面的代码更快一点。

# Make a data set of random sequences
rows = 46482
cols = 56
binsize = 20
founder.cols = 1:8
sample.cols = setdiff(1:cols,founder.cols)
data = as.data.frame( matrix( sample( c("A","C","T","G"), 
                                      rows * cols, replace=TRUE ), 
                              ncol=cols ) )

# Split the data into bins
binlevels = gl(n=ceiling(rows/binsize),k=20,length=rows)
databins = split(data,binlevels)

# A function for making a similarity matrix
compare_cols = function(i,j,mat) mean(mat[,i] == mat[,j])
compare_group_cols = function(mat, group1.cols, group2.cols) {
  outer( X=group1.cols, Y=group2.cols, 
        Vectorize( function(X,Y) compare_cols(X,Y,mat) ) )
}

# Apply the function to each bin
mats = lapply( databins, compare_group_cols, sample.cols, founder.cols )

# And just to check. Random sequences should match 25% of the time. Right?
hist( vapply(mats,mean,1), n=30 ) # looks like this is the case
于 2013-11-12T03:30:53.183 回答