哭着求ID转换放过我--非模式生物转录组测序分析的辛酸史

上次我的推文里面说到，这次要学习ID转换，说到ID转换，虽然我只是一个新手，也没多少项目经验，但我可以说，我在ID转换上吃的苦，是成吨的。为了让大家意识到ID转换的基本情况和重要性，我们必须要讲一下转录组测序（RNA-seq）分析的那些事儿。如果不想看，可以直接跳到（三）ID转换

前几期的推文如下：

整天说要做数据挖掘，咱先把R安装上好么？保姆级R语言安装指导。

（一）没有人可以代替我去走弯路

今天可能会串线RNA-seq的linux上游分析。

想当初我故意说漏嘴，让老师知道我学了生信，然后把组里的数据给我处理，我就跌进了一个新手遇到就可能会弃坑的大坑，从此开始了hard模式。那时候曾老师说“怕什么，有我指导你，你还担心？”

讲真，现在我是真的怕了。

我的想法：拿到几个人类的数据，跑一跑曾老师上课时候教的流程，总之代码样样都是有的，只要重复就可以了。

实际情况：

工具从linux强行更换为OS：原来在曾老师那里学习的时候，已经有备好的服务器，所以从来没操心过工具，可是一换到自己的项目，突然发现没有服务器，只有一台iMAC，配置还不高，系统还不一样，曾经一度说要安装成Linux，各路大神惊讶得以为我被盗号了，其实我是真的不知道原来OS跟linux是类似的。
两系统代码有差异：因为工具从linux服务器换成OS系统，很多代码就开始不一样了，比如linux中用的wget，在OS中是curl -O。这构建环境变量，安装软件，每一步都走得磕磕绊绊。
数据类型愁死人：最大的难度，还是数据的限制。原来学习的时候，处理的都是挑好的数据，人类的样本+QC过关+有生物学重复，如果数据是这种类型，那直接复制粘贴代码运行即可。可是我的数据却是非重复样本的人类+非模式生物数据。就好像有了一种我换了泳衣戴了泳帽拿了游泳圈，然后你跟我说，今天我们爬山！

（二）非模式生物数据引发的血案

参考基因组+注释文件的下载：由于这部分之前并没有接触，都是老师已经安排好的，自己下载的时候是十脸懵逼，啊？什么是注释文件？怎么查到下载地址？ensemble，UCSC，NCBI到底挑哪个？

由于没有学好背景知识，也没见过猪跑，就只是下载参考基因组就折腾了很多次，加上非模式生物的参考基因组并不是每个平台都有的，以下就是我总结的经验：（A）平台推荐ensemble> NCBI> UCSC；（B）在哪里下载的参考基因组，就要在哪里下载注释文件，必须是同平台，不然容易串；（C）找下载链接可不是容易时，就算找到了该物种的数据，也不知道挑哪一个，不知道从哪里听来的toplevel，我就把toplevel的链接发给了大神，以至于他从1.0G的压缩包，直接解压出50G的内容。咳咳，想看我怎么坑他的，详情请看：

坑队友系列：ftp协议下载ensembl参考基因组，大概要点击原文才有得看了。

内容太多文字描述太费劲，因此我做了一个流程图（觉得PPT做得好看的请点赞or留言）。PPT制作插件大推荐~让PPT从此高大上起来！矬矬的屁屁踢千篇一律，神级的<em>ppt</em>毫不费力

（三）ID转换的多种方法

先看一下我们的数据，ensembl_id是从Linux的上游分析得到的，但是这些编号真是摸不着头脑，我们熟悉的只是基因名字，比如TP53，TGF-β1，EGFR等等，因此，为了方便大家沟通交流，我们还是要把平台id（例如，ensembl 平台的id，ensembl_id）转换我们常说的基因名，也就是gene symbol（大家统一使用的官方名字）

比如：朱一龙（官方名字，唯一的正统名，就像symbol）

粉丝叫他：亲儿子（那是粉丝平台自己给的名字，在粉丝内部流通，相当于ensembl id等各种平台自己给的命名，但真正拿出来沟通还是要说名字）

1. Method 1 直接通过gtf文件得到

在FeatureCounts计数的时候，改一个参数，将gene_id改为gene _name，可直接输出gene_name而不是gene_id。但是鉴于后续我们还需要使用ensembl id去做GO/KEGG等注释分析，所以ensembl id还是不可少的。但这个方法涉及到Linux操作，在这里我们就不讲了，感兴趣的小伙伴可以自己去查看一下FeatureCounts的帮助文档。

2. 使用R包org.Hs.eg.db

这个方法是曾老师教程里面提到的，这里我们用data_1作为例子

## 安装并加载R包
BiocManager::install('org.Hs.eg.db')
library(AnnotationDbi)
library(org.Hs.eg.db)

k<-keys(org.Hs.eg.db,keytype = 'ENSEMBL')  # 30292 elements
head(k,10)
?keys  # 注释数据库对象及其子代、方法等，请查看使用方法
########################## 以上内容若是看起来困难，那就不要理解了，先放着

## 提取出list，就是目前已有的ensembl_id和gene name的对应表，约有3万个基因
list<-select(org.Hs.eg.db,keys = k,
             columns = c('ENTREZID','SYMBOL'),
             keytype = 'ENSEMBL')
head(list,5)   ## 简单查看一下list的前五行内容，熟悉list内容（如下图）

## 把ensembl id变成列名
rownames(data_1)<- data_1$ensembl_id

## 筛选出我们数据在list中存在的基因
gene_id<- data_1[rownames(data_1) %in% list$ENSEMBL,]

## ID转换
colnames(gene_id) # 先把ensembl这一列的列名改成和list一样的名字
colnames(gene_id)<- c('ENSEMBL','control_1','control_2','control_3',
                      'experiment_1','experiment_2','experiment_3')
colnames(gene_id)
gene_id<- merge.data.frame(gene_id,list,by ='ENSEMBL' )

以上每一个代码都需要认真理解，因为将来会用到很多次

3. 使用biomart进行ID转换

上面那个R包是针对人类数据的ID转换的，一到非模式生物，就是永不了，后来在各位大佬的指点下，使用biomart完成了非模式生物的ID转换。这里我们用data_2进行模拟（注意data_1=data_2，这是我们上一篇推文中使用不同方法读取的同一个数据）

## 安装包并查看数据
BiocManager::install('biomaRt')
library(biomaRt)
## 尝试用biomart http://www.biotrainee.com/thread-1789-1-1.html 

listMarts()  ## 查看目前提供的数据库

##  biomart                  version
# 1 ENSEMBL_MART_ENSEMBL     Ensembl Genes 94
# 2 ENSEMBL_MART_MOUSE       Mouse strains 94
# 3 ENSEMBL_MART_SNP         Ensembl Variation 94
# 4 ENSEMBL_MART_FUNCGEN     Ensembl Regulation 94

##
my_mart<- useMart('ENSEMBL_MART_ENSEMBL')  # formal class mart

## 查看数据集
datasets<- listDatasets(my_mart)
datasets
dim(datasets)  ## 138行，3列，也就是说，有138个数据集

#######################以上这部分内容如果看不懂，也无所谓，可以不去理解，就这么跑下来也可以。

准备工作完成后，开始正式调试

## 查找目标基因组 确定物种正确http://asia.ensembl.org/Macaca_fascicularis/Info/Index 
# 人类的
human_dataset<- useDataset('hsapiens_gene_ensembl',mart = my_mart) # 1.1 mb

##
data_2_value<- data_2$ensembl_id  
attr1<- c('ensembl_gene_id','chromosome_name',
          'hgnc_symbol','hgnc_id')  ## 出来的表格会包含这4个列名

data_2_gene<- getBM(attributes = attr1,  ## 四个需要输出的列名
                    filters = 'ensembl_gene_id', ## 通过ensembl_gene_id 来筛选数据
                    values = data_2_value, ## 需要筛选的内容
                    mart = human_dataset)  ## 用于筛选的数据库

得到的数据如下：

整合数据

## 照例，还是把data_2中的ensembl_id这一列的名字，改得和需要合并的数据一致
colnames(data_2) # 先把ensembl这一列的列名改成和list一样的名字
colnames(data_2)<- c('ensembl_gene_id','control_1','control_2','control_3',
                     'experiment_1','experiment_2','experiment_3')
colnames(data_2)
data_2<- merge.data.frame(data_2,data_2_gene,by ='ensembl_gene_id' )

通过检查数据发现，使用biomart，所有的ID都转换成功了，而R包ENSG00000256904的ID转换失败了，其最主要的原因是list中没有这个ensembl_id的数据。所以，我们要尽量使用新版本的R包或者更加强大的数据库来进行ID转换。

好了，终于写完了。。我可以去睡午觉了。。。拜拜

点击原文，看我是怎么坑大神的。

本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。