2016-03-23 19 views
1

私は別のテーブルに基づいてあるテーブルから値を取得する関数を実装しようとしています。実際のデータフレームは50,000を超える観測値を持つため、このネストされたforループの実装は効果的ではありません。私は過去数日間、それを見て、うまくいくものを見つけようとしましたが、できなかったものを探しました。私のデータは特別な順序ではありません(個人、セグメントなど)ので、物事が順不同であっても仕事ができる必要があります。ここでR - ネストされたループと遅いパフォーマンス

はで動作するように私のデータのおもちゃの例は以下のとおりです。

region_map <- data.frame(Start = c(721290, 1688193), End= c(1688192, 2926555)) 
individual <- c("Ind1","Ind2","Ind3","Ind4") 
segment <- data.frame(SampleID = c("Ind1","Ind1","Ind2","Ind2","Ind3","Ind3","Ind4","Ind4","Ind4"), 
         Start = c(721290, 1688194, 721290, 1688200, 721290, 2926600, 721290, 1688193, 690), 
         End = c(1688192, 2926555,1688190, 2900000, 2926555, 3000000, 1500000, 2005000, 500000), 
         State = c(1,2,2,5,4,2,2,6,5)) 

そして、ここでは、私が何をしようとしているの簡単な例です:すなわち

Generate.FullSegmentList <- function(segments, individuals, regionmap){ 
    FullSegments <- data.frame() 
    for(region in 1:nrow(regionmap)){ 

      for(ind in individuals){ 
       # If there is not a segment within that region for that individual 
       if(nrow(
        segments[segments$start >= regionmap$Start[region] & 
            segments$End <= regionmap$End[region] & 
            segments$SampleID == ind , ] 
       ) == 0){ 
        Temp <- data.frame(SampleID = ind, 
             Start = regionmap$Start[region], 
             End = regionmap$End[region], 
             State = 3 
        ) 
       } 
       # If there is a segment within that region for that individual 
       if(nrow(
        segments[segments$Start >= regionmap$Start[region] & 
            segments$End <= regionmap$End[region] & 
            segments$SampleID == ind , ] 
       ) == 1){ 
        Temp <- data.frame(SampleID = segments$SampleID, 
             Start = regionmap$Start[region], 
             End = regionmap$End[region], 
             State = segments$State[segments$Start >= regionmap$Start[region] & 
                    segments$SampleID == ind ] 
        ) 
       } 
       FullSegments <- list(FullSegments, Temp)    
      } 
    } 
    FullSegments 
} 

、私が見てする必要があります各地域(約53,000)で値域(State、存在しなければ3の値を与える)を各individualの領域に割り当て、すべての個体ごとにすべての領域で新しいdata.frameを作成します。これを行うには、地域と重複していて、それをテーブルに追加するsegment(これらは〜25,000件あります)がある地域と個人をループしています。ここで

は、上記のおもちゃのデータからの出力は与えるものです:

SampleID  Start  End  State 
Ind1   721290  1688192  1 
Ind1   1688193  2926555  2 
Ind2   721290  1688192  2 
Ind2   1688193  2926555  5 
Ind3   721290  1688192  4 
Ind3   1688193  2926555  4 
Ind4   721290  1688192  2 
Ind4   1688193  2926555  6 
それが実行するのに非常に長い時間がかかりますことを除いて、私は、それを必要とする正確にどのような作品であるとして、この関数は、(使用して

system.time、私は実行するために3ヶ月以上かかるだろう)。私はこれを行うためのより良い方法がなければならないことを知っています。私はapply関数を実装しようとしましたが、data.frameの代わりにリストを使うためにいくつかの質問がありました。私はまた、これを単純化するためのdata.tableとplyrオプションがあることを知りました。私はこれらを試しましたが、ifステートメントでネストされたループで動作するようになっていませんでした。

私はこの複合体を何か書いたのはこれが初めてであるため、回答の説明をいただければ幸いです。私が関連していると思います

質問:

ループのネストされた上の他の多くの質問が適用機能(例えばapply(df, 1, function(x){ mean(x) })行うためによく働くやって計算を伴いますしかし、私はそれをdata.frameからdata.frameへの値のマッピングに採用することはできませんでした。

答えて

2

領域とセグメントの開始および終了座標のような「の範囲の整数」上で動作IRanges Bioconductorパッケージ。

source("https://bioconductor.org/biocLite.R") 
biocLite("IRanges") 

読み込み、それをして、パッケージをインストールして、関心の範囲の表現を作成

library(IRanges) 
r <- with(region_map, IRanges(Start, End)) 
s <- with(segments, IRanges(Start, End)) 

結果は、これまであなたが間の重複を見つけることに興味を持っている

> r 
IRanges object with 2 ranges and 0 metadata columns: 
      start  end  width 
     <integer> <integer> <integer> 
    [1] 721290 1688192 966903 
    [2] 1688193 2926555 1238363 
> s 
IRanges object with 9 ranges and 0 metadata columns: 
      start  end  width 
     <integer> <integer> <integer> 
    [1] 721290 1688193 966904 
    [2] 1688194 2926555 1238362 
    [3] 721290 1688190 966901 
    [4] 1688200 2900000 1211801 
    [5] 721290 2926555 2205266 
    [6] 2926600 3000000  73401 
    [7] 721290 1500000 778711 
    [8] 1688193 2005000 316808 
    [9]  690 500000 499311 

です'クエリ'セグメントと '件名' region_map

olaps <- findOverlaps(s, r) 

これは重複の何百万人にも拡大されます

> olaps 
Hits object with 9 hits and 0 metadata columns: 
     queryHits subjectHits 
     <integer> <integer> 
    [1]   1   1 
    [2]   1   2 
    [3]   2   2 
    [4]   3   1 
    [5]   4   2 
    [6]   5   1 
    [7]   5   2 
    [8]   7   1 
    [9]   8   2 
    ------- 
    queryLength: 9/subjectLength: 2 

を与えます。

あなたは、すべての地域のすべての個人の状態に興味があると言いました。あなたのコードからは、地域に属さない個人が状態3を持つように見えます。私は、我々が

idx <- matrix(c(subjectHits(olaps), 
       match(segments$SampleID[queryHits(olaps)], individual)), 
       ncol=2) 

を発見したオーバーラップに基づいてマトリクスに二列のインデックスを作成し、状態を更新するためにインデックス行列を使用する3

state <- matrix(3, nrow(region_map), length(individual), 
       dimnames=list(NULL, individual)) 

すべての状態でマトリックスを作成

state[idx] <- segments$State[queryHits(olaps)] 

これは、実際には、それぞれの領域x個々の組み合わせの状態 - 希望する結果を要約します。 1つの可能性のある問題は、同じ個体の2つのセグメントが単一の領域に重なり、セグメントの状態が異なる場合です。 1つの状態だけが割り当てられます。 data.frameとして、それをキャスト

> state 
    Ind1 Ind2 Ind3 Ind4 
[1,] 1 2 4 2 
[2,] 2 5 4 6 

、例えば、

data.frame(SampleID=colnames(state)[col(state)], 
      Start=region_map[row(state), "Start"], 
      End=region_map[row(state), "End"], 
      State=as.vector(state)) 
+0

これは私のために働くもので、実際のデータを理解して修正することができます。私は染色体情報も持っていたので、私のデータにはGenomicRangesパッケージを使用しなければならなかった。すべてを理解するのにはしばらく時間がかかりましたが、非常に徹底的で役立つ説明に感謝します! –

+0

ああ、私はsystem.timeを使ってこの時間を計りました:user:0.46、system:0.06、elapsed:0.51。かなり素晴らしい。 –

+1

@GaiusAugustusあなたが生産的な一日を過ごしたように聞こえます。あなたの質問がBioconductor関連のものであれば、[Bioconductor support site](https://support.bioconductor.org)に投稿する方が良いでしょう –

0

私はあなたが「この複合体」を必要としないとは思わない。あなたはいくつかの結合であなたがしているすべてを行うことができます。この場合、data.tableを使用します。

あなたは回答の説明を求めましたが、これについてはdata.table homepageの方向に向けるよりもうまくできません。 set*:=コマンドが何を行い、どのように「更新によって参照を」更新するかを理解することが重要です。

データをdata.tableに設定します。

library(data.table) 

dt_individual <- data.table(SampleID = individual) 
dt_region <- data.table(region_map) 
dt_segment <- data.table(segment) 

だけで(すべてのデータを表示するには、重複領域

setkey(dt_join, SampleID, Start, End) 
setkey(dt_segment, SampleID, seg_Start, seg_End) 

foverlaps(dt_join, 
      dt_segment, 
      type="any") 

# SampleID seg_Start seg_End State Start  End 
# 1:  Ind1 721290 1688192  1 721290 1688192 
# 2:  Ind1 1688194 2926555  2 1688193 2926555 
# 3:  Ind2 721290 1688190  2 721290 1688192 
# 4:  Ind2 1688200 2900000  5 1688193 2926555 
# 5:  Ind3 721290 2926555  4 721290 1688192 
# 6:  Ind3 721290 2926555  4 1688193 2926555 
# 7:  Ind4 721290 1500000  2 721290 1688192 
# 8:  Ind4 1688193 2005000  6 1688193 2926555 

を見つけることfoverlaps機能を使用して一緒に

## Change some column names of `dt_segment` so we can identify them after the joins 
setnames(dt_segment, c("Start", "End"), c("seg_Start", "seg_End")) 

## create a 'key_col' to join all the individuals to the regions 
dt_join <- dt_individual[, key_col := 1][ dt_region[, key_col := 1], on="key_col", allow.cartesian=T][, key_col := NULL] 
# SampleID Start  End 
# 1:  Ind1 721290 1688192 
# 2:  Ind2 721290 1688192 
# 3:  Ind3 721290 1688192 
# 4:  Ind4 721290 1688192 
# 5:  Ind1 1688193 2926555 
# 6:  Ind2 1688193 2926555 
# 7:  Ind3 1688193 2926555 
# 8:  Ind4 1688193 2926555 

今それをすべての参加すなわち収まる両方のものを地域内とそうでない地域内で)cartesianを参加させてからassiあなたは

dt_join[dt_segment, on="SampleID", nomatch=0, allow.cartesian=T] 
+0

私はこれで少し混乱しています。 1)あなたはInd3のために4つのアウトプットを持っています。実際のデータでは、すべてのセグメントが> 1の領域に収まるでしょう。2)与えられた値(私のデータでは値= 3)?私はdata.tableパッケージを使用しましたが、この複合体のために決して使用しませんでした。 –

+0

明確にするために、私の出力には、各個人の地域ファイルからその地域内の状態(その地域に属するセグメントによって識別される)を持つ地域ごとに1行があることに注意してください。出力(例:2行目)は重複しない領域とリストされた状態を持ちます。 –

+1

@GaiusAugustus - 私は 'foverlaps'を使う答えを変更しました。 – SymbolixAU

1

を望むようにそれ以外の領域内のものとされたものとGN値は、あなたはnrow(some-subset-of-your-data)を読んで、あなたのコード内の行の多くを持っています。 sum(the-conditions)に切り替えた場合、パフォーマンスが急速に向上します。たとえば:

入れます:

nrow(segments[segments$start >= regionmap$Start[region] & 
            segments$End <= regionmap$End[region] & 
            segments$SampleID == ind , ]) == 0 

この方法

sum(segments$start >= regionmap$Start[region] & 
            segments$End <= regionmap$End[region] & 
            segments$SampleID == ind) == 0 

中に、Rは、メモリ内の自分のサブセット化、データフレームごとに格納されません。

さらに、この操作をブール値として保存するので、ループごとに1回だけ呼び出す必要があります。

isEmpty <- sum(segments$start >= regionmap$Start[region] & 
            segments$End <= regionmap$End[region] & 
            segments$SampleID == ind) == 0 

if(isEmpty){ 
### do something 
} else if(!isEmpty) { 
### do something else 
} 
関連する問題