R中的CoverageHeatmap（生物导体）功能问题

Question

我有 2 组成对比对，其中查询基因组 1 (q1) 与参考基因组比对，查询基因组 2 (q2) 与同一参考基因组比对。因此，我与参考基因组中的坐标系进行了比对。对齐方式是 Granges 对象的形式。

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

因此，我制作了一个 q1 的 GRanges 对象，其断点位于中心。例如，如果 q1 相对于参考基因组在支架 1 处 bp 833 有一个断点，然后在其任一侧的 500 处取一个窗口，则 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 对象，但所有 seqlength 的长度均为 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 设置相同的 seqlevel 和 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 窗口，则执行以下操作：

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