上期我们已经从GTEx下载提取了Liver的119个样本的表达谱文件，本期转换miRNA ID并将TCGA的样本与GTEx的进行合并.
PS:文末有福利.
我们已经得到的GTEx的表达谱文件部分内容如下：
而我们之前得到的TCGA的表达谱格式是这样的：
为了使两个数据库的ID对应起来，我们得找名称对应关系才行．这个可以在人类基因命名委员会（HGNC）（http://www.genenames.org） 找到．这个库包含了人类基因所有的命名，别名和ID，可以说是天下命名，无出其右．
可以直接在下载页面 （http://www.genenames.org/cgi-bin/statistics） 底部下载所有的命名关系的一个大表，但这次我们只需要处理miRNA的ID转换，那只下载miRNA的就行：
下载txt格式的文件后打开是这样的：
上色的三列是本次ID转换我们所需要的，我们将根据上面这个表把GTEx的symbol和ensemble ID转为和GTEx中的alias，但我们仔细看HGNC中的ensembel_gene_id与GTEx里的有点不同,GTEx里的每个ID后面还有小数点,而HGNC的并没有,这是由于ensemble版本造成的,只需看小数点前面的ID就可以了,小数点前面ID的一致就是同一个miRNA.这个问题也会在TCGA基因表达谱ID转换时出现,同理只需管小数点前面的ID就可以了.
Ｒ实现转ID的过程比较方便，Ｒ代码如下：
# 转ID:
rm(list=ls())
library(data.table)
library(stringr)
# GTEx的liver表达谱路径：
GTEx_liver_path<-"./GTEx_Analysis_v6p_RNA_Liver2.tsv"
# 读入表达谱文件：
GTEx_reads<-fread(GTEx_liver_path,sep="\t",header =T)
# 重命名:
names(GTEx_reads)[1:2]<-c("ensembl_gene_id","symbol")
# 去掉ensemble小数点:
GTEx_reads[,ensembl_gene_id:=str_split_fixed(GTEx_reads$ensembl_gene_id,"[.]",2)[,1]]
# HGNC下载miRNA表达谱文件路径:
HGNC_mir_path<-"../HGNC/RNA_micro.txt"
HGNC_mir<-fread(HGNC_mir_path,sep="\t",header =T)
# 提取HGNC_mir中需要的3列:
HGNC_mir<-HGNC_mir[,c("ensembl_gene_id","symbol","alias_symbol")]上面的代码完成了读入文件和统一列名和去除GTExensemb_gene_id中的小数点.但在合并时我们有三种选择:
# 只匹配ensembl_gene_id进行合并:
GTEx_reads_merged_id<-merge(GTEx_reads,HGNC_mir,by=c("ensembl_gene_id"))
# 只匹配symbol进行合并:
GTEx_reads_merged_symbol<-merge(GTEx_reads,HGNC_mir,by=c("symbol"))
# 匹配ensembl_gene_id和symbol进行合并:
GTEx_reads_merged_both<-merge(GTEx_reads,HGNC_mir,by=c("ensembl_gene_id","symbol"))只匹配
ensembl_gene_id进行合并;只匹配
symbol进行合并;匹配
ensembl_gene_id和symbol进行合并;
三种方法都尝试了下发现差异很大:
# 提取三次合并结果:
a1<-GTEx_reads_merged_id[,alias_symbol]
# 1172条
a2<-GTEx_reads_merged_symbol[,alias_symbol]
# 1249条
a3<-GTEx_reads_merged_both[,alias_symbol]
# 1113条结果发现通过两个条件(ensemb_gene_id和symbol)合并得到的结果最少(1113条),只通过symbol得到的最多(1249条).但不管是1113还是1249条,都远少于miRNA总数1881,ID转换率低于66.4%.我搜了几条发现这应该是由于版本差异导致的.GTE应该用的是GRCh37.p13,而HGNC上面是最新版的GRCh38.p10.我们通过一个两个例子来说明.
同一个ensemble ID:
ENSG00000211563在不同版本gencode的结果:
| Source | Ensembl_gnee_id | Symbol |
|---|---|---|
| Gencode_v19 | ENSG00000211563.2 | MIR338 |
| Gencode_v23 | ENSG00000211563.3 | AC115099.1 |
| Gencode_v26 | ENSG00000211563.4 | MIR3065 |
| HGNC | ENSG00000211563 | MIR3065 |
| GTEx | ENSG00000211563.2 | MIR338 |
同一个
miRNAhsa-mir-338 在不同版本gencode的结果:
| Source | Symbol | Ensembl_gnee_id |
|---|---|---|
| Gencode_v19 | MIR338 | ENSG00000211563.2 |
| Gencode_v23 | None | None |
| Gencode_v26 | MIR338 | ENSG00000283604.1 |
| HGNC | MIR338 | ENSG00000283604 |
| GTEx | MIR338 | ENSG00000211563.2 |
可以看出同一个ensembl_gene_id在三个不同版本的Gencode所对应的Symbol中都是不一样的,而同一个Symblol在不同版本的Gencode中的Ensembl_gene_id也不同.
但还是可以发现:
HGNC与最新的Gencode保持一致,而GTEx与Gencode_v19版本保持一致,也就是说,ID转换率低是由于新旧版本不匹配导致的,如果要在新旧版本间选个标准,最好还是选择新版本,因为新旧版本的差异是由于新版本对就版的更正所导致的.也就是说新版本对”基因”的注释更加准确,但转换率会比较低.
刚才提到GTEx用的是Gencode_v19版本的,我自己提取了Gencode_v19里所有的miRNA(共3055条),然后以ensembl_gene_id与GTEx的相匹配,结果匹配到了2839条,转换率为92.9%.但我们都知道miRNA一共才不到2000条,这3055条是包含可信度较低的,transcript_status "NOVEL"以及level 3的miRNA.为了使得转换得到的miRNA更准确,假阳性更低.也为了避免同一个ensembl_gene_id匹配到两个不同的miRNA,或者同一个symbol匹配到两个不同的ensembl_gene_id,我们还是以最新版本也就是HGNC版为标准,并且以ensembl_gene_id和symbol两个都匹配为条件进行ID转换.最后对列名排序后输出文件:
# 获取列名:
a_names<-names(GTEx_reads_merged_both)
# 重新排序列名:
a_names_sorted<-a_names[c(1,2,length(a_names),3:(length(a_names)-1))]
setcolorder(GTEx_reads_merged_both,a_names_sorted)
# 输出文件:
fwrite(GTEx_reads_merged_both,"./GTEx_Analysis_v6p_RNA_Liver_miRNA_HGNCmerged.tsv",sep = "\t",col.names = T)上图第三列即是merge过后得到的我们一般常见的miRNA 名称.基因ID转换与此类似.以后做RNA-seq时我们再进行示范.
点击左下角”阅读原文”即可下载最新的SCI影响因子汇总表.
更多原创精彩内容敬请关注生信杂谈：