一个基因上面有多个探针最后只能选一个吗
最近学员提出来了一个蛮古老的表达量芯片数据集的讨论，因为 它是做了这个PPARα的基因敲除，但是学员在分析表达量矩阵做差异的时候发现PPARα本身其实并没有统计学显著的差异表达。
数据集是：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE8292
学员选取的样品是control条件下，而不是 a single dose of WY14643 (400 μl of 10 mg/ml WY14643) 药物处理的 ：
GSM205769 liver_wildtype_WY14643_6h_rep1
GSM205770 liver_wildtype_WY14643_6h_rep2
GSM205771 liver_wildtype_WY14643_6h_rep3
GSM205772 liver_wildtype_WY14643_6h_rep4
GSM205778 liver_PPARa-knockout_WY14643_6h_rep1
GSM205779 liver_PPARa-knockout_WY14643_6h_rep2
GSM205780 liver_PPARa-knockout_WY14643_6h_rep3
GSM205781 liver_PPARa-knockout_WY14643_6h_rep4
GSM205782 liver_PPARa-knockout_WY14643_6h_rep5
同样的control条件下，我们比较 野生型和PPARα的基因敲除两个分组， 理论上肯定是表达量有统计学显著的差异。
首先整理基因芯片表达量矩阵和探针注释以及分组信息
如下所示的标准代码：
rm(list = ls()) ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
getOption('timeout')
options(timeout=10000)
library(AnnoProbe)
library(GEOquery)
gset <- geoChina("GSE8292")
gset
gset[[1]]
a=gset[[1]] #
dat=exprs(a) #a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
dim(dat)#看一下dat这个矩阵的维度
dat[1:4,1:4] #查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列
boxplot(dat[,1:4],las=2)
dat=log2(dat)
boxplot(dat[,1:4],las=2)
library(limma)
dat=normalizeBetweenArrays(dat)
boxplot(dat[,1:4],las=2)
pd=pData(a)
#通过查看说明书知道取对象a里的临床信息用pData
## 挑选一些感兴趣的临床表型。
library(stringr)
# pd=pd[1:43,]
# dat=dat[,1:43]
phe=pd
# 这里一定要人工干预，每个数据集项目的分组不一样
# 是主观判断，自己选择
# 比如这里我们仅仅是关心 没有药物处理的
kp=grepl('control',phe$title)
table(kp)
phe=phe[kp,]
dat=dat[,kp]
group_list=ifelse(grepl('wildtype',phe$title),'WT','KO')
table(group_list)
dat[1:4,1:4]
dim(dat)
a
ids=idmap('GPL1261','soft')
head(ids)
ids=ids[ids$symbol != '',]
cg=dat[ids[ids$symbol=='Ppara',1],]
boxplot(cg[1,] ~ group_list)
boxplot(cg[2,] ~ group_list)
boxplot(cg[3,] ~ group_list)
save(dat,group_list,ids,file = 'before_remove_dup_id.Rdata')
可以看到这个 PPARα 基因其实有3个对应的探针 ：
PPARα 基因其实有3个对应的探针
我们授课提到的默认流程是，多个探针就选取表达量最大的探针作为这个基因的代表即可，所以如下所示：
dat=dat[rownames(dat) %in% ids$ID,]
ids=ids[match(rownames(dat),ids$ID),]
head(ids)
colnames(ids)=c('probe_id','symbol')
ids$probe_id=as.character(ids$probe_id)
rownames(dat)=ids$probe_id
dat[1:4,1:4]
ids=ids[ids$probe_id %in% rownames(dat),]
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第一次出现的信息
dat['Actb',]
dat['Gapdh',]
save(dat,group_list,phe,file = 'step1-output.Rdata')
这个时候默认进行表达量矩阵质量控制，进行差异分析，代码很常规就不列出来了。主成分图，热图，火山图，都挺好的：
火山图
但是可以看到：
source('step2-check.R')
source('step3-DEG.R')
load('deg.Rdata')
deg['Ppara',]
这个基因虽然在敲除的组里面确实是表达量下降了，但是居然统计学不显著，在边缘疯狂试探 ：
> deg['Ppara',]
logFC AveExpr t P.Value adj.P.Val
Ppara -0.4464635 10.28399 2.079013 0.06913026 0.3137238
但是，如果我们前面并不进行基因的探针过滤步骤，直接进行差异分析，如下所示：
rm(list = ls()) ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load(file = 'before_remove_dup_id.Rdata')
library(limma)
design=model.matrix(~factor( group_list ,levels = c('WT','KO')))
fit=lmFit(dat,design)
fit=eBayes(fit)
## 上面是limma包用法的一种方式
options(digits = 4) #设置全局的数字有效位数为4
#topTable(fit,coef=2,adjust='BH')
deg = topTable(fit,coef=2,adjust='BH', n=Inf)
deg[ids[ids$symbol=='Ppara',1],]
可以看到，这个基因的3个探针，都显示出来了：
> deg[ids[ids$symbol=='Ppara',1],]
logFC AveExpr t P.Value adj.P.Val B
1439675_at -0.4465 10.284 -2.068 7.030e-02 0.331133 -4.956
1449051_at -1.0463 9.057 -3.883 4.126e-03 0.086350 -2.107
1457721_at -1.7000 7.427 -8.500 1.918e-05 0.003281 3.427
如果多个探针就选取表达量最大的探针作为这个基因的代表即可，我们针对这个Ppara 就会挑选到1439675_at这个探针，以至于它虽然下调，但是并不达标。
但是如果选择了 1457721_at这个探针，作为Ppara基因的代表，它就比较好的下调啦。
几个思考
选取表达量最大的探针作为这个基因的代表合理吗？
PPARα的基因敲除意味着表达量芯片或者转录组测序里面，它表达量都会下降吗？
学徒作业
找到同一个基因敲除的表达量芯片和转录组测序数据，一般来说只能是从明显基因里面找啦，下载其对应的表达量芯片和转录组测序数据做差异分析，看看作者敲除的基因是否确实有表达量下降的情况发生！
腾讯云开发者
