1

我有两个 Granges 对象，我希望它们被合并以便合并两个 Granges，即使元数据在两个对象中都不存在。

> t
GRanges object with 2 ranges and 5 metadata columns:
   seqnames         ranges strand |       pvalue      qvalue meth.diff           gc.X  gc.score
      <Rle>      <IRanges>  <Rle> |    <numeric>   <numeric> <numeric>      <GRanges> <numeric>
[1]       MT   [ 708,  708]      + | 2.898639e-04 0.007018699 0.2231039     MT:706-710        80
[2]       MT   [1147, 1147]      - | 6.043240e-05 0.003882324 0.2243177   MT:1146-1150        80
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths



> s
GRanges object with 1 range and 6 metadata columns:
  seqnames     ranges strand |       pvalue      qvalue meth.diff      gc.X     gc.name  gc.score
       <Rle>  <IRanges>  <Rle> |    <numeric>   <numeric> <numeric> <GRanges> <character> <numeric>
[1]       MT [708, 708]      + | 0.0002898639 0.007018699 0.2231039  MT:708:+  rs28412942         0
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths

我想在最后：

> combined
GRanges object with 2 ranges and 5 metadata columns:
   seqnames         ranges strand |       pvalue      qvalue meth.diff           gc.X  gc.score  gc.name
      <Rle>      <IRanges>  <Rle> |    <numeric>   <numeric> <numeric>      <GRanges> <numeric>
[1]       MT   [ 708,  708]      + | 2.898639e-04 0.007018699 0.2231039     MT:706-710        80   rs28412942
[2]       MT   [1147, 1147]      - | 6.043240e-05 0.003882324 0.2243177   MT:1146-1150        80   NA
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths

我尝试的所有方法都要求具有相同数量的列。我可以尝试创建所需的列并填充 NA，但在我看来有点矫枉过正，我确定存在一种方法，但我找不到它：/

谢谢 ！

4

2 回答 2

1

这很容易使用findOverlaps

library(GenomicRanges)  
m <- findOverlaps(s, t)

# Add gc.name to subject GRanges (i.e. left join)
mcols(t)$gc.name <- NA
mcols(t)[subjectHits(m), "gc.name"] = mcols(s)[queryHits(m), "gc.name"]
t
#GRanges object with 2 ranges and 6 metadata columns:
#      seqnames       ranges strand |    pvalue     qvalue meth.diff        gc.X
#         <Rle>    <IRanges>  <Rle> | <numeric>  <numeric> <numeric> <character>
#  [1]       MT [ 708,  708]      + | 0.3361535 0.06058539 0.4743142           a
#  [2]       MT [1147, 1147]      - | 0.4637233 0.19743361 0.3010486           b
#       gc.score     gc.name
#      <numeric> <character>
#  [1]        80  rs28412942
#  [2]        80        <NA>
#  -------
#  seqinfo: 1 sequence from an unspecified genome; no seqlengths

样本数据

set.seed(2018)
t <- GRanges(
    seqnames = c("MT", "MT"),
    IRanges(start = c(708, 1147), end = c(708, 1147)),
    strand = c("+", "-"),
    pvalue = runif(2),
    qvalue = runif(2),
    meth.diff = runif(2),
    gc.X = letters[1:2],
    gc.score = rep(80, 2))

set.seed(2018)
s <- GRanges(
    seqnames = "MT",
    IRanges(start = 708, end = 708),
    strand = "+",
    pvalue = runif(1),
    qvalue = runif(1),
    meth.diff = runif(1),
    gc.X = letters[3],
    gc.name = "rs28412942",
    gc.score = 0)
于 2018-08-06T22:22:42.863 回答
1

如果有多个重叠，@Mautirus Evers的答案表格将失败。findOverlaps正确识别两个重叠，但由于输入查询只有一行，第二部分mcols仅分配第二个重叠。

“失败”测试：

set.seed(2018)
t <- GRanges(
  seqnames = "MT",
  IRanges(start = 708, end = 708),
  strand = "+",
  pvalue = runif(1),
  qvalue = runif(1),
  meth.diff = runif(1),
  gc.X = letters[1],
  gc.score = rep(80, 1))
t

# GRanges object with 1 range and 5 metadata columns:
#       seqnames    ranges strand |            pvalue            qvalue          meth.diff        gc.X  gc.score
#          <Rle> <IRanges>  <Rle> |         <numeric>         <numeric>          <numeric> <character> <numeric>
#   [1]       MT       708      + | 0.336153471143916 0.463723270455375 0.0605853858869523           a        80
#   -------
#   seqinfo: 1 sequence from an unspecified genome; no seqlengths

set.seed(2018)
s <- GRanges(
  seqnames = "MT",
  IRanges(start = c(708, 708), end = c(708, 708)),
  strand = c("+", "+"),
  pvalue = runif(2),
  qvalue = runif(2),
  meth.diff = runif(2),
  gc.X = letters[3:4],
  gc.name = c("rs28412942", "rs28412942dupl"),
  gc.score = 0)
s

# GRanges object with 2 ranges and 6 metadata columns:
#       seqnames    ranges strand |            pvalue             qvalue         meth.diff        gc.X        gc.name  gc.score
#          <Rle> <IRanges>  <Rle> |         <numeric>          <numeric>         <numeric> <character>    <character> <numeric>
#   [1]       MT       708      + | 0.336153471143916 0.0605853858869523 0.474314193474129           c     rs28412942         0
#   [2]       MT       708      + | 0.463723270455375  0.197433612542227 0.301048603840172           d rs28412942dupl         0
#   -------
#   seqinfo: 1 sequence from an unspecified genome; no seqlengths

library(GenomicRanges)  
m <- findOverlaps(s, t, type = "any", select = "all")
m

#Hits object with 2 hits and 0 metadata columns:
#      queryHits subjectHits
#      <integer>   <integer>
#  [1]         1           1
#  [2]         2           1
#  -------
#  queryLength: 2 / subjectLength: 1

# Add gc.name to subject GRanges (i.e. left join)
mcols(t)$gc.name <- NA
mcols(t)[subjectHits(m), "gc.name"] <- mcols(s)[queryHits(m), "gc.name"]
t

# GRanges object with 1 range and 6 metadata columns:
#       seqnames    ranges strand |            pvalue            qvalue          meth.diff        gc.X  gc.score        gc.name
#          <Rle> <IRanges>  <Rle> |         <numeric>         <numeric>          <numeric> <character> <numeric>    <character>
# [1]       MT       708      + | 0.336153471143916 0.463723270455375 0.0605853858869523           a        80 rs28412942dupl
# -------
# seqinfo: 1 sequence from an unspecified genome; no seqlengths

我最终使用plyranges join_overlap_left_directed（或 join_overlap_left），然后过滤我需要保留的任何列。

plyranges 解决方案：

set.seed(2018)
t <- GRanges(
  seqnames = c("MT",  "MT"),
  IRanges(start = c(708, 800), end = c(708, 800)),
  strand = c("+", "-"),
  pvalue = runif(2),
  qvalue = runif(2),
  meth.diff = runif(2),
  gc.X = letters[1:2],
  gc.score = rep(80, 2))

set.seed(2018)
s <- GRanges(
  seqnames = "MT",
  IRanges(start = c(708, 708), end = c(708, 708)),
  strand = c("+", "+"),
  pvalue = runif(2),
  qvalue = runif(2),
  meth.diff = runif(2),
  gc.X = letters[3:4],
  gc.name = c("rs28412942", "rs28412942dupl"),
  gc.score = 0)

library(plyranges)
t <- join_overlap_left_directed(t, s)
t

# GRanges object with 3 ranges and 11 metadata columns:
#       seqnames    ranges strand |          pvalue.x           qvalue.x       meth.diff.x      gc.X.x gc.score.x          pvalue.y           qvalue.y       meth.diff.y      gc.X.y        gc.name gc.score.y
#          <Rle> <IRanges>  <Rle> |         <numeric>          <numeric>         <numeric> <character>  <numeric>         <numeric>          <numeric>         <numeric> <character>    <character>  <numeric>
#   [1]       MT       708      + | 0.336153471143916 0.0605853858869523 0.474314193474129           a         80 0.336153471143916 0.0605853858869523 0.474314193474129           c     rs28412942          0
#   [2]       MT       708      + | 0.336153471143916 0.0605853858869523 0.474314193474129           a         80 0.463723270455375  0.197433612542227 0.301048603840172           d rs28412942dupl          0
#   [3]       MT       800      - | 0.463723270455375  0.197433612542227 0.301048603840172           b         80              <NA>               <NA>              <NA>        <NA>           <NA>       <NA>
#   -------
#   seqinfo: 2 sequences from an unspecified genome; no seqlengths

PS我很抱歉将其添加为答案，但由于我的声誉低下，我无法发表评论。

于 2022-01-12T22:54:28.853 回答