5

私はここに質問を投稿しました：2番目のファイルの範囲に入る1つのファイルの数に基づいて2つのファイルをマージすることについてのRの一致した範囲のマージ。これまでのところ、これを実現するためにコードをつなぎ合わせることに成功していません。私が抱えている問題は、使用しているコードがファイルを1行ずつ比較していることです。これは問題です。1。）1つのファイルが他のファイルよりもはるかに長い、2。）同じ行の範囲だけでなく、長いファイルのすべての範囲ペアをスキャンするために、短いファイルの行が必要です。 。

私は元の質問に投稿された関数を使用しており、最初のファイルのすべての行を2番目のファイルの各行と比較するより一般的なループに適用する方法があるはずですが、私はしていません。 tはまだそれを理解しました。何か提案があれば、よろしくお願いします。

****編集済み。

データの性質は次のとおりです。ほとんどの範囲は一意ですが、各範囲は必ずしも一意ではありません。それらも同じサイズではなく、完全に他の中に収まるものもあります。findIntervalしたがって、「降順ではない」順序に分類するために範囲を並べ替えることができないため、エラーが発生します。

各データフレームの最初の6行は次のとおりです。

file1test <- data.frame(SNP=c("rs2343", "rs211", "rs754", "rs854", "rs343", "rs626"), BP=c(860269, 369640, 861822, 367934, 706940, 717244))


file2 <- data.frame(Gene=c("E613", "E92", "E49", "E3543", "E11", "E233"), BP_start=c(367640, 621059, 721320, 860260, 861322, 879584), BP_end = c(368634, 622053, 722513, 879955, 879533, 894689))

したがって、ご覧のとおり、5行目の範囲は4行目の範囲内にあり、最初のファイルの2つのSNPは4行目の範囲内にありますが、2行目の範囲内にあるのは1つだけです。

SNPを含む最初のファイルには、約400行しかありません。ただし、範囲を含む2番目のファイルには約20Kがあります。出力として生成したいのは、最初のファイル（SNP）の行と、2番目のファイルのBP範囲に含まれるBPを含むデータフレームです。SNPが2つの範囲に分類される場合、2回表示されます。

4

3 に答える 3

7

BioconductorのGenomicRangesパッケージは、このタイプの操作用に設計されています。たとえば、read.delimを使用してデータを読み込み、次のようにします。

con <- textConnection("SNP     BP
rs064   12292
rs319   345367
rs285   700042")
snps <- read.delim(con, head=TRUE, sep="")

con <- textConnection("Gene    BP_start    BP_end
E613    345344      363401
E92     694501      705408
E49     362370      368340") ## missing trailing digit on BP_end??
genes <- read.delim(con, head=TRUE, sep="")

次に、それぞれから「IRanges」を作成します

library(IRanges)
isnps <- with(snps, IRanges(BP, width=1, names=SNP))
igenes <- with(genes, IRanges(BP_start, BP_end, names=Gene))

（座標系に注意してください。IRangesは開始と終了が範囲に含まれることを期待します。また、end = start -1の場合、end> = startは0幅の範囲を期待します）。次に、遺伝子（「subject」）と重複するSNP（IRangesの用語では「query」）を見つけます

olaps <- findOverlaps(isnps, igenes)

2つのSNPが重複しています

> queryHits(olaps)
[1] 2 3

そしてそれらは第1と第2の遺伝子と重複しています

> subjectHits(olaps)
[1] 1 2

クエリが複数の遺伝子と重複している場合は、queryHitsで繰り返されます（その逆も同様です）。次に、データフレームを次のようにマージできます。

> cbind(snps[queryHits(olaps),], genes[subjectHits(olaps),])
    SNP     BP Gene BP_start BP_end
2 rs319 345367 E613   345344 363401
3 rs285 700042  E92   694501 705408

通常、遺伝子とSNPには染色体と鎖（鎖が重要ではないことを示す「+」、「-」、または「*」）情報があり、これらのコンテキストでオーバーラップを実行する必要があります。'IRanges'インスタンスを作成する代わりに、'GRanges'（ゲノム範囲）を作成すると、その後の簿記が自動的に処理されます。

library(GenomicRanges)
isnps <-
    with(snps, GRanges("chrA", IRanges(BP, width=1, names=SNP), "*")
igenes <-
    with(genes, GRanges("chrA", IRanges(BP_start, BP_end, names=Gene), "+"))
于 2012-08-09T23:27:47.367 に答える
7

あなたが求めているのは。だと思いますconditional join。これらはSQLで簡単であり、このsqldfパッケージを使用すると、SQLを使用してRでデータフレームを簡単にクエリできます。

比類のないSNPをどのように処理するかに応じてバージョンを選択するだけです。

内部参加バージョン：

> sqldf("select * from file1test f1 inner join file2 f2 
+       on (f1.BP > f2.BP_start and f1.BP<= f2.BP_end) ")

出力：

     SNP     BP  Gene BP_start BP_end
1 rs2343 860269 E3543   860260 879955
2  rs754 861822 E3543   860260 879955
3  rs754 861822   E11   861322 879533
4  rs854 367934  E613   367640 368634
>

左参加バージョン：

> sqldf("select * from file1test f1 left join file2 f2 
+       on (f1.BP > f2.BP_start and f1.BP<= f2.BP_end) ")

出力：

     SNP     BP  Gene BP_start BP_end
1 rs2343 860269 E3543   860260 879955
2  rs211 369640  <NA>       NA     NA
3  rs754 861822 E3543   860260 879955
4  rs754 861822   E11   861322 879533
5  rs854 367934  E613   367640 368634
6  rs343 706940  <NA>       NA     NA
7  rs626 717244  <NA>       NA     NA
>

=BPがBP_startまたはBP_endと完全に一致する場合に、BPがどのグループに分類されるかが重要な場合は、配置する場所に注意する必要があることに注意してください。

于 2012-08-10T04:09:21.760 に答える
0

とを使用して、ベースの範囲でマージできます。一致するもののみを返すには（内部結合）：applywhich

do.call(rbind, apply(file1test, 1, function(x) {
  if(length(idx <- which(x["BP"] >= file2$BP_start & x["BP"] <= file2$BP_end)) > 0) {
    cbind(rbind(x), file2[idx,])
  }
}))
#      SNP     BP  Gene BP_start BP_end
#x  rs2343 860269 E3543   860260 879955
#1   rs754 861822 E3543   860260 879955
#2   rs754 861822   E11   861322 879533
#x1  rs854 367934  E613   367640 368634

またはfile1testからのすべて（左結合）：

do.call(rbind, apply(file1test, 1, function(x) {
  if(length(idx <- which(x["BP"] >= file2$BP_start & x["BP"] <= file2$BP_end)) > 0) {
    cbind(rbind(x), file2[idx,])
  } else {
    cbind(rbind(x), file2[1,][NA,])
  }
}))
#      SNP     BP  Gene BP_start BP_end
#x  rs2343 860269 E3543   860260 879955
#x1  rs211 369640  <NA>       NA     NA
#1   rs754 861822 E3543   860260 879955
#2   rs754 861822   E11   861322 879533
#x2  rs854 367934  E613   367640 368634
#x3  rs343 706940  <NA>       NA     NA
#x4  rs626 717244  <NA>       NA     NA
于 2019-07-15T13:23:51.830 に答える