2017-09-11 8 views
1

私のタイトルは私が望むほど明確ではありません。あなたはvariantTable[i,]$nonAFR_weightedは、複数の列(EAS_MAF、AMR_MAF、AFR_MAF、EUR_MAF、SAS_MAF)間の条件に基づいて計算されて見ることができるようにR:お互いに依存する複数の条件を使用する場合のループの回避

for(i in 1:length(variantTable[,1])){ 
    #N stores counts of numbers of population totals of populations that contain the variant in question 
    N = 0 
    #NF.pop stores frequencies 
    NF.EAS = 0 
    NF.AMR = 0 
    NF.EUR = 0 
    NF.SAS = 0 
    if(variantTable[i,]$EAS_MAF > 0){ 
    NF.EAS = EAScount * variantTable[i,]$EAS_MAF 
    N = N + EAScount 
    } 
    if(variantTable[i,]$AMR_MAF > 0){ 
    NF.AMR = AMRcount * variantTable[i,]$AMR_MAF 
    N = N + AMRcount 
    } 
    if(variantTable[i,]$EUR_MAF > 0){ 
    NF.EUR = EURcount * variantTable[i,]$EUR_MAF 
    N = N + EURcount 
    } 
    if(variantTable[i,]$SAS_MAF > 0){ 
    NF.SAS = SAScount * variantTable[i,]$SAS_MAF 
    N = N + SAScount 
    } 
    variantTable[i,]$nonAFR_N <- N 

    variantTable[i,]$nonAFR_weighted <- (NF.EAS + NF.AMR + NF.EUR + NF.SAS)/N 
} 

:以下は、ループを使用して書かれたいくつかのコードです。

私は、ループが、特に私のデータセットが900000行で構成されているという事実を考慮して、Rでこれを行う最も速い方法ではないことを知っています。

私はifelseを使って作業を始めましたが、メソッドを適用しましたが、このような状況でそれらを使用する方法がわかりません。私は単純に1つの行を取り、その行の値を計算し、次にapplyメソッドを使用する関数を作成しようとしましたが、入力が何であるかわからないのでこれは機能しませんでした。

このような問題を最善に進めるにはどのようなアドバイスが必要ですか?

編集:ここに私のデータのdputがある:

> dput(head(variantTable)) 
structure(list(CHROM = c("1", "1", "1", "1", "1", "1"), POS = c(69224L, 
69428L, 69486L, 69487L, 69496L, 69521L), ID = c("rs568964432", 
"rs140739101", "rs548369610", "rs568226429", "rs150690004", "rs553724620" 
), REF = c("A", "T", "C", "G", "G", "T"), ALT = c("T", "G", "T", 
"A", "A", "A"), AF = c(0.000399361, 0.0189696, 0.000199681, 0.000399361, 
0.000998403, 0.000399361), AC = c(2L, 95L, 1L, 2L, 5L, 2L), AN = c(5008L, 
5008L, 5008L, 5008L, 5008L, 5008L), EAS_AF = c(0, 0.003, 0.001, 
0, 0, 0), AMR_AF = c(0.0029, 0.036, 0, 0, 0.0014, 0.0029), AFR_AF = c(0, 
0.0015, 0, 0.0015, 0.003, 0), EUR_AF = c(0, 0.0497, 0, 0, 0, 
0), SAS_AF = c(0, 0.0153, 0, 0, 0, 0), consequence = c("nonsynonymous SNV", 
"nonsynonymous SNV", "synonymous SNV", "nonsynonymous SNV", "nonsynonymous SNV", 
"nonsynonymous SNV"), gene = c("OR4F5", "OR4F5", "OR4F5", "OR4F5", 
"OR4F5", "OR4F5"), refGene_id = c("NM_001005484", "NM_001005484", 
"NM_001005484", "NM_001005484", "NM_001005484", "NM_001005484" 
), AA_change = c("('D', 'V')", "('F', 'C')", "('N', 'N')", "('A', 'T')", 
"('G', 'S')", "('I', 'N')"), X0.fold_count = c(572L, 572L, 572L, 
572L, 572L, 572L), X4.fold_count = c(141L, 141L, 141L, 141L, 
141L, 141L), EAS_MAF = c(0, 0.003, 0.001, 0, 0, 0), AMR_MAF = c(0.0029, 
0.036, 0, 0, 0.0014, 0.0029), AFR_MAF = c(0, 0.0015, 0, 0.0015, 
0.003, 0), EUR_MAF = c(0, 0.0497, 0, 0, 0, 0), SAS_MAF = c(0, 
0.0153, 0, 0, 0, 0), nonAFR_AF = c(0.0029, 0.104, 0.001, 0, 0.0014, 
0.0029), nonAFR_N = c(309227, 1128036, 262551, 0, 309227, 309227 
), nonAFR_weighted = c(0.0029, 0.0261704282487438, 0.001, NaN, 
0.0014, 0.0029)), .Names = c("CHROM", "POS", "ID", "REF", "ALT", 
"AF", "AC", "AN", "EAS_AF", "AMR_AF", "AFR_AF", "EUR_AF", "SAS_AF", 
"consequence", "gene", "refGene_id", "AA_change", "X0.fold_count", 
"X4.fold_count", "EAS_MAF", "AMR_MAF", "AFR_MAF", "EUR_MAF", 
"SAS_MAF", "nonAFR_AF", "nonAFR_N", "nonAFR_weighted"), row.names = c(NA, 
6L), class = "data.frame") 

集団数(EAScount、AMRcountなど)のような、以前に定義されています

EAScount <- length(variantTable$EAS_MAF[variantTable$EAS_MAF>0]) 
AMRcount <- length(variantTable$EAS_MAF[variantTable$AMR_MAF>0]) 
AFRcount <- length(variantTable$EAS_MAF[variantTable$AFR_MAF>0]) 
EURcount <- length(variantTable$EAS_MAF[variantTable$EUR_MAF>0]) 
SAScount <- length(variantTable$EAS_MAF[variantTable$SAS_MAF>0]) 

私が探しています出力variantTable $ nonAFR_nとvariantTable $ nonAFR_weightedの計算です。正しい計算では、以下のような例があります。

> variantTable[2,] 
    CHROM POS   ID REF ALT  AF AC AN EAS_AF AMR_AF AFR_AF EUR_AF SAS_AF  consequence 
2  1 69428 rs140739101 T G 0.0189696 95 5008 0.003 0.036 0.0015 0.0497 0.0153 nonsynonymous SNV 
    gene refGene_id AA_change X0.fold_count X4.fold_count EAS_MAF AMR_MAF AFR_MAF EUR_MAF SAS_MAF 
2 OR4F5 NM_001005484 ('F', 'C')   572   141 0.003 0.036 0.0015 0.0497 0.0153 
    nonAFR_AF nonAFR_N nonAFR_weighted 
2  0.104 1128036  0.02617043 
+0

あなたは(他には何もしないので、私は、あなたがifelseを必要とわからない)このために適用することができます。あなたはデータの再現可能な例を提供できますか?あなたはちょうど 'dput(head(variantTable))'と結果を投稿することができます – csgroen

+0

こんにちは、私はdputを含むように私の投稿を編集しました。ありがとう。 – spiral01

+0

心配はいりません。また、ループの前にEAScount、AMRcount、EURcount、SAScountが定義されていますか? – csgroen

答えて

4

これは動作しますか?

library(dplyr) 
variantTable %>% mutate(
    NF.EAS = EAScount * EAS_MAF, 
    NF.AMR = AMRcount * AMR_MAF, 
    NF.EUR = EURcount * EUR_MAF, 
    NF.SAS = SAScount * SAS_MAF, 
    nonAFR_N = EAScount * (EAS_MAF>0) + AMRcount * (AMR_MAF>0) + EURcount * (EUR_MAF>0) + SAScount * (SAS_MAF>0), 
    nonAFR_weighted = (NF.EAS + NF.AMR + NF.EUR + NF.SAS)/nonAFR_N) %>% 
    select(-c(NF.EAS,NF.AMR,NF.EUR,NF.SAS)) 

のmutateが追加またはテーブルに列を変更し、それが$表記せずに列名の使用を可能にします。あなたのif構造体は必要ありません。なぜなら、ゼロを乗算すれば、デフォルトのゼロ値が得られるからです。そのため、スクリプトの最初の部分を単純化することができます。算術演算子と一緒に使用する場合

ブールは整数01に強制されているので、私はどちらかNを計算するために、すべてのif構造を必要としませんでした。

最後の列はかなり簡単です。

これらの操作はすべてベクトル化されています。これは、離散値ではなく直接加算、減算、掛け算、分割された列を意味します。

また、より効率的で読みやすく:

EAScount <- sum(variantTable$EAS_MAF>0) 
AMRcount <- sum(variantTable$AMR_MAF>0) 
AFRcount <- sum(variantTable$AFR_MAF>0) 
EURcount <- sum(variantTable$EUR_MAF>0) 
SAScount <- sum(variantTable$SAS_MAF>0) 
+0

ありがとうございます、これは非常に迅速に機能しました。私はあなたの答えで何が起こっているのか完全にはわからないのですか?特に%>%? – spiral01

+0

これはpipe演算子と呼ばれ、 'dplyr'のチュートリアルを見つけて、ウェブサイトでいくつかの質問を読んでみると、Rで作業を続ける予定がある場合には時間がかかるでしょう。少し解明するために私の答えを編集します –

+0

私はそれが今より少し明確になることを願っています –

関連する問題