У меня есть 2 набора попарных выравниваний, где геном запроса 1 (q1) выровнен с эталонным геномом, а геном запроса 2 (q2) выровнен с тем же эталонным геномом. Поэтому у меня есть оба выравнивания с системой координат в эталонном геноме. Выравнивания имеют форму объектов GRange.
Я хотел бы спроецировать точки останова q2 на q1, выровняв точки останова q1 в центре, и искать любую кластеризацию точек останова q2 вокруг точек останова q1, все в системе координат эталонного генома.
Поэтому я создаю объект GRange из q1 с точками останова в центре. Например, если есть точка останова в q1 относительно эталонного генома на каркасе 1, п.н. 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
Затем я создаю объект GRanges точек останова на q2, но все длины последовательностей имеют длину 1. Я пересекаю его с объектом q1 GRanges, так что q2 получает только точки, которые можно спроецировать на q1.
Для функции CoverageHeatmap требуется:
windows:
набор GRange равной длины.
track:
объект GRange или 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 в качестве объекта GRange и 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)
