r-田庄-(左)加入



我有两个Granges对象,我希望将它们合并,以便在两个对象中都不存在元数据的情况下组合两个GRange。

> 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填充,但在我看来这有点过头了,我确信有一个方法存在,但我找不到它:/

谢谢!

使用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)

如果有多个重叠,答案表单@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

附言:我很抱歉将此添加为答案,但由于我的声誉不佳,我无法发表评论。

相关内容

  • 没有找到相关文章

最新更新