FPKM值基因表达量的计算、基因ID转gene symbol的例子

高通量测序数据一般公司都会提供两种矩阵，一种是Row counts，前面说过的用于差异基因的筛选。第二种是FPKM值，可以理解为转录组基因的表达矩阵，可以用于做热图和基因表达变化的比较。但是数据挖掘中，我们只能下载到counts矩阵，所以要得到基因的表达量还需要通过计算。

这里我们使用biomaRt包举个例子，由counts值计算FPKM。

一、FPKM的计算

加载counts矩阵：

setwd("E:/生物信息学")ma <- read.csv("GSE169758_markdup.featco.2.counts.csv",header = T,row.names = 1)

下载安装R包：

BiocManager::install("biomaRt")library(biomaRt)

链接到ensembl数据库，并选择合适的物种，因为我们下载的是人的数据，所以dataset = "hsapiens_gene_ensembl"，小鼠的则选择dataset = “mmusculus_gene_ensembl”。计算基因长度。

ensembl <- useMart("ensembl") ensembl = useDataset("hsapiens_gene_ensembl",mart=ensembl)test <- getBM(attributes=c('ensembl_gene_id', 'start_position',                           'end_position','ensembl_transcript_id',                           'transcript_length'),mart = ensembl)head(test)#  ensembl_gene_id start_position end_position ensembl_transcript_id transcript_length#1 ENSG00000210049            577          647       ENST00000387314                71#2 ENSG00000211459            648         1601       ENST00000389680               954#3 ENSG00000210077           1602         1670       ENST00000387342                69#4 ENSG00000210082           1671         3229       ENST00000387347              1559#5 ENSG00000209082           3230         3304       ENST00000386347                75#6 ENSG00000198888           3307         4262       ENST00000361390               956

将得到的test文件排序，并去除重复的ID：

test <- test[order(test$ensembl_gene_id,test$transcript_length,decreasing = T),]g <- test[!duplicated(test$ensembl_gene_id),]g <- g[,c(1,5)]head(g)dim(g)summary(g)

取基因长度文件和count文件交集的基因：

ng=intersect(rownames(ma),g$ensembl_gene_id) length(ng)lengths=g[match(ng,g$ensembl_gene_id),2]names(lengths) <- g[match(ng,g$ensembl_gene_id),1]head(lengths)

根据最终选取的counts矩阵计算每个样本的文库大小：ma <- ma[names(lengths),]

total_count <- colSums(ma)head(total_count)#Mcc1     Mcc2     Mcc3     Mcc4     Mcc5     Mcc6 #39315944 15971423 16166354 25798072 33085950 1772396

计算FPKM，并保存文件：

ma_fpkm <- t(do.call( rbind,                      lapply(1:length(total_count),                             function(i){                               10^9*ma[,i]/lengths/total_count[i]                               #lengths向量自动遍历                             }) ))ma_fpkm[1:4,1:4]#[,1]      [,2]      [,3]        [,4]#ENSG00000000003  0.01340093  0.000000  0.000000  0.01021143#ENSG00000000005  0.00000000  0.000000  0.000000  0.00000000#ENSG00000000419 28.87258994 24.594527 23.186029 27.78157422#ENSG00000000457  2.23382616  1.677457  1.696455  2.31666061write.csv(ma_fpkm, file = "ma_fpkm.csv")

二、基因ID转gene symbol

从前面的结果我们很清楚的看到，每个基因都是用ID表示的，类似于ENSG00000000003，但看这个名称，根本不知道它是什么，我们通常也不习惯这样表示基因，还是习惯于Gene symbol。接下来我们将这个FPKM文件进行基因ID注释。

首先加载人基因数据库：

library(stringr)BiocManager::install("org.Hs.eg.db")library(org.Hs.eg.db)

准备三个文件，一个是gene ID，一个是gene symbol，还有一个是我们需要注释的文件，准备一列ensembl_id：

gs <- toTable(org.Hs.egSYMBOL)head(gs)#gene_id symbol#1       1   A1BG#2       2    A2M#3       3  A2MP1#4       9   NAT1#5      10   NAT2#6      11   NATPge <- toTable(org.Hs.egENSEMBL)head(ge)#gene_id      ensembl_id#1       1 ENSG00000121410#2       2 ENSG00000175899#3       3 ENSG00000256069#4       9 ENSG00000171428#5      10 ENSG00000156006#6      12 ENSG00000196136colnames(ma_fpkm) <- colnames(ma)ma_fpkm <- as.data.frame(ma_fpkm)ma_fpkm$ensembl_id <- row.names(ma_fpkm)

接着就是将这几个文件一一比对merge，就可以得到注释的gene symbol：

b <- merge(ma_fpkm,ge,by="ensembl_id",all.x=T)c <- merge(b,gs,by="gene_id",all.x=T) c=c[order(c$ensembl_id),]c=c[!duplicated(c$ensembl_id),]c=c[match(ma_fpkm$ensembl_id,c$ensembl_id),]head(c)#            Pan1      Pan2        Pan3        Pan4      Pan5        Pan6   symbol#29657  4.55462494  2.612894  3.31715074  3.80483778  4.487819  2.49863368   TSPAN6#27673  0.00000000  0.000000  0.00000000  0.02566601  0.000000  0.00000000     TNMD#33278 27.28248482 18.442794 22.18114506 23.02063714 25.941927 20.66173442     DPM1#26324  1.56335564  2.044088  2.33643810  2.39752223  1.953669  2.30114296    SCYL3#25673  4.02781699  4.294728  5.45416997  4.60184848  3.717569  5.34177442 C1orf112#15388  0.01591371  0.000000  0.01808745  0.01172831  0.000000  0.01215142      FGRwrite.csv(c,file = "genesymbol_fpkm.csv")

这样我们的问题就解决了！

接下来我们会着重于转录组数据的可视化和一些后续的分析！欢迎关注分享和交流讨论！

FPKM值基因表达量的计算、基因ID转gene symbol的例子相关推荐

matlab算幅值裕量,matlab计算相角裕度
问题: 如图所示的系统, 画出当K=45时的伯德图, 并确定增益裕度和相位裕度. 计算使系统稳定的最大K值, 并用劳斯阵列验证其结果. R + ? - K 1 (...... (0 1) 从而可 ...
uniport ID 转换为 gene symbol（ID转换）
网站 ID转换网址 input 选择 uniprot accession output 选择 Gene symbol 一般一次转换不能超过1万个基因ID,数据量太大就转换不出来,如果一次转换太多,需要 ...
Entrez ID 和 gene symbol 和 Ensembl ID转换
简单粗暴网页版 bioDBnet - Biological Database Network db2db模块可以非常方便地实现各种ID之间的转换,当然 bioDBnet的功能有很多,可以试试. 手动迁 ...
基因 ID 匹配利器
一.背景对于每个生物信息分析的人来说,ID 匹配(映射)是一项非常常见,但又很繁琐的任务.假设,我们有一个来自上游分析的 gene symbol 或报告的 ID 列表,然后我们的下一个分析却需要使用 ...
BiomaRt 包进行基因ID转换
R中两个常用的基因ID转换包:BiomaRt和clusterProlifer 1. clusterProlifer 总体来说clusterProlifer用起来是比较方便的,通过现成的org.*.e ...
R语言基因表达量转换（TPM、FPKM、RPKM）
基因表达量一般以TPM或FPKM为单位来展示. TPM,Transcripts Per Kilobase Million 计算公式: TPMi=(Ni/Li)*1000000/sum(Ni/Li+-- ...
2020.08.14【RNA-Seq流程】丨将HTseq生成的基因COUNT值转换为FPKM值
2021.04.22更新:2.0版本 count转换FPKM介绍公式介绍 R语言脚本公式介绍下方公式是FPKM的综述文章写的,参数并不好理解,但是对数量级表示很明白 https://harold ...
GTEx:基因型和基因表达量关联数据库
欢迎关注"生信修炼手册"! GTEx全称如下 Genotype-Tissue Expression 该项目对来自人体多个组合和器官的样本,同时进行了转录组测序和基因分型分析,构建了 ...
R语言survival包coxph函数构建cox回归模型、ggrisk包ggrisk函数可视化Cox回归的风险评分图、使用风险得分的中位数计算最佳截断值cutoff（基于LIRI基因数据集）
R语言survival包coxph函数构建cox回归模型.ggrisk包ggrisk函数可视化Cox回归的风险评分图.使用风险得分的中位数计算最佳截断值cutoff(基于LIRI基因数据集) 目录

FPKM值基因表达量的计算、基因ID转gene symbol的例子

FPKM值基因表达量的计算、基因ID转gene symbol的例子相关推荐

最新文章

热门文章