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

生信技能树

发布于 2020-04-22 16:44:01

1.7K00

代码可运行

文章被收录于专栏：生信技能树生信技能树

运行总次数：0

代码可运行

前些天我发现了乳腺癌领域的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第一次出现的信息

这个代码非常好用，你一定要学习哦！

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2020-04-22，如有侵权请联系 cloudcommunity@tencent.com 删除

硬件开发