多个探针对应同一个基因取最大值的代码进化历史

生信技能树

发布于 2019-10-25 00:44:26

2.6K0

发布于 2019-10-25 00:44:26

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

我的GEO芯片数据分析教程本来就是为粉丝写的，基本上就是生信菜鸟团QQ群的诸位问什么，我就临时搜索整理讲解那个知识点，非常融洽，目录如下：

第一讲：GEO，表达芯片与R
第二讲：从GEO下载数据得到表达量矩阵
第三讲：对表达量矩阵用GSEA软件做分析
第四讲：根据分组信息做差异分析
第五讲：对差异基因结果做GO/KEGG超几何分布检验富集分析
第六讲：指定基因分组boxplot指定基因list画热图
第七讲：根据差异基因list获取string数据库的PPI网络数据
第八讲：PPI网络数据用R或者cytoscape画网络图
第九讲：网络图的子网络获取
第十讲：hug genes如何找

最近全国巡讲的学员又问到了多个探针对应同一个基因取最大值类似的问题，我们的斯老师找到了我三年前的博客：多个探针对应一个基因，取平均值或者最大值我看到里面的留言很有趣：

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

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

因为芯片技术已经被时代抛弃，所以我们这里也不继续深究了，我感兴趣的是我的代码进化路程

　第一版，使用split结合 sapply

下面代码写于2016年6月左右

library('hgu95av2.db')
library(CLL)
data(sCLLex)
sCLLex=sCLLex[,1:8] ## 样本太多，我就取前面8个
group_list=sCLLex$Disease
exprSet=exprs(sCLLex)
exprSet=as.data.frame(exprSet)
exprSet$probe_id=rownames(exprSet)
head(exprSet)
probe2symbol=toTable(hgu95av2SYMBOL)
dat=merge(exprSet,probe2symbol,by='probe_id')
results=t(sapply(split(dat,dat$symbol),
　　　　　　　　　function(x) colMeans(x[,1:(ncol(x)-1)])))

这个代码写完我就忘记了这回事。

第二版，使用by函数

下面代码写于2017年6月左右，这个时候因为是临时授课，其实忘记了自己一年前写过这个代码，所以很粗糙的又写了一次：

 table(rownames(exprSet) %in% ids$probe_id)
  dim(exprSet)
  exprSet=exprSet[rownames(exprSet) %in% ids$probe_id,]
  dim(exprSet)

  ids=ids[match(rownames(exprSet),ids$probe_id),]
  head(ids)
  exprSet[1:5,1:5]
 if(F){
   tmp = by(exprSet,ids$symbol,
            function(x) rownames(x)[which.max(rowMeans(x))] )
   probes = as.character(tmp)
   dim(exprSet)
   exprSet=exprSet[rownames(exprSet) %in% probes ,]
   dim(exprSet)

   rownames(exprSet)=ids[match(rownames(exprSet),ids$probe_id),2]
   exprSet[1:5,1:5]
 }

具体的代码注释，可以看我以前学徒的笔记：分组计算描述性统计量函数—by()函数

第三版，使用duplicated和order函数

写完第二个版本的时候，这个生信人的20个R语言习题已经布置给了一百多个学员和学徒，而根据他们的反馈，我这个by函数，计算耗时非常可怕，我仔细检查后，又写了一个代码：

https://github.com/jmzeng1314/5years/blob/master/learn-R/tasks/3-r-20-codes.R

  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)

还找实习生写了代码注释，这样大家理解起来就不困难了。更多数据挖掘视频都在B站：