计算单个基因上exon的有效长度
2019-03-06 本文已影响1人
热衷组培的二货潜
哇,是时候要整理一下以前的学习记录了,以前傻里傻气不知道通过平台来记录代码,都是放在一个R脚本
里面结果导致今天别人问一个同样的问题,居然花了好久才找到自己以前的记录,赶紧记下来进行分类。
gff3文件格式参考
library(rtracklayer)
anno <- import.gff3("E://B512/Project/data/all.gff3")
exons <- anno[anno@elementMetadata$type %in% c("exon"),]
tmp <- split(exons,as.character(exons$Parent))
Gene_length <- sum(width(reduce(tmp)))
as.data.frame(Gene_length) -> PerTranscript_exon_length
write.table(PerTranscript_exon_length,"PerTranscript_exon_length.txt",quote =F)
小知识:也可以通过
featureCounts
指定exon
然后计算基因Counts
,也可以得到基因或者转录本的有效的exon
长度。