1

我有一个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 组成，这些 s 是沿seqnames（即染色体）的线性区间，由start和定义end。这些exons 具有取决于strand（对于每个 s 的所有exons都相同transcript_id）的方向/方向性：如果strand == "-"then 实际上end是外显子的第一个位置并且start是它的最后一个位置。如果strand == "+"thenstart和end分别是第一个和最后一个位置。

这些CDS线是exons 的子集。通常，在 each 中保留第一个和最后几个exons transcript_id，每个exon都有相同的CDS间隔（就坐标而言）。但是，例外情况是：

  1. CDS间隔是 s 的子集，exon这意味着它们可以在exon(s) 内开始和/或结束。可以是第一个CDS区间的第一个位置在最接近的那个之后，和/或最后一个区间exon的最后一个位置在最接近的那个之前。也有可能 a将有一个（因此也是一个）满足这些定义。CDSexontranscript_idexonCDS
  2. transcript_id's 其中所有exons 具有相同的CDS间隔
  3. transcript_id's 其中exon不具有相同CDS间隔的 s 只是前几个
  4. transcript_id's 其中exon不具有相同CDS间隔的 s 只是最后几个

我正在寻找一种快速方法function，它将返回CDS相对于 s 已组合的坐标transcript_id的exon局部坐标。本质上，start结果是直到第一个区间+第一个区间宽度data.frame的总和-它是（或如果）的子集，并且是宽度的总和。exonCDSstartCDSstartexonendstrand == "-"endCDS startCDS

到目前为止，这是我正在做的事情，但速度很慢：

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

4

0 回答 0