Альтернативы для петли в R?

Question

Я получил 2 файла, которые я хотел бы объединить, используя R.

head(bed)
chr8 41513235 41513282 ANK1.Exon1
chr8 41518973 41519092 ANK1.Exon2

Первый дает интервалы и их имена. (Хромосома, от, до, имя)

head(coverage)
chr1 41513235 20
chr1 41513236 19
chr1 41513237 19

Второй дает покрытия для отдельных баз. (Хромосома, положение, охват)

Теперь я хочу получить имя каждого Экзона, написанное рядом с каждой Позицией. Это приведет к некоторым позициям без "Экзона", которые я хочу удалить впоследствии.

Я разобрался, как делать то, что я хочу. Однако для циклов требуется 3 и около 15 часов вычислительного времени. Поскольку циклы for не лучшая практика в R, я хотел бы знать, если кто-нибудь знает лучший способ, чем:

coverage <- cbind(coverage, "Exon")
coverage[,4] <- NA

for(i in 1:nrow(bed)){
 for(n in bed[i,2]:bed[i,3]{
  for(m in 1:nrow(coverage)){
   if(coverage[m,2]==n){
    file[m,4] <- bed[i,4]
   }
  }
 }
}

na.omit(coverage)

Поскольку все три позиции находятся в интервале "ANK1.Exon1", вывод должен выглядеть следующим образом:

head(coverage) 
chr1 41513235 20 ANK1.Exon1 
chr1 41513236 19 ANK1.Exon1 
chr1 41513237 19 ANK1.Exon1 

Answer 1

Самый быстрый способ выполнить то, что я искал, был:

library("sqldf")
res <- sqldf("select * from coverage f1 inner join bed f2
on(f1.position >=f2.'from' and f1.position <=f2.'to')")

Время расчета сократилось до секунд. Чтобы получить точный результат, как указано выше, датафрейм был дополнительно уменьшен.

res <- cbind(res[1:4],res[8])

Спасибо за вашу помощь.

Редактировать: для больших наборов данных, в которых могут быть одинаковые позиции в более чем одной хромосоме, полезно добавить:

res <- sqldf("select * from coverage f1 inner join bed f2
on(f1.position >=f2.'from' and f1.position <=f2.'to' and f1.Chromosome = f2.Chromosome)")

Answer 2

Этот алгоритм является линейным, если bed а также coverage входы отсортированы, и bed вход не перекрывает целые

> coverage <- read.table("coverage")
> bed <- read.table("bed")
> 
> coverage <- cbind(coverage, "Exon")
> coverage[,4] <- NA
> 
> i_coverage <- 1
> i_bed <- 1
> 
> while(i_coverage <= length(coverage[,1]) && i_bed <= length(bed[,1])) {
+   if(coverage[i_coverage, 2] < bed[i_bed, 2]){
+     i_coverage <- i_coverage + 1
+   }else{
+     #then coverage[i_coverage, 2] >= bed[i_bed, 2]
+     if(coverage[i_coverage, 2] <= bed[i_bed, 3]){
+       coverage[i_coverage,4] <- as.character(bed[i_bed, 4])
+       i_coverage <- i_coverage + 1
+     }else{
+       i_bed <- i_bed + 1
+     }
+   }
+ }

ты получаешь:

> print(coverage)
V1       V2 V3     "Exon"
1 chr1 41513235 20 ANK1.Exon1
2 chr1 41513236 19 ANK1.Exon1
3 chr1 41513237 19 ANK1.Exon1

Answer 3

Использование GenomicRanges:

library("GenomicRanges")

#data
x1 <- read.table(text="chr1 41513235 41513282 ANK1.Exon1
chr1 41518973 41519092 ANK1.Exon2")

x2 <- read.table(text="chr1 41513235 20
chr1 41513236 19
chr1 41513237 19")

#Convert to Granges object:
g1 <-  GRanges(seqnames=x1$V1,
               IRanges(start=x1$V2,
                       end=x1$V3),
               Exon=x1$V4)


g2 <-  GRanges(seqnames=x2$V1,
               IRanges(start=x2$V2,
                       end=x2$V2),
               covN=x2$V3)
#merge
mergeByOverlaps(g1,g2)

#output
# DataFrame with 3 rows and 4 columns
#                            g1       Exon                          g2      covN
#                     <GRanges>   <factor>                   <GRanges> <integer>
# 1 chr1:*:[41513235, 41513282] ANK1.Exon1 chr1:*:[41513235, 41513235]        20
# 2 chr1:*:[41513235, 41513282] ANK1.Exon1 chr1:*:[41513236, 41513236]        19
# 3 chr1:*:[41513235, 41513282] ANK1.Exon1 chr1:*:[41513237, 41513237]        19