R 中 CoverageHeatmap (Bioconductor) 函数的问题

Question

我有 2 组成对比对，其中查询基因组 1 (q1) 与参考基因组比对，查询基因组 2 (q2) 与同一参考基因组比对。因此，我与参考基因组中的坐标系进行了两种比对。路线采用 GRanges 对象的形式。

我想将 q2 的断点投影到 q1，方法是在中心对齐 q1 的断点，并在 q1 断点周围寻找任何 q2 断点的聚类，所有这些都在参考基因组坐标系中。

因此，我创建了一个 q1 的 GRanges 对象，其断点位于中心。例如，如果 q1 中有一个断点，相对于脚手架 1 的参考基因组，bp 833，然后在其两侧的 500 上取一个 window，q1 GRanges 对象将有一个元素：

GRanges object with 1 range and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]       S1  333-1333      *
  -------
  seqinfo: 576 sequences from an unspecified genome; no seqlengths

然后我在q2上构造一个断点的GRanges对象，但是所有seqlengths的长度都是1。我把它和q1的GRanges对象相交，这样q2只得到可以投影到q1上的点。

CoverageHeatmap 函数需要：

windows:
一组等长的 GRanges

曲目：
指定覆盖率的 GRanges 或 RleList 对象

当我调用 CoverageHeatmap 函数时，我总是收到此错误和警告消息：

Error: subscript contains out-of-bounds ranges
In addition: Warning message:
In e1 == Rle(e2) :
  longer object length is not a multiple of shorter object length
Called from: S4Vectors:::.subscript_error("subscript contains out-of-bounds ", 
    "ranges")

我尝试了很多方法来尝试使这项工作有效，但仍然收到相同的错误和警告消息。这是我的代码（包括当我尝试使用 q2 作为 GRanges 对象和 RleList 的函数时）

## BP Pairwise comparison, using 3rd genome as co-ordinate reference
# q1 is used as the centre point reference, with q2 bps projected on to it. 
# gr_ref_q1 is the pw alignment between the reference and query genome 1
# gr_ref_q2 is the pw alignment between the reference and query genome 2
# We construct two GRanges objects to feed into CoverageHeatMaps
library(schoolmath)
library(heatmaps)
library(IRanges)

bp_3gen_v2 <- function(gr_ref_q1, gr_ref_q2, win){

  # Failsafes (check ref genome is the same, etc)
  if(!(is.even(win))){stop("win should be an even number")}

  ## Construct g1_rco (1st GRanges object)
  # IRanges object
  q1_starts1 <- start(ranges(gr_ref_q1)) - (win*0.5)
  q1_starts2 <- end(ranges(gr_ref_q1)) - (win*0.5)
  q1_starts <- c(q1_starts1, q1_starts2)
  q1_ends1 <- start(ranges(gr_ref_q1)) + (win*0.5)
  q1_ends2 <- end(ranges(gr_ref_q1)) + (win*0.5)
  q1_ends <- c(q1_ends1, q1_ends2)
  q1_ir_ob <- IRanges(start = q1_starts, end = q1_ends) 

  # GR object
  g1_vec_seq <- as.vector(seqnames(gr_ref_q1))
  gr1_seqnames <- c(g1_vec_seq, g1_vec_seq)
  g1_rco <- GRanges(seqnames = gr1_seqnames, ranges = q1_ir_ob, 
                    seqinfo = seqinfo(gr_ref_q1))

  # Remove negative ranges from GR object
  g1_rco <- g1_rco[!(start(ranges(g1_rco)) < 0)] 



  ## Construct g2_rco (2nd GRanges object)
  # IRanges object
  q2_starts <- start(ranges(gr_ref_q2))
  q2_ends <- end(ranges(gr_ref_q2))
  q2_bps <- c(q2_starts, q2_ends)
  q2_ir_ob <- IRanges(start = q2_bps, end = q2_bps)
  # GR object
  g2_vec_seq <- as.vector(seqnames(gr_ref_q2))
  gr2_seqnames <- c(g2_vec_seq, g2_vec_seq)
  g2_rco <- GRanges(seqnames = gr2_seqnames, ranges = q2_ir_ob,
                    seqinfo = seqinfo(gr_ref_q2))

  # Try removing anywhere in g2_rco that is not present in g1_rco
  # find intersection of seqnames 
  g_inter <- intersect(g1_vec_seq, g2_vec_seq)
  # apply to g2_rco to remove out of bound scaffols
  g2_rco <- g2_rco[seqnames(g2_rco) == g_inter]
  # now to remove out of bound ranges (GRanges object)
  g2_red <- intersect(g1_rco, g2_rco)
  # And try as RleList object
  g2_red_rle <- coverage(g2_red)


  # Heatmap
  heat_map <- CoverageHeatmap(windows = g1_rco, track = g2_red_rle)

Answer 1

为了避免这些问题并实现您的需要，最简单的解决方案是让两个 GRanges 具有相同的 seqlevels 和 seqlenghts。如果您知道这个供您参考，请提供它，如果不知道，请尝试这个：

第一个示例数据集：

library(heatmaps)
gr1 = GRanges(seqnames=c(1,2,3),
IRanges(start=c(1,101,1001),end=c(500,600,1500)))

gr2 = GRanges(seqnames=c(2,2,3,3),
IRanges(start=c(1,301,1,1201),end=c(2500,4800,3500,9700)))

然后我们做一个组合范围来得到水平和长度：

combined= range(c(gr1,gr2))

seqlevels(gr1) = as.character(seqnames(combined))
seqlevels(gr2) = as.character(seqnames(combined))
seqlengths(gr1) = end(combined)
seqlengths(gr2) = end(combined)

然后可以通过以下方式轻松获得热图：

CoverageHeatmap(gr1,coverage(gr2))

或者，如果您只想查看在 gr2 中具有某些值的 gr1 windows，则执行：

CoverageHeatmap(gr1[countOverlaps(gr1,gr2)>0],coverage(gr2))

R 中 CoverageHeatmap (Bioconductor) 函数的问题

Problem with CoverageHeatmap (Bioconductor) function in R

r

bioinformatics

heatmap

bioconductor

genomicranges