如何将 data.frame 对象强制为 R 中的 Genomic Ranges 对象？答案

【问题标题】：How to coerce data.frame objects to Genomic Ranges objects in R?如何将 data.frame 对象强制为 R 中的 Genomic Ranges 对象？
【发布时间】：2018-08-27 00:00:21
【问题描述】：

我在 R 中有几个床文件作为 data.frame 对象。现在我想按元素找到两个精益床文件之间的重叠。

为了澄清我的问题，我需要在第一个床文件中逐行进行测试（已经在数据框对象中），因此仅将数据框的一行作为查询，然后将其提供给间隔树，其中保存第二床文件的间隔树（但需要先强制 Granges 对象）。

我的数据看起来像：

   idx  chrom     start    End    name        score  p-value
    1   chr1      32727    32817  MACS_peak_1  8.69 1.150748e-11
    2   chr1      52489    52552  MACS_peak_2  4.26 2.347418e-11
    3   chr1      65527    65590  MACS_peak_3  4.19 2.386635e-11
    4   chr1      65773    65904  MACS_peak_4  2.02 4.950495e-11
    5   chr1      66001    66117  MACS_peak_5  5.66 1.766784e-11
    6   chr1     115700   115769  MACS_peak_6 10.30 9.708738e-12
    7   chr1     136389   136452  MACS_peak_7  4.26 2.347418e-11
    8   chr1     235352   235415  MACS_peak_8  4.26 2.347418e-11
    9   chr1     235636   235700  MACS_peak_9  5.66 1.766784e-11
    10  chr1     432895   432958 MACS_peak_10  4.26 2.347418e-11


f1 <- function(bed.1, bed.2){
  query<- GRanges()
  subject = bed.2
  for(i in 1: length(bed.1)){
    query<-bed.1[i]
    o <- GenomicRanges::findOverlaps(query, subject, minoverlap = 2L, algorithm="intervaltree")
    hitfrom_<-query[queryHits(o)]
    hitTo_<-subject[subjectHits(o)]
    pint <-pintersect(hitfrom_, hitTo_)
    return(pint)
  }
}

这是我的代码，如何在 bed.1 中迭代一组 GRanges 对象并调用 findOverlap() 函数来查找重叠的 GRanges 的位置。这段代码没有给我想要的结果。有人帮我吗？？谢谢

【问题讨论】：

你能提供一个你的 2 个“data.frame”的例子吗？
您是否有理由尝试在 R 中执行所有这些操作？ bedtools 旨在查找床文件之间的重叠。
@heathobrien：我知道这个过程是由 bedtools 完成的。我想要做的是在 R 环境中实现相同的东西。有没有最简单的方法在 R 中进行这种分析？？？
@heathobrien：示例床文件，我已附在原始线程中。
我猜你正在从“IRanges”包中寻找类似find[or count]Overlaps(IRanges(df1$start, df1$End), IRanges(df2$start, df2$End)) 的东西？

标签： r bioinformatics genome

【解决方案1】：

我猜你不需要逐行操作。

text1 <- "idx  chrom     start    End    name        score  p-value
1   chr1      32727    32817  MACS_peak_1  8.69 1.150748e-11
2   chr1      52489    52552  MACS_peak_2  4.26 2.347418e-11
3   chr1      65527    65590  MACS_peak_3  4.19 2.386635e-11
4   chr1      65773    65904  MACS_peak_4  2.02 4.950495e-11
5   chr1      66001    66117  MACS_peak_5  5.66 1.766784e-11
6   chr1     115700   115769  MACS_peak_6 10.30 9.708738e-12
7   chr1     136389   136452  MACS_peak_7  4.26 2.347418e-11
8   chr1     235352   235415  MACS_peak_8  4.26 2.347418e-11
9   chr1     235636   235700  MACS_peak_9  5.66 1.766784e-11
10  chr1     432895   432958 MACS_peak_10  4.26 2.347418e-11"

bed1 <- read.table(text=text1, head=T, as.is=T)

library(GenomicRanges)

bed1.gr <- GRanges(bed1$chrom, IRanges(bed1$start, bed1$End))
bed2 <- data.frame(chr=c("chr1", "chr1"),
                   start=c(30000, 130000),
                   end=c(60000, 200000), stringsAsFactors = FALSE)
bed2.gr <- GRanges(bed2$chr, IRanges(bed2$start, bed2$end))
op <- findOverlaps(bed1.gr, bed2.gr)

op.df <- data.frame(que=queryHits(op), sub=subjectHits(op),
                    stringsAsFactors = FALSE)

bed1$que <- 1:nrow(bed1)
bed2$sub <- 1:nrow(bed2)

bed.n <- merge(bed1, op.df, by="que", all=T)
bed.n <- merge(bed.n, bed2, by="sub", all=T)
bed.n$que <- NULL
bed.n$sub <- NULL
bed.n
#    idx chrom start.x    End         name score      p.value  chr start.y   end
# 1    1  chr1   32727  32817  MACS_peak_1  8.69 1.150748e-11 chr1   30000 6e+04
# 2    2  chr1   52489  52552  MACS_peak_2  4.26 2.347418e-11 chr1   30000 6e+04
# 3    7  chr1  136389 136452  MACS_peak_7  4.26 2.347418e-11 chr1  130000 2e+05
# 4    5  chr1   66001  66117  MACS_peak_5  5.66 1.766784e-11 <NA>      NA    NA
# 5    6  chr1  115700 115769  MACS_peak_6 10.30 9.708738e-12 <NA>      NA    NA
# 6    3  chr1   65527  65590  MACS_peak_3  4.19 2.386635e-11 <NA>      NA    NA
# 7    4  chr1   65773  65904  MACS_peak_4  2.02 4.950495e-11 <NA>      NA    NA
# 8    9  chr1  235636 235700  MACS_peak_9  5.66 1.766784e-11 <NA>      NA    NA
# 9   10  chr1  432895 432958 MACS_peak_10  4.26 2.347418e-11 <NA>      NA    NA
# 10   8  chr1  235352 235415  MACS_peak_8  4.26 2.347418e-11 <NA>      NA    NA

【讨论】：

【解决方案2】：

让我们考虑以下可重现的示例：

a <- GRanges(
  seqnames=Rle(c("chr1", "chr2", "chr3", "chr4"), c(3, 2, 1, 2)),
  ranges=IRanges(seq(1, by=9, len=8), seq(7, by=9, len=8)),
  rangeName=letters[seq(1:8)], score=sample(1:20, 8, replace = FALSE))

b <- GRanges(
  seqnames=Rle(c("chr1", "chr2", "chr3","chr4"), c(4, 3, 1, 1)),
  ranges=IRanges(seq(2, by=5, len=9), seq(4, by=5, len=9)),
  rangeName=letters[seq(1:9)], score=sample(1:20, 9, replace = FALSE))

然后，按元素重叠两个 GRange 对象：

ov <- as(findOverlaps(a,b), "List")

ov 将重叠命中索引向量作为压缩整数列表对象返回。

【讨论】：

【解决方案3】：

为此目的有一个名为makeGRangesFromDataFrame的函数：

makeGRangesFromDataFrame(df,
keep.extra.columns=FALSE,
ignore.strand=FALSE,
seqinfo=NULL,
seqnames.field=c("seqnames", "seqname",
"chromosome", "chrom",
"chr", "chromosome_name",
"seqid"),
start.field="start",
end.field=c("end", "stop"),
strand.field="strand",
starts.in.df.are.0based=FALSE)

如果您的起始列的名称不是上述默认值，您必须使用start.field="Start" 或您的列被调用的任何名称来调用它。 end.field、strand.field 等也一样。

【讨论】：

【解决方案4】：

假设bed.1 和bed.2 都是数据帧，那么这将起作用：

library(GenomicRanges)


colnames(bed.1)[c(2,4)] <- c("seqnames", "end") #GRanges insists on these column names
colnames(bed.2) <- colnames(bed.1)

gr1 <- GRanges(bed.1) #will work with the right column names
gr2 <- GRanges(bed.2)

#if you want the set of ranges overlapped by at least one range from each:
int <- intersect(gr1, gr2) #GRanges output
df.int = as.data.frame(int) #data frame of new ranges

【讨论】：

【解决方案5】：

另一种做法是直接从bed文件本身开始，通过rtracklayer的导入功能直接导入到GRanges对象中：

source("https://bioconductor.org/biocLite.R")
biocLite("rtracklayer")
library(rtracklayer)
import("example.bed", format="bed")

这为您提供了您的 GRanges 对象，而无需通过数据帧中间。

【讨论】：