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

Question

我有两组成对比对，其中查询基因组1（q1）与参考基因组比对，查询基因组2（q2）与相同的参考基因组比对。因此，我在参考基因组中都与坐标系比对。对齐方式采用GRanges对象的形式。

我想通过对齐中心的q1断点将q2的断点投影到q1上，并在参考基因组坐标系中寻找所有围绕q1断点的q2断点聚类。

因此，我将q1的GRanges对象的断点放在中间。例如，如果相对于参考基因组在1号支架bp处q1有一个断点，则为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对象，但是所有seqlengths的长度都为1。我将其与q1 GRanges对象相交，以便q2仅获得可以投影到q1上的点。

CoverageHeatmap函数需要：

windows：

一组等长的GRanges

track：

指定coverage的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)

我有两组成对比对，其中查询基因组1（q1）与参考基因组比对，查询基因组2（q2）与相同的参考基因组比对。因此，我有两个对齐方式...

Answer 1

为了避免这些问题并实现您所需要的，最简单的解决方案是两个GRanges具有相同的seqlevels和seqlenghts。如果您知道此信息以供参考，请提供它；否则请尝试执行此操作：