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 に答える