1. 什么是差异表达基因

在不同组织中表达发生明显变化的基因

是导致细胞状态发生变化的关键基因

是表达谱分析的主要对象

2. 寻找差异表达基因的两种方法

倍数变化阀值(一般设置为2倍)

具体方法：

找出所有基因的表达变化率

按照表达变化率排序

上调两倍或者下调两倍算作差异表达基因

适合条件：实验重复数极少

差异基因数目比例(top5%, 即最上调的2.5%，最下调的2.5%)

用假设检验来做

3. R 语言实现

该研究者试图通过Infliximab治疗溃疡性结肠炎来比较有应答与无应答组间的基因表达变化情况

3.1 材料准备

1.下载“GSE12251” 23个项目

image.png

制作一个target.txt 文件

相当于记录临床信息，12 个为应答组，11个为无应答组

image.png

下载“anonotation.csv”

去 affy 官网下载，下载下来的是"HG-U133_Plus_2.na36.annot.csv", 需要整理如下

image.png

3.2 R 代码

library(affy)

library(limma)

##import phenotype data

phenoData = read.AnnotatedDataFrame('target.txt')

pheno = pData(phenoData)

View(pheno)

#import Annotaion

anno = read.csv("Annotation.csv",head=T)

View(head(anno))

##RMA normalization

#eset.rma = RMA(Data)

eset.rma

datExpr = exprs(eset.rma)

#补充缺失值

library(impute)

#KNN法计算缺失值

imputed_gene_exp = impute.knn(datExpr,k=10,rowmax = 0.5,

colmax=0.8,maxp =3000, rng.seed=362436069)

datExpr2 = imputed_gene_exp$data

write.table(datExpr2,file="Expdata.txt",sep="\t")

获得一个校正过后的差异表达数据

image.png

#######

#样本分组

Group = factor(pheno$group,levels=c('Responder','Nonresponder'))

design = model.matrix(~0+Group)

colnames(design)

design

#线性模型拟合

fit

#构建比对模型，比较两个条件下的表达数据

contrast.matrix

levels=design)

#####################

####################

library(xlsx)

library(futile.logger)

#比对模型进行差值计算

fit2

#贝叶斯检验

fit2

#找出差异基因检验结果并输出符合条件的结果

用responser - Nonresponder 看正表达和负表达

diff = topTable(fit2,adjust.method="fdr",coef=1,p.value=0.05,

lfc=log(2,2),number=5000,sort.by = 'logFC')

#diff

diffGene = annoGene.Symbol[match(rownames(diff),anno$ID_REF)]

diff$ID_REF = rownames(diff)

diff = diff[,c(8,7,1:6)]

diff = diff[diff$Gene != '---',]

#output

write_xlsx2(diff, 'DEG.xlsx', sheetName=colnames(contrast.matrix)[1], colNames=True, rowNames=False, append=True)

最后找到差异表达基因列表，供后续分析

image.png

4. 整合数据分析

上面的实例是从一个批次的样本，只需要归一化和背景校正即可。

如果在公共数据库中获取多个批次的数据样本，并希望消除这些批次间的差异（即 batch effects），那么应该采取哪些方法来实现这一目标？

批间差一般由不同人员，不同环境，不同时间，不同机器，不同批次引起的

##sva--combat

library(sva) # 神奇，这个包竟然是实验室的师兄写的

library(pamr)

batch = pheno[,c('batch')]

batch

modcombat = model.matrix(~1,data=pheno)

combat_edata = ComBat(dat=datExpr2,batch=batch,

mod=modcombat,par.prior=T,

prior.plot=T)

write.table(combat_edata,file="Expdata.txt",sep="\t")

image.png