专栏首页生信技能树多个探针对应同一个基因取最大值的代码进化历史

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

我的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站:

表达芯片的公共数据库挖掘系列推文感兴趣的也可以去看看;

本文分享自微信公众号 - 生信技能树(biotrainee),作者:生信技能树

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-10-22

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 这个WGCNA作业终于有学徒完成了!

    共画了3张热图,最后一张热图展示如下图,与原文对比'Ligamentocyte'和'Chondrocyte'相比较其他组是高表达的。

    生信技能树
  • TCGA的28篇教程-GTEx数据库-TCGA数据挖掘的好帮手

    这个时候我们就需要想办法加大正常组织测序样本量,既然TCGA数据库没有,我们就从其他数据库着手。

    生信技能树
  • GEO数据挖掘-第一期-胶质母细胞瘤(GBM)

    lncRNAs PVT1 and HAR1A are prognosis biomarkers and indicate therapy outcome for...

    生信技能树
  • WCF服务发布到IIS时候,只能根据hostname访问,不能根据IP地址访问的解决办法

    本文转载:http://www.cnblogs.com/deerbox/archive/2013/05/13/3076248.html

    跟着阿笨一起玩NET
  • Jenkins(一)

    Jenkins2.x支持pipeline as code,可以通过代码来「描述」部署流水线。

    zx钟
  • HTML5设计原理(上)

    今天我想跟大家谈一谈HTML5的设计。主要分两个方面:一方面,当然了,就是HTML5。我可以站在这儿只讲HTML5,但我并不打算这样做,因为如果你想了解HTML...

    RP道貌不岸然
  • Newbe.Mahua 1.11 支持热更新

    从此版本开始,支持插件热更新的 Newbe.Mahua SDK 将给开发者带来更加丝滑的体验。

    newbe36524
  • Math-Model(五)正交分解(QR分解)

    矩阵的正交分解又称为QR分解,是将矩阵分解为一个正交矩阵Q和一个上三角矩阵的乘积的形式。

    Pulsar-V
  • Face++人脸识别及身份证、银行卡扫描

    年前忙着赶项目,也没时间更新,现在告一段落,因为是贷款类项目,涉及到审批验证等信息,不可避免的使用到了人脸识别、身份证验证、银行卡扫描等相关技术,这里就来聊聊这...

    iOSSir
  • php中ini_set的函数修改php.ini的参数

    对于使用虚拟空间的站长来说,PHP的一些配置是很难更改的,不过PHP给我们提供一个ini_set函数,可以临时修改PHP配置文件php.ini的设置,无需打开此...

    周俊辉

扫码关注云+社区

领取腾讯云代金券