【已知基因名得到基因长度】计算FPKM、TPM用

很多数据从网上下过来只有gene symbol和counts数，想要计算FPKM就要得到基因长度。

得到基因长度可以这样去做，首先去ensembl下载对应物种的gtf（不知道如何下，可以自行搜索）。

##导入数据与计算
library(GenomicFeatures) #导入包，没有安装这个包可以搜一下如何安装
txdb <- makeTxDbFromGFF("Mus_musculus.GRCm39.105.gtf",format="gtf")##这里以小鼠的gtf文件为例，将数据导入为TxDb对象
exons_gene <- exonsBy(txdb, by = "gene") 
exons_gene_lens <- lapply(exons_gene,function(x){sum(width(reduce(x)))}) #计算总外显子长度
exons_gene_lens[1:10]

输出的exon_gene_lens如下图

得到的基因长度其实个列表，如上图的列表为例。第一个元素是3262，元素名是ENSMUSG00000000001。利用下面的代码可以转化为dataframe。

##转换为dataframe
gene_length <- sapply(exons_gene_lens,function(x){x})
id_length <- as.data.frame(gene_length)
id_length

id_length结果如下图,行名为geneid，第一列为gene length，这里的gene length是全外显子长度总和。

 接下来在对应上gene symbol就行了。

上ensembl BioMart下载geneid和genename对应的文件。

 

选择以TSV的格式导出 ，得到如下文件 

 将这个文件导入到R里，这里注意要去掉列名里的Gene stable ID和Gene name的空格，然后全选复制。

GN <- read.table("clipboard",header = T,fill = T) #有些geneid没有对应的gene name所以要写参数fill=T
GN #输出确认下

GN是这样的

 再合并一下

rownames(GN) <- GN$GenestableID #用geneid作为行名
GN<- GN[sort(rownames(GN)),] #根据行名进行排序
id_length <- id_length[sort(rownames(id_length)),] ##根据行名进行排序；因为只有一例，这里会变成向量，但没关系
id_GN_length <- cbind(GN,id_length) #合并一下
colnames(id_GN_length) <- c("gene_id","gene_name","length") #改一下行名
id_GN_length #得到最终的dataframe

 最终得到的dataframe是这样的

后面和自己的data合并再一起，就不写了，就是一些涉及数据框的基础操作了。

本文部分涉及GenomicFeatures的代码参考自http://www.bioinfo-scrounger.com