Нахождение локальных координат в линейных интервалах
У меня есть data.frame
с interval
координаты (экзонов и экзонов кодирующих последовательностей, транскриптов генов):
df <- data.frame(seqnames=rep("chr1",24),
start=c(3670552,3670552,3421702,3421702,3214482,3216025,4807823,4807914,4808455,4808455,4828584,4828584,4830268,4830268,4832311,4832311,4837001,4837001,4839387,4839387,4840956,4840956,4844963,4844963),
end=c(3671498,3671348,3421901,3421901,3216968,3216968,4807982,4807982,4808486,4808486,4828649,4828649,4830315,4830315,4832381,4832381,4837074,4837074,4839488,4839488,4841132,4841132,4846739,4845013),
strand=c(rep("-",6),rep("+",18)),
exon_number=c(1,1,2,2,3,3,1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9),
width=c(947,797,200,200,2487,944,160,69,32,32,66,66,48,48,71,71,74,74,102,102,177,177,1777,51),
type=c('exon','CDS','exon','CDS','exon','CDS','exon','CDS','exon','CDS','exon','CDS','exon','CDS','exon','CDS','exon','CDS','exon','CDS','exon','CDS','exon','CDS'),
#exon_number=c(1,1,2,2,3,3,1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8),
transcript_id=c(rep("a",6),rep("b",18)),stringsAsFactors = F)
> head(df)
seqnames start end strand exon_number width type transcript_id
1 chr1 3670552 3671498 - 1 947 exon a
2 chr1 3670552 3671348 - 1 797 CDS a
3 chr1 3421702 3421901 - 2 200 exon a
4 chr1 3421702 3421901 - 2 200 CDS a
5 chr1 3214482 3216968 - 3 2487 exon a
6 chr1 3216025 3216968 - 3 944 CDS a
Его организация выглядит следующим образом:
каждый transcript_id
состоит из набора exon
s, которые являются линейными интервалами вдоль seqnames
(которые являются хромосомами), определяемые как start
а также end
, Эти exon
имеют ориентацию / направленность, которая зависит от strand
(что одинаково для всех exon
с каждого transcript_id
): если strand == "-"
тогда на самом деле end
это первая позиция экзона и start
это его последняя позиция. Если strand == "+"
затем start
а также end
являются первой и последней позициями соответственно.
CDS
линии являются подмножеством exon
s. Обычно оставляют первые и последние несколько exon
в каждом transcript_id
каждый exon
имеет идентичный CDS
интервал (с точки зрения координат). Тем не менее, исключения:
CDS
интервалы являются подмножествомexon
s означает, что они могут начинаться и / или заканчиваться в течениеexon
(Ы). Вполне возможно, что первая позиция первойCDS
интервал после того из ближайшегоexon
и / или последняя позиция последнейCDS
интервал до того из ближайшегоexon
, Также возможно, чтоtranscript_id
будет иметь одинexon
(и, следовательно, одинCDS
отвечая этим определениям.transcript_id
где всеexon
с идентичнымиCDS
интервалыtranscript_id
гдеexon
ы, которые не имеют идентичныеCDS
интервалы только первыеtranscript_id
гдеexon
ы, которые не имеют идентичныеCDS
интервалы только последние
Я ищу быстро function
который вернет локальные координаты CDS
относительно координат transcript_id
чья exon
с были объединены. По существу, start
в результате data.frame
это сумма exon
ширина интервала до первого CDS
интервал + start
из первых CDS
- start
из exon
это подмножество (или end
если strand == "-"
) а также end
является CDS
start
+ сумма CDS
ширина.
Вот что я делаю до сих пор, но это медленно:
res.df <- do.call(rbind,lapply(unique(df$transcript_id),function(t){
transcript.df <- df %>% dplyr::filter(transcript_id == t)
up.cds.df <- dplyr::filter(transcript.df,type == "exon",exon_number < min(dplyr::filter(transcript.df,type == "CDS")$exon_number))
down.cds.df <- dplyr::filter(transcript.df,type == "exon",exon_number > max(dplyr::filter(transcript.df,type == "CDS")$exon_number))
cds.df <- dplyr::filter(transcript.df,exon_number >= min(dplyr::filter(transcript.df,type == "CDS")$exon_number,exon_number <= max(dplyr::filter(transcript.df,type == "CDS")$exon_number)))
first.cds.exon.number <- (cds.df %>% dplyr::filter(type == "CDS") %>% dplyr::filter(exon_number == min(exon_number)))$exon_number
if(transcript.df$strand[1] == "+"){
cds.start <- dplyr::filter(cds.df,type == "CDS",exon_number == first.cds.exon.number)$start-dplyr::filter(cds.df,type == "exon",exon_number == first.cds.exon.number)$start+1+sum(up.utr.df$width)
} else{
cds.start <- dplyr::filter(cds.df,type == "exon",exon_number == first.cds.exon.number)$end-dplyr::filter(cds.df,type == "CDS",exon_number == first.cds.exon.number)$end+1+sum(up.utr.df$width)
}
return(data.frame(transcript_id=t,start=cds.start,end=cds.start+sum(dplyr::filter(cds.df,type == "CDS")$width),stringsAsFactors = F))
}))
> res.df
transcript_id start end
1 a 151 2092
2 b 92 782
Любой совет, как ускорить это? возможно использование dplyr