从 R/Bioconductor 中的 IRanges 对象中提取值答案

【问题标题】：Extracting values from IRanges objects in R/Bioconductor从 R/Bioconductor 中的 IRanges 对象中提取值
【发布时间】：2012-04-12 02:47:39
【问题描述】：

我已使用 import.bw()（来自 rtracklayer 包）将 UCSC 可对齐性轨道导入 R，但无法访问我需要的值。

例如：我想提供一条染色体和一个碱基并返回该位置的值。

我的对象叫 al100：

> al100
RangedData with 21591667 rows and 1 value column across 25 spaces
            space               ranges   |       score
         <factor>            <IRanges>   |   <numeric>
1            chr1       [10001, 10014]   | 0.002777778
2            chr1       [10015, 10015]   | 0.333333343
3            chr1       [10016, 10026]   | 0.500000000
4            chr1       [10027, 10031]   | 1.000000000

我想要一个指定染色体和位置并取回分数的函数。如果我想要一个或两个值，这很简单，但是当我要查找 700 万个值时，循环将不起作用；每次查询 4/5 秒，大约是 10 个月，这不是一个选项。

例如，chr1，位置 10011 将返回值 0.002777778（其中 x 是包含染色体和位置列表的单独对象）

到目前为止，我发现的唯一方法是询问我的位置是否等于或大于范围的起点和/或等于或等于或小于范围的终点。不是很好。

score(al100["chr1"])[ which( start(al100["chr1"]<=x$POS[1])) & end(al100["chr1"]<=x$POS[1]))   ]

【问题讨论】：

底部的代码块是需要 4/5 秒才能运行的代码块吗？该查询中是否存在错误？看起来您正在寻找小于您在 x (end...<=x...) 中指定位置的结尾。括号去掉了吗？或者start 函数真的接受布尔向量吗？

标签： r bioconductor

【解决方案1】：

对于一个可重现的例子

library(rtracklayer)
example(import.bw)
gffRD

给予

> head(gffRD, 3)
RangedData with 3 rows and 7 value columns across 1 space
                                  space       ranges |     type       source
                               <factor>    <IRanges> | <factor>     <factor>
1 Escherichia_coli_K-12_complete_genome [ 337, 2799] |      CDS glimmer/tico
2 Escherichia_coli_K-12_complete_genome [2801, 3733] |      CDS glimmer/tico
3 Escherichia_coli_K-12_complete_genome [3734, 5020] |      CDS glimmer/tico
     phase   strand        note     shift     score
  <factor> <factor> <character> <numeric> <numeric>
1       NA        +          NA        NA  5.347931
2       NA        +          NA        NA 11.448764
3       NA        +          NA        NA  6.230648

定义感兴趣区域

roi <- GRanges("Escherichia_coli_K-12_complete_genome",
               IRanges(c(337, 3734), width=1))

然后使用findOverlaps 在gffRD 和roi 之间进行映射

olaps <- findOverlaps(gffRD,roi)
df <- DataFrame(seqnames=seqnames(roi)[subjectHits(olaps)],
                 start=start(roi)[subjectHits(olaps)],
                 Score=score(gffRD)[queryHits(olaps)])

olaps 包含有关哪些查询匹配哪些主题的信息

> olaps
Hits of length 2
queryLength: 14
subjectLength: 2
  queryHits subjectHits 
   <integer>   <integer> 
 1         1           1 
 2         3           2

数据框是

> df
DataFrame with 2 rows and 3 columns
                               seqnames     start     Score
                                  <Rle> <integer> <numeric>
1 Escherichia_coli_K-12_complete_genome       337  5.347931
2 Escherichia_coli_K-12_complete_genome      3734  6.230648

【讨论】：