R语言基因ID转换

基因ID转换可分为在不同物种间的Gene ID转换和同一物种不同基因ID类型的转换。今天分享的主要内容是针对同一物种不同基因ID类型的转换。

我们常见的基因名称有：Gene ID(Entrez ID)、Gene Symbol和Ensemble ID等。

Gene ID(又称Entrez ID):是由 NCBI (National Center for Biotechnology Information) 提供的一种独特的标识符，用于标识基因、蛋白质、文献等生物学信息资源中的条目。Entrez ID 是一个数字编号，用于数据库中每个记录的唯一标识。其中一个 Entrez ID 可能对应有多个Gene symbol。

Gene Symbol：Gene Symbol(基因符号)是用来表示基因的标准化名称或缩写，它是生物学和医学研究中基因的常用标识符。基因符号通常由国际基因命名委员会(如 HGNC、MGI)制定，并且是基因的标准命名，用于避免混淆和提供一致性。Gene Symbol 有3大特点，包括标准化、广泛应用和与数据库相关联。

1.标准化：基因符号通常是简短的字母或数字组合，具有一致的格式，便于记忆和使用。

2.广泛应用：在基因研究、文献、数据库中广泛使用。基因符号通常是研究人员用来表示和识别基因的首选方法。

3.与数据库相关联：每个基因符号都可以与其他相关信息（如基因的功能、位置、蛋白质产品等）关联。

Ensembl ID：是 Ensembl 数据库中每个基因、转录本和蛋白质等生物学实体的唯一标识符。Ensembl ID 是一个标准化的标识符，用于在生物信息学研究中引用基因、转录本、外显子等生物学元素。常见的Ensembl ID包括Ensembl Gene ID、Ensembl Transcript ID和Ensembl Protein ID。

根据上述的解释说明，我们可以发现每种ID的应用场景都会有所不同。在大多数的情况下，我们在数据库下载的基因表达谱可能包含探针ID(可根据自身携带的平台文件进行匹配并进行去重)和Ensemble ID等。虽然Gene ID的表现形式和含义是丰富多彩的，但我们在进行一系列生信的分析的时候用到都是Gene Symbol。所以我们就需要进行不同ID之间的相互转换。

Count转换FPKM\RPKM、TPM、CPM

在生物信息学和数据分析中，count 数据类型通常指的是一个变量或数据集，其中包含对某些事件、物种、基因、类别等的计数。这种数据通常是整数型，表示发生次数或频率。

基因表达计数（RNA-Seq 数据）：在转录组学中，count通常指的是基因在样本中的表达量，以基因的转录本（如基因、外显子）为单位，表示每个基因在每个样本中的 reads数量（例如，RNA-Seq 数据中的 reads counts）。

在 RNA-Seq 中，count 数据通常指每个基因或转录本在样本中的读取数。常见的表达量计算方法之一是计数矩阵（Count Matrix），其中每一列代表一个样本，每一行代表一个基因，单元格中的数字即为该基因在该样本中的读取计数（reads count）。

Count 数据通常表示对某些事件、类别、物种或基因的计数，广泛应用于 RNA-Seq、微生物群落分析和统计学中。

这种数据类型通常是整数型，表示事件的发生次数或某种生物学单位的丰度。处理 count 数据时，通常需要进行标准化、转换以及差异分析等步骤。

处理 Count 数据的常用方法:

1.标准化：由于不同样本的测序深度可能不同，因此对 count 数据进行标准化是常见的步骤。例如，TPM（Transcripts Per Million）、FPKM（Fragments Per Kilobase of transcript per Million mapped reads）等方法可用于标准化 RNA-Seq 数据。

2.转换：在分析中，count 数据通常需要进行某种形式的转换。例如，使用 log2 转换，以减小高计数的影响或使数据符合正态分布。

3.统计模型：对于 RNA-Seq 和其他生物数据，常用的统计方法包括负二项回归（如 DESeq2 和 edgeR）来分析count 数据的差异。

上述我们分别阐述了Gene ID转换和Count转换成FPKM\RPKM、TPM、CPM。接下来我们用一个急性胰腺炎的实例进行相关分析的展示。

准备文件：

1. 1.  基因长度文件

2. 2.  急性胰腺炎的临床和表达数据

基因长度文件的获取(人)

首先，我们访问https://www.gencodegenes.org/human/，我们点击第一个的GTF按钮。之后解压在当前文件夹中就行。

下载急性胰腺炎的表达矩阵

点击http就可以下载了该表达矩阵。

接下来我们就使用代码演示，Gene ID转换和Count转换FPKM\RPKM、TPM、CPM

setwd("C:\\Users\\Administrator\\Desktop\\count_fpkm")

rm(list=ls())

#下载有基因长度的信息

library(GenomicFeatures) #加载R包

## 1.读取GTF文件

txdb <- makeTxDbFromGFF("gencode.v47.annotation.gtf") #TxDb：R语言中用于储存gtf文件的一种格式

## 2.提取外显子信息

exonic <- exonsBy(txdb, by="gene") #可以提取基因外显子部分，计算counts数据比对的为基因外显子序列

## 3.外显子长度求和

exonic.gene.sizes <- sum(width(reduce(exonic))) #reduce可以删除外显子重复序列部分

## 4.结果整理，获得ENSEMBL及长度

mygeneids <- data.frame(gene_id=names(exonic.gene.sizes),width=exonic.gene.sizes)

#ENSEMBL ID转换成Gene SYMBOL

library(org.Hs.eg.db) #人org.Hs.eg.db

library(AnnotationDbi)

mygeneids$gene_id <- gsub("\\..*", "",mygeneids$gene_id) #去除版本号

#匹配上一列gene SYMBOL

mygeneids$gene_symbol <- mapIds(org.Hs.eg.db,

keys=mygeneids$gene_id,

column="SYMBOL",

keytype="ENSEMBL",

multiVals="first") #名称转换由ENSEMBL转换为SYMBOL

rownames(mygeneids) <- c(1:nrow(mygeneids))

#导出包含ENSEMBL ID、Gene SYMBOL和width的文件

write.csv(mygeneids,"gene_id_width_gene_symbol.csv")

#导出包含Gene SYMBOL和width的文件

library(dplyr)

gene_length <- mygeneids %>%

distinct(gene_symbol,.keep_all = T) %>%

filter(!is.na(gene_symbol)) %>%

select(-gene_id)

gene_length <- cbind(gene_length$gene_symbol,gene_length$width)

colnames(gene_length) <- c("gene_symbol","width")

#基因-长度文件

gene_length <- data.frame(gene_length)

write.csv(gene_length,"gene_length.csv")

####GEO数据下载####

rm(list=ls())

#今天需要转换的示例文件是急性胰腺炎的数据

#临床数据的获取

library(GEOquery)

#有关R包的加载

library(affy)

library(limma)#处理均值

library(Biobase)

#数据下载

options(timeout=1000000)

gset <- getGEO('GSE194331', destdir=".",

AnnotGPL = T,     ## 注释文件

getGPL = T)       ## 平台文件

gset[[1]]

cli <- pData(gset[[1]])#临床信息

#导出临床数据

write.csv(cli,"GSE194331_cli.csv")

####ID转换####

#读取表达矩阵

rm(list=ls())

exp<-read.table("GSE194331_exp.txt", header = TRUE, sep = "\t")

gene_id_width_gene_symbol<-read.csv("gene_id_width_gene_symbol.csv")

gene_id_width_gene_symbol<-gene_id_width_gene_symbol[,-1]

exp_gene_id<-merge(exp,gene_id_width_gene_symbol,by.x="GeneName",by.y="gene_id")

exp_gene_id1<-exp_gene_id

exp_gene_id1 <-exp_gene_id1[, c(ncol(exp_gene_id1), 1:(ncol(exp_gene_id1)-1))]

exp_gene_id1<-exp_gene_id1[,-2]

#去重

matrix0 <- order(rowMeans(exp_gene_id1[,-1]),decreasing = T)

#调整表达谱的基因顺序

expr_ordered=exp_gene_id1[matrix0,]

#对于有重复的基因，保留第一次出现的那个，即行平均值大的那个

keep=!duplicated(expr_ordered$gene_symbol)#$Name是指基因名那一列

#得到最后处理之后的表达谱矩阵

expr_max=expr_ordered[keep,]

expr_max<-na.omit(expr_max)

rownames(expr_max)<-expr_max[,1]

expr_max<-expr_max[,-1]

exp_gene_id2<-expr_max

exp_gene_id2$width <- as.numeric(exp_gene_id2$width)

exp_gene_id2 <- as.matrix(exp_gene_id2) #转换为矩阵

counts<-exp_gene_id2

# counts to CPM

cpm <- apply(X =subset(counts, select = c(-width)),

MARGIN =2,

FUN =function(x){

x/sum(as.numeric(x)) * 10^6

})

# counts to RPKM\FPKM

geneLengths <- as.vector(subset(counts, select = c(width))) # 创建基因长度的向量，为下一步标准化做准备 

rpkm <- apply(X = subset(counts, select = c(-width)),

MARGIN = 2,

FUN = function(x) {

10^9 * x / geneLengths / sum(as.numeric(x))

}) # 计算rpkm

# counts to TPM

rpk <- apply( subset(counts, select = c(-width)), 2,

function(x) x/(geneLengths/1000)) #求基因长度归一化值

tpm <- apply(X = rpk,

MARGIN = 2,

FUN = function(x){

x / sum(as.numeric(x)) * 10^6

})#使用 rpk 值根据样本大小进行标准化

# FPKM\RPKM to TPM

fpkmToTpm <- function(fpkm) {

exp(log(fpkm) - log(sum(fpkm)) + log(1e6))

}## rpkm转换为tpm

#保存当前R语言环境中所有的变量

save.image("ID转换.RData")