筛选差异基因，倍数法，t检验，SAM

一、 实验内容

利用SAM筛选差异基因

比较筛选差异基因的方法：倍数法，t检验，SAM

二、 实验目的

掌握筛选差异基因的方法，为了探究倍数法，t检验，SAM这三种筛选差异基因的方法异同和优越性，为我们后期自己做项目寻求最优的方法

原理

倍数法(FC)

FC=Xi/Xc

T检验

H0:Gene i表达没有差异

P<阈值,拒绝H0

SAM（图1）

图1

三、 实验数据工具及步骤

实验数据

从GEO数据库下载GSE52327_series_matrix.txt,样本均是乳腺癌，样本分类按照表达醛脱氢酶（ALDH）在原代人乳房的细胞群含量，实验组是ALDH+，对照组是ALDH-

label.txt,

GPL570.txt平台信息具体三列,probe ID, Gene Symbol, ENTREZ_GENE_ID

工具

R.3.6.3(limma,samr)

步骤

SAM

[1] 下载数据并导入

[2] 根据平台文件，探针对应基因，一对多，一对空，多对多去掉，多对一取均值

[3] 根据注释信息分组并添加标签，ALDH+（8个），ALDH-（8个）,（1,0）

[4] 利用SAM法筛选差异基因,得到分析结果

三种方法比较

[1] 根据上面探针对应基因后的表达谱用FC和t-test进行差异表达分析

[2] 设定不同的阈值,|logFC|<1,2 p<0.1,0.01,0.05,0.001,分别得到差异基因个数

[3] 比较三种方法,利用Venn图展示出来(|logFC|<2,p<0.1)

[4] 画火山图,将差异基因可视化(y:log2FC,x:adjust_p)

[5] 热图以t-test为例

#数据处理，探针对基因

setwd("F:\\实验\\案例分析\\生物信息学案例分析\\实验一")

probe_exp<-read.table("GSE52327_series_matrix.txt",header=T,sep="\t",row.names=1) #读基因表达文件 54675*16

probeid_geneid<-read.table("GPL570_probe_geneid.txt",header=T,sep="\t") #读探针文件

label<-read.table("label.txt",as.is=T,header=T,sep="\t")

probe_name<-rownames(probe_exp)

loc<-match(probeid_geneid[,1],probe_name)#按照probe_exp的probe进行匹配（长度与probeid_geneid[,1]相同，返回平台文件probe在表达谱probe_name位置）

probe_exp<-probe_exp[loc,]#确定能匹配上的probe表达值（把探针对基因（空）删除了）

raw_geneid<-as.numeric(as.matrix(probeid_geneid[,3]))

index<-which(!is.na(raw_geneid)) #找出有geneid的probeid并建立索引

geneid<-raw_geneid[index] #提取有geneid的probe

exp_matrix<-probe_exp[index,] #找到每个geneid的表达值

geneidfactor<-factor(geneid)

gene_exp_matrix<-apply(exp_matrix,2,function(x) tapply(x,geneidfactor,mean)) #多个探针对应1个基因的情况，取平均值

rownames(gene_exp_matrix)<-levels(geneidfactor) #geneid作为行名

geneid<-rownames(gene_exp_matrix)

gene_exp_matrix2<-cbind(geneid,gene_exp_matrix)

write.table(gene_exp_matrix2,file="geneid_exp.txt",sep="\t",quote=F,row.names=F)

#分组

index1<-which(label=="aldh status: ALDH+")#8

label[index1,]<-1

index2<-which(label=="aldh status: ALDH-")#8

label[index2,]<-0

#数据准备

label<-as.numeric(as.matrix(label))

exp<-gene_exp_matrix

gid<-rownames(gene_exp_matrix)

gid<-as.numeric(gid)

if(!require("BiocManager"))

install.packages("BiocManager",update = F,ask = F)

BiocManager::install("samr")

library("samr")

###SAM

source("detec.slab2.r")

source("samr.compute.delta.table.zhang.r")

rownames(exp)<-gid

colnames(exp)<-label

label[label!=1]<-2

data=list(x=exp,y=label, geneid=rownames(exp),genenames=paste("g",as.character(1:nrow(exp)),sep=""), logged2=TRUE)

samr.obj<-samr(data, resp.type="Two class unpaired", nperms=1000)

source("num_sig_zhang.R") #计算检验统计量

ensemble_zhang<-num_sig_zhang(data,samr.obj,fdr=c(0.001,0.01,0.05,0.10))

#得到不同阈值下:差异基因个数,检验统计量的参数,d(i),基因表达情况，写出

write.csv(ensemble_zhang[[1]],file="br_sig_num.csv")

write.csv(ensemble_zhang[[2]],file="br_sig_gene.csv")

write.csv(samr.obj$tt,file="br_d_stat.csv")

write.csv(ensemble_zhang[[3]],file="br_cut_inf.csv")

2.

#FC&T检验

setwd("F:\\实验\\案例分析\\生物信息学案例分析\\实验一")

exp<-read.table("geneid_exp.txt",header=T,sep="\t",row.names=1)

label<-read.table("label.txt",as.is=T,header=T,sep="\t")

library(limma)

rt<-exp

index1<-which(label=="aldh status: ALDH+")#8

index2<-which(label=="aldh status: ALDH-")#8

hi_exp<-rt[,index1]

lo_exp<-rt[,index2]

rt<-cbind(hi_exp,lo_exp)

#differential

class<-c(rep("ALDH+",8),rep("ALDH-",8))

design<-model.matrix(~factor(class))

colnames(design)<-c("hi","lo")

fit<-lmFit(rt,design)

fit2<-eBayes(fit)

allDiff=topTable(fit2,adjust='fdr',coef=2,number=200000)

write.table(allDiff,file="limmaTab.txt",sep="\t",quote=F)

#通过不同阈值,得到差异基因列表

foldChange=1

padj=0.05

#UP log2FC>=1

#down log2FC<=-1

foldChange=1

padj=0.05

diff_FC<-allDiff[with(allDiff, (logFC>foldChange|logFC<(-foldChange))),] #1518

diffup<-allDiff[with(allDiff, ((logFC>foldChange))),] #910

diffdown<-allDiff[with(allDiff,((logFC<(-foldChange)))),] #608

diff_p<-allDiff[with(allDiff, (adj.P.Val

dim(diff_p)

dim(diffdown)

dim(diffup)

dim(diff_FC)

write.table(diff_p,file="diff_p.txt",sep="\t",quote=F)

write.table(diff_FC,file="diff_FC.txt",sep="\t",quote=F)

write.table(diffup,file="up.txt",sep="\t",quote=F)

write.table(diffdown,file="down.txt",sep="\t",quote=F)

#Venn图，查看三种方法结果的交叠情况

gid<-rownames(rt)

str(ensemble_zhang[[2]])

diffgene_sam<-gid[which(ensemble_zhang[[2]]!=0)]

a<-as.numeric(rownames(diff_FC))

b<-as.numeric(rownames(diff_p))

library("VennDiagram")

venn.diagram(list(FoldChange_2=a,Ttest_0.05=b,SAM_0.1=diffgene_sam),resolution = 300,

imagetype = "tiff",

alpha=c(0.5,0.5,0.5),fill=c("red","yellow","blue"),

cat.fontface=4,fontfamily=3,

main="FoldChange&T-test&SAM",

main.cex = 2, main.fontface = 2, main.fontfamily = 3,

filename = "Venn.tif")

#volcano

pdf(file="vol.pdf")

xMax=max(-log10(allDiff$adj.P.Val)) #x轴范围为P.value最值

yMax=max(abs(allDiff$logFC)) #y轴范围为foldchange最值

plot(-log10(allDiff$adj.P.Val),allDiff$logFC,xlab="adj.P.Val",ylab="logFC",main="volcano") #把所有的foldchange值和P.value打黑点

diffSub1=subset(allDiff,allDiff$adj.P.Val<0.05&(allDiff$logFC)>1)

diffSub2=subset(allDiff,allDiff$adj.P.Val<0.05&(allDiff$logFC)<(-1))

points(-log10(diffSub1$adj.P.Val),diffSub1$logFC,pch=20,col="red",cex=0.4)

points(-log10(diffSub2$adj.P.Val),diffSub2$logFC,pch=20,col="green",cex=0.4)

#打红点绿点

abline(h=0,lty=2,lwd=3) #画中线

dev.off()

###热图

setwd("F:\\实验\\案例分析\\实验五")

exp<-read.table("geneid_exp.txt",header=T,sep="\t",row.names=1)

label<-read.table("label.txt",header=T,sep="\t")

index1<-which(label==1)#93

index2<-which(label==0)#54

exp1<-exp[,index1]

exp2<-exp[,index2]

exp<-cbind(exp1,exp2)

loc<-match(a,rownames(exp))

exp11<-exp[loc,]#13*4113

annotation_col= data.frame(sample= factor(rep(c("aldh status: ALDH+", "aldh status: ALDH-"), c(8,8)))) #列名

rownames(annotation_col) = colnames(exp)

ann_colors = list(value = c("aldh status: ALDH+"= "#7570B3","aldh status: ALDH-"="#E574AE"))

pheatmap(exp1, annotation_col = annotation_col,

annotation_colors = ann_colors,color=colorRampPalette(colors=c("navy", "white","firebrick3"))(50),

main="heatmap of diffgene by t-test")

实验结果：

1.SAM 结果:得到4个表格如图2

图 2 阈值在在p的0.001,0.01,0.05,0.10,得到的差异基因个数

检验统计量的参数,分别p的阈值是0.001,0.01,0.05,0.10，Delta是d(i)-dE(i),med false pos是错误率中位数，第三列90%置信区间错误率，called是筛选出的差异基因数，第五列错误发现率中位数，第六列90%置信区间错误发现率，最后两列统计量阈值（最大最小）

d(i)值

基因表达情况(正数上调,负数下调)

方法比较

三种方法比较

A:Venn图,查看基因交叠情况,|log2FC|>2, T-test p<=0.1/0.05, SAM p<=0.1 B:条形图 C~F: 热图，三种方法以及三种方法的交集得出的差异基因表达情况

火山图(差异基因可视化,红色上调,绿色下调)

结果分析:

从上面所有统计结果来看,三种方法中,FC法筛选出的差异基因最多,所以更没有参考性,SAM和t检验结果有较高的一致性,但在运算上,SAM算法要更加优秀,筛选出的差异基因也会更准确,但是运算量大,当然,FC值可以用来判断基因表达上调还是下调.所以,在时间允许下筛选差异基因首选SAM法