getGeneLengthFromGTF.R

## Please see https://www.biostars.org/p/196548/#196757
## Transcript Length
library(GenomicRanges)
library(rtracklayer)

## reading your mm10 gtf file and read only exons
GTFfile <- gzfile(description = "E:/gencode.vM13.annotation.gtf.gz", open = "r") ## your mm10 gtf file
GTF <- import.gff(GTFfile, format="gtf", genome="mm10", feature.type="exon") 
grl <- reduce(split(GTF, elementMetadata(GTF)$gene_id))
reducedGTF <- unlist(grl, use.names=T)
elementMetadata(reducedGTF)$gene_id <- rep(names(grl), elementNROWS(grl))
elementMetadata(reducedGTF)$widths <- width(reducedGTF)

## calculation of length
calc_length <- function(x) {
  sum(elementMetadata(x)$widths)
}
output <- t(sapply(split(reducedGTF, elementMetadata(reducedGTF)$gene_id), calc_length))
output <- data.frame(t(output))
output$id = rownames(output)
output$id = gsub(pattern = "(.*)[.](.*)", replacement = "\\1", x = output$id)
output = output[,c(2,1)]
colnames(output)[2] = "Transcript.Length"
head(output)

write.csv(output, file = "gencode.vM13.geneLenght.csv")

 

转载于:https://www.cnblogs.com/qiniqnyang/p/7518991.html