在RNA-seq下游分析中经常遇到需要将基因表达矩阵行名的ensembl_id ( gene_id ) 转换为gene symbol ( gene_name )的情况,而在转换时经常会出现多个ensembl_id对应与一个gene symbol的情形,此时就出现了重复的gene symbol。重复的gene symbol当然是不能作为基因表达矩阵行名的,此时就需要我们去除重复的gene symbol。
这里博主使用R语言,在表达谱数据中重复基因--取平均/取最大值
基因名去重复的一点思考:
这两种思路的差别在于,第一种只取表达量最高的基因,认为只有这个基因有意义,其余表达量靠后的相同基因不重要。第二种则是合并3所有具有相同基因名的基因,考虑到了所有基因的表达情况。实际分析中,由于我们一般差异分析是对不同样品中的同一个基因进行的比较,因此这两种方法实际差别并不大,按自己需求选择即可。
以博主在GEO数据库下载并合并好但并没有去除基因名有重复的表达谱数据为例。
rm(list = ls())
setwd("D:/R.result/2.He/ww2023.8.22_injury")#设置工作路径
library(limma)
library(pheatmap)
inputFile="Exp.txt"
rt=read.table(inputFile, header=T, sep="\t", check.names=F)#读取数据
这是GEO直接下载并合并的数据,不清楚是否有重复值,但无论是高通量数据和芯片数据,基本都会无可避免,出现重复基因名和对应的表达矩阵,所以为确保分析结果准确性,分析前应对数据做检测,去除重复基因以及对应表达矩阵。
第一种方法,博主觉得比较靠谱,对于相同的基因,我们应该挑选行平均值大的那一整行,不能打乱随便取。
#计算行平均值,按降序排列
matrix0 <- order(rowMeans(rt[,-1]),decreasing = T)
#调整表达谱的基因顺序
expr_ordered=rt[matrix0,]
#对于有重复的基因,保留第一次出现的那个,即行平均值大的那个
keep=!duplicated(expr_ordered$Name)#$Name是指基因名那一列
#得到最后处理之后的表达谱矩阵
expr_max=expr_ordered[keep,]
接下来就是可以继续分析啦,当然博主建议先把这个矩阵保留输出,再进行分析...继续上代码~
expr_max=as.matrix(expr_max)#对上述矩阵转化
rownames(expr_max) <- expr_max[,1]
exp=expr_max[,2:ncol(expr_max)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
lea=rbind(id=colnames(data),data)
write.table(lea,file = "GSE.txt",sep = "\t",quote = F,col.names = F)#输出名为GSE.txt的文件
这个就是对于相同的基因,挑选行平均值大的那一整行,不打乱随便取。
第二种,取平均值,博主觉得这个可能有点不太行,但是结果和第一种差不多,大家有兴趣可以试试,话不多说,上代码~
rm(list = ls())
setwd("D:/R.result/2.He/ww2023.8.22_injury")#设置工作路径
library(limma)
library(pheatmap)
inputFile="Exp.txt"
rt=read.table(inputFile, header=T, sep="\t", check.names=F)#读取数据
#1)利用aggregate函数,对相同的基因名按列取平均
expr_mean=aggregate(.~Name,mean,data=rt)
这就是一种取平均值的方法
2)对于重复的基因名字,取表达值最大的哪一行
其实aggregate也可以对相同的基因使用max函数取最大值,但是这样处理是有问题的。例如同一个基因出现了三次,那么会有三行数据。如果使用aggregate+max,对于每一个样本,他会从三个值中挑选最大的那个值最为这个样本的表达值,这样做是不科学的。
#利用aggregate函数,对相同的基因名按列取取最大值
expr_max=aggregate(.~Name,max,data=rt)
这是第二种取最大值方法
后续保存对应矩阵的话,按照第一种方法后续操作即可~
今天就到这里啦~