使用不匹配的gtf版本信息会导致近一半ID无法转化

太多人有这样的疑问，为什么自己进行ID转换的时候，成功率很低，今天就为你解惑。

当我们遇到了这样一个表达矩阵（其实就是从tcga数据库下载，通过ucsc的xena浏览器啦）：

dim(exprSet)
exprSet=exprSet[apply(exprSet,1, function(x) sum(x>1) > 2),]
dim(exprSet)

> dim(exprSet) 
[1] 45821   122
> exprSet[1:4,1:4]

TCGA-A1-A0SP-01A TCGA-E2-A1LS-01A TCGA-BH-A0B9-01A TCGA-E2-A1LL-01A
ENSG00000000003.13             4283             6734             8676             2456
ENSG00000000005.5                 6                6               13               21
ENSG00000000419.11             2508             2255             6767             1713
ENSG00000000457.12              958             4363             3373              557


看到基因名字是 ENSG00000000003.13  这样的ID， 就需要转换：

常规的做法是：

library(stringr)
surGenes=str_split(rownames(exprSet),'[.]',simplify = T)[,1]
library(clusterProfiler)
library(org.Hs.eg.db)
df <- bitr(unique(surGenes), fromType = "ENSEMBL",
           toType = c( "ENTREZID" ,'GENENAME','SYMBOL'),
           OrgDb = org.Hs.eg.db)
head(df)


但使用不匹配的gencode的gtf版本信息会导致近一半ID无法转化 !

'select()' returned 1:many mapping between keys and columns
Warning message:
In bitr(unique(surGenes), fromType = "ENSEMBL", toType = c("ENTREZID",  :
  49.13% of input gene IDs are fail to map...
> head(df)


那我们就要想如何进行完美的基因信息匹配，首先需要去ucsc的xena浏览器里面下载到encode.v22.annotation.gene.probeMap 文件，这个是他们的表达量矩阵的基因的id的注释信息。

代码如下：

a=read.table('gencode.v22.annotation.gene.probeMap',header = T)
head(a)
ids=a[match(rownames(exprSet),a$id),1:2]
head(ids)
# 可以达到完美匹配
> head(ids)
                      id     gene
58775 ENSG00000000003.13   TSPAN6
58774  ENSG00000000005.5     TNMD
54795 ENSG00000000419.11     DPM1
3850  ENSG00000000457.12    SCYL3
3845  ENSG00000000460.15 C1orf112
910   ENSG00000000938.11      FGR


但这还远远不够我们惊奇的发现一个有一些基因的id是对应同一个基因的名字，如下所示:

tmp=table(ids$gene)
head(ids[ids$gene %in% names(tmp[tmp==2]),])
ids[ids$gene %in% head(names(tmp[tmp==2])),]


输出结果：

                      id     gene
29654 ENSG00000150076.21    CCDC7
55841 ENSG00000160200.16      CBS
20157  ENSG00000213204.7 C6orf165
29651 ENSG00000216937.10    CCDC7
7190   ENSG00000241962.8  C2orf15
30312  ENSG00000242338.5   BMS1P4
33473  ENSG00000248671.6  ALG1L9P
33474  ENSG00000254978.2  ALG1L9P
30309  ENSG00000271816.1   BMS1P4
20156  ENSG00000272514.4 C6orf165
7191   ENSG00000273045.4  C2orf15
55140  ENSG00000274276.3      CBS


对于这样的基因有多种处理方式， 我们这里选择删除表达量低的基因即可。

代码如下:

colnames(ids)=c('probe_id','symbol')  
ids=ids[ids$symbol != '',]
dat=exprSet
ids=ids[ids$probe_id %in%  rownames(dat),]
dat[1:4,1:4]   
dat=dat[ids$probe_id,]

ids$median=apply(dat,1,median) #ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids
ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果s
dat=dat[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名
dat[1:4,1:4]  #保留每个基因ID第一次出现的信息


最后就得到了一个纯粹的表达量矩阵，如下:

> dat[1:4,1:4]   
       TCGA-A1-A0SP-01A TCGA-E2-A1LS-01A TCGA-BH-A0B9-01A TCGA-E2-A1LL-01A
ZZZ3               4583             2047             7294             1984
ZZEF1              2144              851             4761             1698
ZYX               10216            20907            14932            15101
ZYG11B             2017              966             3507             1320


之后的操作你懂的。。。

尽情的玩耍吧！