我已经使用import.bw() (来自rtracklayer包)将UCSC alignability track导入到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])) ]发布于 2012-03-29 02:20:09
作为一个可重现的例子
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.230648https://stackoverflow.com/questions/9908716
复制相似问题