3

ここに省略されたデータセットがあります： 

SNP chr       BP log10   PPA
rs10068  17 56555 1.16303 0.030
rs10032  17 56561 26.364 0.975
rs10354  17 34951 4.3212 0.626
rs10043  17 20491 0.00097 0.006
rs10457  17 69572 -0.38403 0.014
rs10465  17 69872 8.19547 0.927

ここで、PPAは関連の事後確率です。私はいくつかの高いlog10値（> 6）を持っているので、これらの領域の周りの信頼区間を決定して、それらがどれだけ大きいか小さいかを自信を持って決定したいと思います。

これを行うには、まず、log10> 6のSNPを特定します。これは、サブセットを使用すると十分に単純です。

newdata <- subset(data, log10 > 6)

ただし、BP 500 +/-リードSNPのBP（log10> 6）を使用して、これらのリードSNPに物理的に近いSNPもこのサブセットに含めたいと思います。ここが私が進むための最良の方法がわからないところです。これは私が取り組むことができるものですか、subsetそれとも最初に元のデータでこれらのリードSNPを特定し、次にそこからサブセット化する必要がありますか？

これらの領域を分離したら、先に進むことができます。

任意の提案をいただければ幸いです。

4

2 に答える 2

6 
s <- read.table(header=T, text="SNP chr       BP log10   PPA
rs10068  17 56555 1.16303 0.030
rs10032  17 56561 26.364 0.975
rs10354  17 34951 4.3212 0.626
rs10043  17 20491 0.00097 0.006
rs10457  17 69572 -0.38403 0.014
rs10465  17 69872 8.19547 0.927")

s $ log10> 6である任意の行までの距離：

outer(s$BP[s$log10 > 6], s$BP, '-')
##       [,1]  [,2]  [,3]  [,4]   [,5]   [,6]
## [1,]     6     0 21610 36070 -13011 -13311
## [2,] 13317 13311 34921 49381    300      0

絶対値が500未満の上記の列は、保持したい列です。

s[apply(outer(s$BP[s$log10 > 6], s$BP, '-'), 2, function(x) any(abs(x) < 500)),]
##       SNP chr    BP    log10   PPA
## 1 rs10068  17 56555  1.16303 0.030
## 2 rs10032  17 56561 26.36400 0.975
## 5 rs10457  17 69572 -0.38403 0.014
## 6 rs10465  17 69872  8.19547 0.927
于 2013-02-18T00:00:16.383 に答える
1

価値のあるものとして、GenomicRangesパッケージを使用したソリューションがあります。このパッケージは、IRangesを使用して間隔データを非常に効率的に構築および処理しますinterval trees。そして、これはすべての異なる染色体も処理できるはずです。

  • 最初のロードデータ：

    # load data
df <- read.table(header=T, text="SNP chr       BP log10   PPA
rs10068  17 56555 1.16303 0.030
rs10032  17 56561 26.364 0.975
rs10354  17 34951 4.3212 0.626
rs10043  17 20491 0.00097 0.006
rs10457  17 69572 -0.38403 0.014
rs10465  17 69872 8.19547 0.927")

  • log10>6のサブセットを取得します

    df.f <- df[df$log10 > 6, ]

  • 初期データとサブセットに範囲をロードして作成します

    require(GenomicRanges)
df.gr <- GRanges(Rle(df$chr), IRanges(df$BP, df$BP))
df.f.gr <- GRanges(Rle(df.f$chr), IRanges(df.f$BP, df.f$BP))

  • log10>6のデータとギャップ=+/-500の他のすべてのSNPのすべての重複を検索します

    olaps <- findOverlaps(df.f.gr, df.gr, maxgap=500)

  • 出力を取得

    # if you just want a whole data.frame with ALL SNPs that have 
# log > 6 or within 500 bases of one with log > 6, then,        
out <- df[subjectHits(olaps), ] 

#       SNP chr    BP    log10   PPA
# 1 rs10068  17 56555  1.16303 0.030
# 2 rs10032  17 56561 26.36400 0.975
# 5 rs10457  17 69572 -0.38403 0.014
# 6 rs10465  17 69872  8.19547 0.927

# In case you want for each SNP that has log > 6, all of the 
# snps that are within 500 bases (either side) apart for 
# this SNP, separately, then,
out.list <- split(out, df$BP[queryHits(olaps)])
# $`56555`
#       SNP chr    BP    log10   PPA
# 1 rs10068  17 56555  1.16303 0.030
# 2 rs10032  17 56561 26.36400 0.975
#     
# $`56561`
#       SNP chr    BP    log10   PPA
# 5 rs10457  17 69572 -0.38403 0.014
# 6 rs10465  17 69872  8.19547 0.927
于 2013-02-18T08:51:48.390 に答える