大鼠表达量芯片数据处理

大鼠表达量芯片数据处理 by 生信技能树

前面咱们组建了交流群：让天下没有难处理的表达量芯片，帮助大家解决了各种稀奇古怪的芯片。

但是绝大部分小伙伴其实是基础知识不牢固，有一些明明是很简单的芯片，仍然是有小伙伴提问，比如：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE10775

Affymetrix公司的芯片第一巨头了，很简单的一个芯片：GPL1355 [Rat230_2] Affymetrix Rat Genome 230 2.0 Array，可以看到是12个样品，分成了4个组：

 GSM271906 Control_1
GSM271907 Control_2
GSM271908 Control_3
GSM271909 NRG+IGF_1
GSM271910 NRG+IGF_2
GSM271911 NRG+IGF_3
GSM271912 NRG_1
GSM271913 NRG_2
GSM271914 NRG_3
GSM271915 IGF_1
GSM271916 IGF_2
GSM271917 IGF_3

我们先看看数据集最开始的文章：Extracellular growth factors and mitogens cooperate to drive mitochondrial biogenesis. J Cell Sci 2009 Dec 15;122(Pt 24):4516-25. PMID: 19920079

简单的看了看，文章做了3次差异分析，所以有3个火山图：

3个火山图

因为是十几年前的数据集和文章，所以图表都很粗糙，我们使用标准代码下载表达量矩阵进行其中一个差异分析即可，首先解压我们的标准代码：

标准代码在交流群：让天下没有难处理的表达量芯片，群公告置顶，感兴趣的自己去拿哈！

然后去里面慢慢修改，主要是下载文件的时候肯定是需要修改 GSE号码啊，每个人分析的数据集不一样：

rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)

getOption('timeout')
options(timeout=10000)

library(AnnoProbe)
library(GEOquery) 
gset <- geoChina("GSE10775")
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 
# 这里一定要人工干预，每个数据集项目的分组不一样
# 是主观判断，自己选择  
group_list=  gsub('_[0-9]','',pd$title)
table(group_list)
 
dat[1:4,1:4]  
dim(dat)
a  
ids=idmap('GPL1355','soft')
head(ids)
ids=ids[ids$symbol != '',]
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['Gapdh',]
save(dat,group_list,phe,file = 'step1-output.Rdata')

然后简单的质量控制，就可以看看表达量矩阵是否合格，如下所示；

cor_top500

可以看到NRG这个分组，跟control的差异不明显，所以如果是NRG去跟control比较，应该是差异基因数量少，另外两个分组就跟control的差异比较大。

文章是3个处理组，分别和对照进行差异分析，我们随便测试一个即可：

library(limma)
design=model.matrix(~factor( group_list ))
fit=lmFit(dat,design)
fit=eBayes(fit)
## 上面是limma包用法的一种方式 
options(digits = 4) #设置全局的数字有效位数为4
#topTable(fit,coef=2,adjust='BH') 
colnames(fit$design)
# 这里的 coef=2 就是选择了 IGF 处理，去和control对比 
# [1] "(Intercept)"               "factor(group_list)IGF"     "factor(group_list)NRG"     "factor(group_list)NRG+IGF"
deg = topTable(fit,coef=2,adjust='BH', n=Inf)

然后顺手画一个火山图：

library(EnhancedVolcano)
EnhancedVolcano(deg,
                lab =  rownames(deg),
                x = 'logFC',pCutoff = 0.01,
                y = 'P.Value')

比原文好看：

顺手画一个火山图

最后的上下调基因通路注释也是很简单：

上下调基因通路注释

这个时候就需要看原文了，文章：Extracellular growth factors and mitogens cooperate to drive mitochondrial biogenesis. J Cell Sci 2009 Dec 15;122(Pt 24):4516-25. PMID: 19920079，为什么有这样的上下调基因及其对应的生物学功能。

GEO数据库里面的表达量芯片数据处理，主要的难点是表达量矩阵获取和探针的基因名字转换，搞定后只需要一定的生物学背景对数据进行合理的分组后就是标准的差异分析，富集分析。主要是参考我八年前的笔记：

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：

数据挖掘（GEO,TCGA,单细胞）2022年暑期班（收官之作），快速了解一些生物信息学应用图表
生信入门课-2022年暑期班（收官之作），你的生物信息学第一课

ixxmu / mp_duty

大鼠表达量芯片数据处理 #2423

大鼠表达量芯片数据处理 by 生信技能树