6

次のRコード行があります。

croppedDNA <- completeDNA[,apply(completeDNA,2,function(x) any(c(FALSE,x[-length(x)]!=x[-1])))]

これが行うことは、DNA 配列のマトリックス (1 行 = 1 配列) 内で普遍的ではない (有益な) サイト (列) を特定し、それらをマトリックスからサブセット化して、新しい「トリミングされたマトリックス」を作成することです。値が同じ列。大きなデータセットの場合、これには約 6 秒かかります。C++ (まだ C++ の初心者) で高速に実行できるかどうかはわかりませんが、試してみるとよいでしょう。私の考えは、Rcpp を使用し、CharacterMatrix の列をループし、同じかどうかを CharacterVector チェックとして列 (サイト) を引き出すことです。それらが同じ場合、その列番号/インデックスを記録し、すべての列について続行します。最後に、それらの列のみを含む新しい CharacterMatrix を作成します。行名と列名を行列の「R バージョン」にあるように維持することが重要です。

私は約2分間書いてきましたが、これまでのところ（まだ終わっていません）：

#include <Rcpp.h>
#include <vector>
using namespace Rcpp;
// [[Rcpp::export]]
CharacterMatrix reduce_sequences(CharacterMatrix completeDNA)
{
  std::vector<bool> informativeSites; 
  for(int i = 0; i < completeDNA.ncol(); i++)
  {
    CharacterVector bpsite = completeDNA(,i);
    if(all(bpsite == bpsite[1])
    {
      informativeSites.push_back(i);
    }
  }
CharacterMatrix cutDNA = completeDNA(,informativeSites);
return cutDNA;
}

私はこれについて正しい道を進んでいますか？もっと簡単な方法はありますか。私の理解では、 std::vector が必要なのは、それらを簡単に成長させることができるためです(保持したい列の数が事前にわからないため)。インデックス付けでは、最後に infoativeSites ベクトルに +1 する必要がありますか (R は 1 から、C++ は 0 からインデックス付けするため)?

ありがとう、ベン W.

4

1 に答える 1

14

サンプルデータ：

set.seed(123)
z <- matrix(sample(c("a", "t", "c", "g", "N", "-"), 3*398508, TRUE), 3, 398508)

OPの解決策：

system.time(y1 <- z[,apply(z,2,function(x) any(c(FALSE,x[-length(x)]!=x[-1])))])
#    user  system elapsed 
#   4.929   0.043   4.976

ベース R を使用した高速バージョン:

system.time(y2 <- (z[, colSums(z[-1,] != z[-nrow(z), ]) > 0]))
#    user  system elapsed 
#   0.087   0.011   0.098

結果は同じです。

identical(y1, y2)
# [1] TRUE

C++ がそれを打ち負かす可能性は非常に高いですが、本当に必要なのでしょうか?

于 2013-05-15T02:52:13.943 に答える