多个探针对应同一个基因到底该如何取舍

前些天我发现了乳腺癌领域的PAM50算法原理探索,在:PAM50的概念及分子分型算法原理 ,其实并不难,然后我注意到他们在挑选50个基因的时候,提到了多个探针对应同一个基因到底该如何取舍。

原文是:For probesets that map to identical Entrez gene names, select the one with highest IQR (for Affy, select mean for Agilent),也就是四分位间距IQR,这个概念主要是在boxplot图表里面显示出来。当然了,不同芯片平台也是有一些细微的差别。

其实没有标准答案的问题

三五年前我的博客:多个探针对应一个基因,取平均值或者最大值 就讨论过这个问题,很多人参与留言:

一代Array探针可以这么做,RNA seq会出现一个gene symbol对应多个isform的数据,(有点类似array的这种情况吧。)我问过俩老师:

  • 一个md Anderson 的老师说他们用最长的CCDS的那个transcript作为这个基因的代表
  • 另一个ucla的老师说他们是将所有的isform表达量加起来作为这个基因的表达量。

因为芯片技术已经被时代抛弃,ngs技术本来就有读成的局限性,不管是谁再问我这样的问题,我都是回答,并没有标准答案。但是我们给出的代码是值得学习的:

我的代码的进化历史

具体详见;多个探针对应同一个基因取最大值的代码进化历史 ,首先是使用split结合 sapply,然后是使用by函数,最后是使用duplicated和order函数。

## 制作好 ids和exprSet,分别是探针注释信息和表达矩阵
identical(ids$probe_id,rownames(exprSet))
dat=exprSet
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第一次出现的信息
dim(dat)

比如,如果你下载CCLE数据库的一千多个细胞系的RNA-seq的counts矩阵,如下:

> a1=read.table('~/Downloads/CCLE_RNAseq_genes_counts_20180929.gct.gz',skip = 2,header = T)
> dim(a1)
[1] 56202 1021
> a1[1:4,1:4] 
 Name Description X22RV1_PROSTATE X2313287_STOMACH
1 ENSG00000223972.4 DDX11L1 12 8
2 ENSG00000227232.4 WASH7P 1340 821
3 ENSG00000243485.2 MIR1302-11 4 1
4 ENSG00000237613.2 FAM138A 6 3

如果你需要把它变成基因名字的表达矩阵,也会遇到一些基因名字重合的问题。

dat=a1[, 3:10] # 随便取几个细胞系,第1,2列是基因名字
rownames(dat)=a1$Name
ids=a1[,1:2] # 第1,2列是基因名字
head(ids)
colnames(ids)=c('probe_id','symbol')

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第一次出现的信息

这个代码非常好用,你一定要学习哦!

文末友情宣传

强烈建议你推荐我们生信技能树给身边的博士后以及年轻生物学PI,帮助他们多一点数据认知,让科研更上一个台阶:

Comments are closed.