github-actions[bot]
commented
3 years ago 使用不匹配的gtf版本信息会导致近一半ID无法转化 by 生信技能树
太多人有这样的疑问，为什么自己进行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
之后的操作你懂的。。。
尽情的玩耍吧！
大概率就是差异分析，基本上看我五年前的《数据挖掘》系列推文 就足够了；
https://mp.weixin.qq.com/s/LpZ-wJmzR-z0rnuOPZwlFA