R中的CoverageHeatmap (生物导体)函数问题

Question

我有两组成对比对，其中查询基因组1 (q1)与参考基因组对齐，查询基因组2 (q2)与同一参考基因组对齐。因此，我和参考基因组中的坐标系都是对齐的。对齐以GRanges对象的形式出现。

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

因此，我使GRanges对象的q1与它的断点在中心。例如，如果q1中有一个相对于参考基因组的断点位于脚手架1，bp 833，那么在其中的任何一边取一个窗口，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

然后，我在GRanges上构造了断点的一个q2对象，但是所有的长度都是1，我将其与q1 GRanges对象相交，这样q2只能获得可以投影到q1上的点。

CoverageHeatmap函数要求：

窗口：

一组等长的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))

或者，如果您只想查看在gr1中具有某些值的gr2窗口，请执行以下操作：

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