生信编程14.多个探针对应一个基因,肿么办?
有一些五六年前的学生们都成长为了各个生物信息学相关公司的小领导,而且他们都有了自己的公众号,知乎号,也算是一番人物。最近他们跟我反馈面试找不到或者说很难直接考核筛选到认真干活的生信工程师,挺有意思的。让我想起来了早在生信技能树论坛创立之初我为了引流,而规划的200个生信工程师面试题。值得继续分享:
当多个探针对应一个基因的时候,可以采取不同的方法,如取最大值,平均值和中位数。
CLL
,里面有CLL慢性淋巴癌基因的表达数据。第二个为hgu95av2.db
,包含36个将芯片探针ID转换为其他主流ID的映射数据。library('CLL')
library('hgu95av2.db')
data(sCLLex)
head(sCLLex)#查看数据
共有22个样本
选取前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)]))
)
分开理解:
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
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)))