r - R中CoverageHeatmap(Bioconductor)函数的问题
问题描述
我有 2 组成对比对,其中查询基因组 1 (q1) 与参考基因组比对,查询基因组 2 (q2) 与同一参考基因组比对。因此,我与参考基因组中的坐标系进行了比对。对齐方式是 Granges 对象的形式。
我想通过在中心对齐 q1 的断点来将 q2 的断点投影到 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:
一组等长的 Grange
track:
指定覆盖范围的 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)
解决方案
为了避免这些问题并达到您的需要,最简单的解决方案是为两个 Grange 设置相同的 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))
推荐阅读
- javascript - 如何使用 Javascript 将菜单连接到显示器?
- ocaml - 如何在 Ocaml 中将列表的元素组合成一个元素?
- c++ - 为什么我的函数以向量数组作为参数没有给出任何值?
- java - 当我尝试使用 setVariable 在类中设置变量时,它没有设置它
- node.js - 每次扫描 fs.watch NodeJS 后保存文件结构?
- postgresql - 如何在 Postgres 中使用带有多个分位数的 percentile_conts
- terraform - 在 vm windows azure 上执行脚本 terraform
- python - 在python中合并json数据
- python - ^ 不支持的操作数类型:'bytes' 和 'bytes' with XOR
- python - 在python中输入数据 - 让它坚持