生信编程14.多个探针对应一个基因,肿么办?

有一些五六年前的学生们都成长为了各个生物信息学相关公司的小领导,而且他们都有了自己的公众号,知乎号,也算是一番人物。最近他们跟我反馈面试找不到或者说很难直接考核筛选到认真干活的生信工程师,挺有意思的。让我想起来了早在生信技能树论坛创立之初我为了引流,而规划的200个生信工程师面试题。值得继续分享:

200个生信工程师面试考题

当多个探针对应一个基因的时候,可以采取不同的方法,如取最大值,平均值和中位数。

我们首先来看取平均值的,需要用到的两个包,第一个为CLL,里面有CLL慢性淋巴癌基因的表达数据。第二个为hgu95av2.db,包含36个将芯片探针ID转换为其他主流ID的映射数据。

library('CLL')
library('hgu95av2.db')
data(sCLLex)
head(sCLLex)#查看数据

共有22个样本

ExpressionSet (storageMode: lockedEnvironment) assayData: 6 features, 22 samples element names: exprs protocolData: none phenoData sampleNames: CLL11.CEL CLL12.CEL ... CLL9.CEL (22 total) varLabels: SampleID Disease varMetadata: labelDescription featureData: none experimentData: use 'experimentData(object)' Annotation: hgu95av2

选取前8个样本对其进行分析

sCLLex=sCLLex[,1:8]
group_list = sCLLex$Disease #以疾病的状态为分组
exprset = exprs(sCLLex)
head(sCLLex)

CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL 1000_at    5.743132  6.219412  5.523328  5.340477  5.229904 1001_at    2.285143  2.291229  2.287986  2.295313  2.662170 1002_f_at  3.309294  3.318466  3.354423  3.327130  3.365113 1003_s_at  1.085264  1.117288  1.084010  1.103217  1.074243 1004_at    7.544884  7.671801  7.474025  7.152482  6.902932 1005_at    5.083793  7.610593  7.631311  6.518594  5.059087

取平均值:

exprset = as.data.frame(exprset)
probe2symbol = toTable(hgu95av2SYMBOL)
dat = merge(exprset, probe2symbol,by='probe_id')
results=t(
  sapply(split(dat,dat$symbol), function(x) colMeans(x[,2:(ncol(x)-1)]))
)
 

重点在于最后一条代码:

results_mean=t(
  sapply(split(dat,dat$symbol), function(x) colMeans(x[,2:(ncol(x)-1)]))
)

分开理解:

  1. split(dat,dat$symbol),这一部分表示的是根据symbol来对数据进行分组,最终得到的是一个列表,列表中的每个元素代表的是同一个基因symbol所对应的探针ID以及基因表达量,如下:

$CAPN7 probe_id CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL 2213 32516_at  4.550316  4.365438  4.096016  4.236659  4.465431 2214 32517_at  5.657783  5.610177  5.277241  5.766585  4.857069 CLL16.CEL CLL17.CEL CLL18.CEL symbol 2213  4.669496  4.682865  5.026866  CAPN7 2214  5.081120  4.872532  5.709810  CAPN7

  1. function(x) colMeans(x[,2:(ncol(x)-1)])) 这个函数接收sapply函数传入的每个列表元素--即同一基因对应的相同探针,然后计算每一列的平均值。

我们可以使用同样的代码取最大值:

results_max=t(sapply(split(dat,dat$symbol), function(x) apply(x[,2:(ncol(x)-1)], 2, max)))

文末友情推荐

(0)

相关推荐

  • EBMT21 Notes: Lymphoma 篇

    CLL  随着新药崛起,CLL中的移植在下降,但是BTK后也亟需新的疗法,目前发现CLL患者中T细胞有些缺陷比如耗竭标志物升高,这可能也是CAR T在CLL中不如DLBCL的原因 Tisa-cel和A ...

  • R语言GEO数据处理(四)

    # 3. id转换 ----------------------------------------------------------------- ##方法一:使用R包转换 index = gse ...

  • 生信编程3.hg38每条染色体的基因、转录本分布

    有一些五六年前的学生们都成长为了各个生物信息学相关公司的小领导,而且他们都有了自己的公众号,知乎号,也算是一番人物.最近他们跟我反馈面试找不到或者说很难直接考核筛选到认真干活的生信工程师,挺有意思的. ...

  • 生信编程直播第七题:写超几何分布检验!

    下载数据 切换到工作目录:cd d/生信技能树-视频直播/第七讲 kegg2gene(第六讲kegg数据解析结果) 暂时不用新的kegg注释数据为了能够统一答案 差异基因list和背景基因list 关 ...

  • 生信编程直播课程优秀学员作业展示1

    题目 人类基因组外显子区域长度 学员:x2yline 具体题目详情请参考生信技能树论坛 题目数据来源为:ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_huma ...

  • 生信编程直播课程优秀学员学习心得及作业展示3

    学习感悟 首先说明一下,我不算是完全从0开始学习,因为生物的知识和python的语言之前都知道一点,但说实话,我的python距离真正的实践还差的很远,也没有常用所以基本忘完. 真的很感谢群主和老师们 ...

  • 生信编程直播课程优秀学员作业展示2

    题目:hg19基因组序列的一些探究 学员:x2yline 具体题目详情请参考生信技能树论坛 数据来源:http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bi ...

  • 生信编程直播第0题-生信编程很简单!

    貌似顺序有点怪异,我们都把1~8题给讲解完了,现在才开始第0题,很明显,这个是临时加入的,因为有很多没有编程基础的朋友也开始跟着我们学习了. 还不知道怎么回事的先查看历史题目: 生物信息学技能面试题( ...

  • 生信编程直播第9题-根据指定染色体及坐标得到参考碱基

    还不知道怎么回事的先查看历史题目: 生信编程直播第0题-生信编程很简单! 生物信息学技能面试题(第1题)-人类基因组的外显子区域到底有多长 生物信息学技能面试题(第2题)-探索人类基因组序列 生物信息 ...

  • 生信编程直播第11题:把文件内容按照染色体分开写出

    问题描述这个需求很常见,因为一般生物信息学数据比较大,比如sam,vcf,或者gtf,bed都是把所有染色体综合在一起的文件.如果想根据染色体把大文件拆分成小的文件呢?比如:ftp://ftp.ncb ...

  • 生信编程直播第12题:json格式数据的格式化

    json数据大家统一用我给的测试数据,自己在浏览器打开下载:http://biotrainee.com/jbrowse/JBrowse-1.12.1/sample_data/json/modencod ...