提取 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)
}
网友评论