科研信息学

如何提取gtf文件中transcript的长度

2020-12-12  本文已影响0人  落寞的橙子

提取 genecode的gtf注释信息
函数如下,调用方式

cout_transcript_length("/you_gtf_dir/your.gtf")
cout_transcript_length<-function(gtf_dir){
library(GenomicRanges)
library(rtracklayer)
suppressMessages(library(tidyverse))
  gtf <- rtracklayer::import(gtf_dir)
  gtf_df=as.data.frame(gtf) 
  gtf_df_exon<-subset(gtf_df,type=="exon")
  transcripts<-as.character(unique(as.character(gtf_df_exon$transcript_id)))
  
  out_tab<-data.frame()
  for (i in transcripts){
    transcripts=i
    rt<-gtf_df_exon %>% filter(transcript_id==i)
    if (nrow(rt)>1){
      lengths=sum(rt$width)
    }else{
      lengths=as.numeric(rt$width)
    }
    out_tab=rbind(out_tab,cbind(transcripts,lengths))
  }
  return(out_tab)
}
上一篇 下一篇

猜你喜欢

热点阅读