#创作灵感#
记录工作实践、项目复盘
写技术笔记巩固知识要点
作为一个R的小白，记录在用metascape分析得到的数据做GO Chord (GO的弦图) 时遇到的困难，提供给大家，愿少走弯路。
背景：我将DEG（差异基因）的列表上传到metascape。下图1是得到的分析结果。
图1.metascape原始分析结果
我从中选择了8个代表性的GOterm(综合log p value，FDR（logq value） 和 enrichment score)。如下图2。
图2.从metascape原始分析结果挑出8个代表性的GOterm
利用上图2，表中的C列和P列的数据，进行数据转换，形成了GO Chord (GO的弦图)的作图文件。如下图3.
图3. GO chord做图文件 无logFC版
如图1所示，大家可以看到metascape分析得到的GO term enrichment 数据中不包含基因的log2FC。用图3这种不包含基因log2FC的数据也可以绘制GO Chord (GO的弦图)，但是基因都是灰色色块表示，不够美观。如下图。
图4.没有log2FC的GO chord图。
做GO Chord (GO的弦图)提供这些基因的log2FC会更好看。
最简单的方法，就是用EXCEL处理。从DEG 文件（测序公司给出的原始的含有所有DEG的log2FC值的文件），如下图5。
图5. DEG文件，显示所有原始数据。
复制图5所示表中的C列external_gene_name和K列log2FC两列的数据，粘贴到图3GO数据表中KL两列，在表中I列中输入VlookUP函数，公式为=VLOOKUP(A2,K:L,2,FALSE)。这样就找到了目的基因的log2FC。最后把不要的数据删掉。如下图6所示。
图6. 用EXCEL找到log2FC。
但由于我要处理的DEG表非常之间的多，用这个方法就很耗时间（蠢）。（然后我就选了看似简单，实际上更耗时间（蠢）的方法）。我知道对大神来说肯定有更简单的R语言，可惜我现在还没学会，只能用笨办法。
最初，我的思路是利用两个表中的相同列external_gene_name比较，用merge函数将图5.DEG 文件中目的基因的log2FC合并到图3.GO Chord作图文件里面进去。
下面是我merge的程序
然而，问题出现了！ 居然找不到所有的拟作图基因(为了美观，我也就选了30个左右的基因来做GO Chord)！很多基因的log2FC返回的是空值，如下图7。
图7，根据Gene_name 无法找到所有的基因，例如kif25基因。
明明利用EXCEL查找功能，这些基因可以在DEG数据中找到！！！！如下图8中找到的kif25基因。
图8，DEG文件中有kif25基因。
相信细心的你已经发现为啥出错了。
对，R语言它严格区分大小写！
我DEG文件中的Gene_name是大写的KIF25！而我Metascape文件中的Gene_name它是小写的kif25！
根据Gene_name来merge的计划是失败了。于是，我转变思路，想将 图3.GO Chord作图文件里的那几十个基因转换成ensembl_gene_ID，然后根据ensembl_gene_ID来merge。
下面是我的程序
rm(list = ls())#清空环境用
library('biomaRt')#没有安装的要先安装
library(tidyverse)#没有安装的要先安装
listMarts()#查看下一步要安装的数据库，如果知道要安装数据库名，如“ensembl”，此步骤可删除
my_mart <-useMart("ensembl")#选择Ensembl数据库
datasets=listDatasets(my_mart)#查看emsebl中不同物种77个数据集，确定下一步要载入的数据集名称，如果知道要安装数据库名，如“drerio_gene_ensembl”，此步骤可删除
view(datasets)#查看emsebl中不同物种77个数据集，确定下一步要载入的数据集名称，如果知道要安装数据库名，如“drerio_gene_ensembl”，此步骤可删除
my_dataset=useDataset("drerio_gene_ensembl",mart = my_mart)#我的是斑马鱼数据，所以用了“drerio_gene_ensembl”数据集
# Unkown ID 转 ensembl测试
Testdata2<- data.frame(UnkownID=c("cep20","bub1ba","danh1" , "eme1",
"hipk1a","kif25","lca5","pfkfb2b","slc2a11a"))#输入数据
view(Testdata2)#检查输入结果用，可删除
getBM(attributes = c("ensembl_gene_id","external_gene_name"),
filters = "external_gene_name",#此处是问题的关键所在，请看后面介绍
values = Testdata2$UnkownID,
mart = my_dataset)
在getBM函数选择filter参数的时候，我选择了最常用的"external_gene_name"，结果是------还是失败！8个问题基因只能转换出两个基因的名字，如下图9。
图9.利用Metascape数据中的基因名字无法转换成ensembl_gene_ID。
到此为止，我才弄明白，我DEG文件中的Gene_name根据的是ensembl的“external_gene_name”，所以是大写的KIF25！而我Metascape文件中的Gene_name，它是小写的kif25，它根据的不是ensembl的“external_gene_name”！
也就是说，它两不是一个妈。之前不能merge的根本原因就在这！
于是，问题来了！！！！！！！！请问Metascape文件中的Gene_name,它根据的到底是是什么类型的基因 ID ？用的哪个的命名？？？？？？？？
我在网上搜了半天，发现有人说Metascape 可以识别的基因 ID 包括「Gene Symbol」、「RefSeq」、「Entrez Gene ID」。
但是就是没人说它输出的是什么类型的基因 ID 。
于是我只能用了一个笨法子。用下面的程序查看ensembl数据库中的filters，如图10。然后一个一个拿来做filter参数尝试。
#查看可供选择的filters
filters=listFilters(my_dataset)
view(filters)
图10.所有的filters
如果你也想看看attributes里面可以输出有多少种类型的数据，你可以用下面程序，如图11。
#查看可供选择的attributes
attributes=listAttributes(my_dataset)
view(attributes)
图11.所有的attributes
经过我的不懈努力和反复尝试，终于成功了，如图12。
答案揭晓了：Metascape文件中的Gene_name,它根据的是-------wikigene_name!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
图12.Metascape数据中的问题基因最终能够转换成ensembl_gene_ID。
后面我的程序就简单了，先将所有类似图3中的“选择8个GOTERM数据.xlsx”文件另存为“Metascape GO result.csv”。保存在桌面“BiomaRt GeneID conversion”文件夹里。再进行转换得到ensembl_gene_id.
下面是我最终的程序
rm(list = ls())#清空环境用
library('biomaRt')#没有安装的要先安装
library(tidyverse)#没有安装的要先安装
listMarts()#查看下一步要安装的数据库，如果知道要安装数据库名，如“ensembl”，此步骤可删除
my_mart <-useMart("ensembl")#选择Ensembl数据库
datasets=listDatasets(my_mart)#查看emsebl中不同物种77个数据集，确定下一步要载入的数据集名称,如果知道要安装数据库名，如“drerio_gene_ensembl”，此步骤可删除
view(datasets)#查看emsebl中不同物种77个数据集，确定下一步要载入的数据集名称，如果知道要安装数据库名，如“drerio_gene_ensembl”，此步骤可删除
my_dataset=useDataset("drerio_gene_ensembl",mart = my_mart)#我的是斑马鱼数据，所以用了“drerio_gene_ensembl”数据集，可根据你的需要修改参数
setwd("C:\\Users\\Admin\\Desktop\\BiomaRt GeneID conversion")#设置读取的文件路径，可根据你的需要修改参数
Metascape GO <- read_csv("Metascape GO result.csv") #读取文件，可根据你的需要修改参数
symbols<- getBM(
attributes = c("ensembl_gene_id","wikigene_name"),#你希望返回的结果，可根据你的需要修改参数，但必须至少有这两项
filters = "wikigene_name",#此处为重点，请你搞清楚自己的基因ID是什么类型的
values = Metascape GO$Gene_ID,#要转换的基因名字数据列，你的文件中有一列名字必须是Gene_ID，可根据你的需要修改参数
mart = my_dataset)
symbols <- symbols[which(symbols$ensembl_gene_id !=""),]#去掉没有找到基因ensembl id的
Metascape GO2 <- right_join(symbols,Metascape GO,by=c("wikigene_name"="Gene_ID"))#合并文件
view(Metascape GO2)#检查结果用，可删除
write.csv(Metascape GO2,file ="Metascape GO result new.csv" )#生成新文件
如果你的数据和我相反，需要将ensembl _gene_id转gene name，或者添加其他信息，可以使用下面程序。
rm(list = ls())#清空环境用
library('biomaRt')#没有安装的要先安装
library(tidyverse)#没有安装的要先安装
listMarts()#查看下一步要安装的数据库，如果知道要安装数据库名，如“ensembl”，此步骤可删除
my_mart <-useMart("ensembl")#选择Ensembl数据库
datasets=listDatasets(my_mart)#查看emsebl中不同物种77个数据集，确定下一步要载入的数据集名称,如果知道要安装数据库名，如“drerio_gene_ensembl”，此步骤可删除
view(datasets)#查看emsebl中不同物种77个数据集，确定下一步要载入的数据集名称，如果知道要安装数据库名，如“drerio_gene_ensembl”，此步骤可删除
my_dataset=useDataset("drerio_gene_ensembl",mart = my_mart)#我的是斑马鱼数据，所以用了“drerio_gene_ensembl”数据集，可根据你的需要修改参数
#ensembl gene id 转gene name 测试
Testdata<- data.frame(ensembl_gene_id=c("ENSDARG00000100115",
"ENSDARG00000017532",
"ENSDARG00000076153",
"ENSDARG00000002037",
"ENSDARG00000071198"))#输入几个你的ensembl gene id
view(Testdata)#检查用，可删除
getBM(attributes = c("ensembl_gene_id","external_gene_name","wikigene_name"),
filters = "ensembl_gene_id",
values = Testdata$ensembl_gene_id,
mart = my_dataset)
# 将文件中 ensembl gene ID转基因name
setwd("C:\\Users\\Admin\\Desktop\\BiomaRt GeneID conversion")#你的文件存的路径，可根据你的文件修改参数
myresult <- read_csv("Myresult.csv") #你保存为csv格式的Myresult文件，可以用EXCEL另存为该格式，可根据你的文件修改参数
symbols<- getBM(
attributes= c("ensembl_gene_id","wikigene_name"),#可添加多种其他类型信息，可根据你的需要修改参数，但必须至少有两项
filters = "ensembl_gene_id",#你想转换的数据类型，可根据你的文件修改参数
values = myresult$Gene_ID,#你想转换所依据的数据，在这里该列内容应该是ensembl _gene_id，你的文件中必须含有Gene_ID这一列名，可根据你的文件具体情况修改本参数
mart = my_dataset)
symbols <- symbols[which(symbols$external_gene_name !=""),]#去掉没有找到基因name的
myresult2 <- right_join(symbols,RNA_seq,by=c("ensembl_gene_id"="Gene_ID"))#合并文件
view(myresult2)#检查结果用，可删除
write.csv(RNA_seq2,file ="Myresult new.csv" )#生成新文件
522
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
