R 中 for 循环的替代方案？

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

第二个是覆盖单个基地。（染色体、位置、覆盖率）

我现在想获得每个位置旁边写的每个外显子的名称。这将导致一些位置没有 "Exon"，我想在之后删除它们。

我找到了一种方法来做我想做的事。但是它需要 3 个 for 循环和大约 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

R 中 for 循环的替代方案？

Alternatives for for loops in R?

code-coverage

r

bioinformatics