在R中使用refGenome阅读GTF格式的文件

install.packages("refGenome")
library(refGenome)

create ensemblGenome object for storing Ensembl genomic annotation data

ens <- ensemblGenome()

read GTF file into ensemblGenome object

read.gtf(ens, "Homo_sapiens.GRCh38.90.chr.gtf")

read.gtf(ens, "Homo_sapiens.GRCh38.90.chr.gtf")
[read.gtf.refGenome] Reading file 'Homo_sapiens.GRCh38.90.chr.gtf'.
[GTF] 2612134 lines processed.
[read.gtf.refGenome] Extracting genes table.
[read.gtf.refGenome] Found 58,243 gene records.
[read.gtf.refGenome] Finished.

看一下他属性

class(ens)
[1] "ensemblGenome"
attr(,"package")
[1] "refGenome"

本次的重点来了create table of genes

my_gene <- getGenePositions(ens)

class(my_gene)
[1] "data.frame"
dim(my_gene)
[1] 58243 26

这26个是变量分别是：

testmygene <- my_gene[1:5,]
View(testmygene)
names(testmygene)
[1] "id" "seqid"
[3] "source" "start"
[5] "end" "score"
[7] "strand" "frame"
[9] "ccds_id" "protein_version"
[11] "protein_id" "exon_version"
[13] "transcript_support_level" "tag"
[15] "exon_number" "gene_id"
[17] "transcript_name" "gene_version"
[19] "gene_name" "gene_source"
[21] "gene_biotype" "transcript_id"
[23] "transcript_source" "exon_id"
[25] "transcript_version" "transcript_biotype"

my_gene %>% group_by(seqid, gene_biotype) %>% summarise(count = n()) -> my_tally
ggplot(my_tally, aes(x =seqid, y = log2(count))) +
  geom_bar(aes(fill = gene_biotype), stat = 'identity', position = 'dodge')