Извлечение значений из объектов IRanges в R/Bioconductor

Я импортировал дорожку выравнивания UCSC в R с помощью import.bw() (из пакета rtracklayer), но у меня возникли проблемы с доступом к нужным мне значениям.

Например: я хочу предоставить хромосому и основание и вернуть значение в этой позиции.

Мой объект называется 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

Мне нужна функция, в которой я указываю хромосому и позицию и возвращаю счет. Это тривиально, если мне нужно одно или два значения, но цикл не будет работать, когда у меня есть 7 миллионов для поиска; при 4/5 секундах на запрос это около 10 месяцев, что не вариант.

Например, chr1, позиция 10011 вернет значение 0,002777778 (где x — отдельный объект, содержащий список хромосом и позиций).

Единственный метод, который я нашел до сих пор, - это спросить, равна ли моя позиция началу или больше и/или равна или равна или меньше конца диапазона. Не очень хорошо.

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

r bioconductor
person CPW    schedule 28.03.2012    source источник
comment
Является ли блок кода внизу тем, который выполняется за 4/5 секунды? И есть ли ошибка в этом запросе? Похоже, вы ищете конец, который меньше указанной вами позиции в x (end...<=x...). А скобки сняты? Или функция start действительно принимает логический вектор?   -  person Jeff Allen    schedule 28.03.2012

Ответы (1)


arrow_upward
1
arrow_downward

Для воспроизводимого примера

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
person Martin Morgan    schedule 28.03.2012