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

Posted on by

前些天我发现了乳腺癌领域的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.