在RNA-seq或芯片数据下游分析中经常遇到需要将基因表达矩阵行名的ensembl_id ( gene_id ) 转换为gene symbol(gene_name)的情况，而在转换时经常会出现多个ensembl_id对应与一个gene symbol的情形，此时就出现了重复的gene symbol。
重复的gene symbol当然是不能作为基因表达矩阵行名的，此时就需要我们去除重复的gene symbol。
gene symbol去重复有一般有两种思路：
1.只保留平均表达量最高的gene symbol
2.合并所有gene symbol（基因表达量进行加和或者取平均）
一、获取ensembl_id与gene symbol的对应文件
- 首先需要得到所需的gtf文件（最好是上游基因计数时所用文件），一般在gencode下载GENCODE - The Mouse GENCODE Release History，本次示例选取小鼠mm10（grcm38）基因组版本，因此选取GENCODE 对应版本M25，选择regions为ALL的GTF文件进行下载即可
image.png
image.png
- 接着需要提取gtf文件中ensembl_id(gene_id)和gene symbol(gene_name)的对应关系，此步在linux或者R中操作都可以，我个人比较喜欢用linux命令，因此示范一下在linux中的操作，最后会得到g2s_vm25_gencode.txt文件
gunzip gencode.vM25.chr_patch_hapl_scaff.annotation.gtf.gz
vim gtf_geneid2symbol_gencode.sh
##################以下为.sh文件内容
gtf="gencode.vM25.chr_patch_hapl_scaff.annotation.gtf"
### gene_id to gene_name
grep 'gene_id' $gtf | awk -F 'gene_id \"' '{print $2}' |awk -F '\"' '{print $1}' >gene_id_tmp
grep 'gene_id' $gtf | awk -F 'gene_name \"' '{print $2}' |awk -F '\"' '{print $1}' >gene_name_tmp
paste gene_id_tmp gene_name_tmp >last_tmp
uniq last_tmp >g2s_vm25_gencode.txt
rm *_tmp
bash gtf_geneid2symbol_gencode.sh
image.png
二、取最高表达量的一个重复gene symbol
以下代码参考修改自ensembl_id转换之两种方法比较 - 简书 (jianshu.com)
整体思路：
构建包含ensembl_id、gene symbol和基因表达中位数的ids对象——将gene symbol按照基因表达从大到小排列——去重复gene symbol行——根据ids的行名保留表达矩阵并将行名转为gene symbol
###环境设置
library(tidyverse) # ggplot2 stringer dplyr tidyr readr purrr tibble forcats
library(data.table) #多核读取文件
head(counts) #counts是需要转换ensembl_id的表达矩阵
##从gtf文件提取信息,获得gencode的基因id对应symbol的ids矩阵
ids <- data.frame(geneid=rownames(counts),
median=apply(counts,1,median)) #计算基因表达中位数，用于之后排序
g2s <- fread('g2s_vm25_gencode.txt',header = F,data.table = F) #载入从gencode的gtf文件中提取的信息文件
colnames(g2s) <- c("geneid","symbol")
table(ids$geneid %in% g2s$geneid) #查看需要转化的geneid在g2s的匹配情况
ids <- ids[ids$geneid %in% g2s$geneid,] #取出在gencode数据库的gtf注释中能找到的geneid
ids$symbol <- g2s[match(ids$geneid,g2s$geneid),2] #match返回其第二个参数中第一个参数匹配的位置,把g2s的geneid按照ids$geneid的顺序一个个取出来，从而得到ids$symbol这一列
ids <- ids[order(ids$symbol,ids$median,decreasing = T),] #把ids$symbol按照ids$median由大到小排序
##去重复
dim(ids); table(duplicated(ids$symbol)) #统计查看重复的symbol
#[1] 56262 3
#FALSE TRUE
#55492 770
ids <- ids[!duplicated(ids$symbol),]#取出不重复的ids$symbol
dim(ids); table(duplicated(ids$symbol))
#[1] 55492 3
#FALSE
#55492
##转化geneid为symbol
counts <- counts[rownames(ids),] #取出表达矩阵中ids有的行
rownames(counts) <- ids[match(rownames(counts),ids$geneid),"symbol"] #根据geneid和symbol进行匹配
head(counts)
三、合并所有重复的gene symbol(推荐)
主要思路为利用aggregate函数根据symbol列中的相同基因合并基因表达矩阵矩阵
library(tidyverse) # ggplot2 stringer dplyr tidyr readr purrr tibble forcats
library(data.table) #多核读取文件
g2s <- fread('g2s_vm25_gencode.txt',header = F,data.table = F) #载入从gencode的gtf文件中提取的信息文件
colnames(g2s) <- c("geneid","symbol")
symbol <- g2s[match(rownames(counts),g2s$geneid),"symbol"] #匹配counts行名对应的symbol
table(duplicated(symbol))
#FALSE TRUE
#55492 770
##使用aggregate根据symbol列中的相同基因进行合并
counts <- aggregate(counts, by=list(symbol), FUN=sum)
counts <- column_to_rownames(counts,'Group.1')
dim(counts)
#[1] 55492 4
id转换前
id转换后
基因名去重复的一点思考：
这两种思路的差别在于,第一种只取表达量最高的基因，认为只有这个基因有意义，其余表达量靠后的相同基因不重要。第二种则是合并所有具有相同基因名的基因，考虑到了所有基因的表达情况，考虑更全面，因此个人更推荐。
实际分析中，由于我们一般差异分析是对不同样品中的同一个基因进行的比较，因此这两种方法实际差别并不大，按自己需求选择即可。
（其实分析时有些人嫌麻烦，还会直接合并symbol和表达矩阵后随缘保留一个重复基因，这种方法这就见仁见智了...）
网友评论