Нахождение локальных координат в линейных интервалах

У меня есть 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 состоит из набора exons, которые являются линейными интервалами вдоль seqnames (которые являются хромосомами), определяемые как start а также end, Эти exonимеют ориентацию / направленность, которая зависит от strand (что одинаково для всех exonс каждого transcript_id): если strand == "-" тогда на самом деле end это первая позиция экзона и start это его последняя позиция. Если strand == "+" затем start а также end являются первой и последней позициями соответственно.

CDS линии являются подмножеством exons. Обычно оставляют первые и последние несколько exonв каждом transcript_idкаждый exon имеет идентичный CDS интервал (с точки зрения координат). Тем не менее, исключения:

  1. CDS интервалы являются подмножеством exons означает, что они могут начинаться и / или заканчиваться в течение exon(Ы). Вполне возможно, что первая позиция первой CDS интервал после того из ближайшего exonи / или последняя позиция последней CDS интервал до того из ближайшего exon, Также возможно, что transcript_id будет иметь один exon (и, следовательно, один CDSотвечая этим определениям.
  2. transcript_idгде все exonс идентичными CDS интервалы
  3. transcript_idгде exonы, которые не имеют идентичные CDS интервалы только первые
  4. transcript_idгде exonы, которые не имеют идентичные CDS интервалы только последние

Я ищу быстро function который вернет локальные координаты CDS относительно координат transcript_id чья exonс были объединены. По существу, start в результате data.frame это сумма exon ширина интервала до первого CDS интервал + start из первых CDS - start из exon это подмножество (или end если strand == "-") а также end является CDSstart + сумма 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

0 ответов

Другие вопросы по тегам