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

前些天我发现了乳腺癌领域的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，帮助他们多一点数据认知，让科研更上一个台阶：

生信爆款入门-全球听（买一得五）（第4期），你的生物信息学入门课
数据挖掘第2期（两天变三周，实力加量），医学生/临床医师首选技能提高课
生信技能树的2019年终总结，你的生物信息学成长宝藏
2020学习主旋律，B站74小时免费教学视频为你领路，还等什么，看啊！！！

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

其实没有标准答案的问题

我的代码的进化历史

文末友情宣传

相关推荐