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

作者: 七禾叶瓣 | 来源:发表于2022-10-16 21:05 被阅读0次

筛选差异基因，倍数法，t检验，SAM
差异表达分析基本原理与方法
DESeq2做差异分析时，组间样本数不一致可以比较么？
热图绘制
差异分析
t检验算法及其在R语言中的实现
#多个数据集差异基因整合#
差异分析
给女朋友写的生统资料_Part18
[生信基础知识]几种常用的差异表达基因识别方法FC,T检验,SA

一、实验内容
利用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<padj)),] #147             
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法

筛选差异基因，倍数法，t检验，SAM
一、实验内容利用SAM筛选差异基因比较筛选差异基因的方法：倍数法，t检验，SAM二、实验目的掌握筛选差异...
差异表达分析基本原理与方法
差异表达分析的本质就是：识别差异表达基因和非差异表达基因筛选方法: 倍数法，t检验法，方差分析，SAM法，信息熵 ...
DESeq2做差异分析时，组间样本数不一致可以比较么？
答案是可以！如果差异基因太多了，可以通过差异倍数进行筛选。参考资料 https://www.biostars.o...
热图绘制
从全部差异基因筛选出一定阈值（FC，P值）筛选以后得到差异基因，再从所有差异基因的表达矩阵中筛选出最终差异基因的表...
差异分析
基因的差异表达，即发现一组在正常样本和患病样本中表达不同的基因。差异基因的筛选方法 1、倍数变化（fold ch...
t检验算法及其在R语言中的实现
定义及分类在进行 t 检验之前让我们先看看它的定义：t 检验法就是在假设检验时利用 t 分布进行概率计算的检验方...
#多个数据集差异基因整合#
加载包，设置差异倍数logFC和P值读取各个数据集的差异基因读取各个数据集的差异基因读取这些差异基因的数据框...
差异分析
差异基因的筛选方法统计检验得到p值比方说你想研究什么基因，可以在全部都是肿瘤组织里面以某个基因分类，分成高表达...
给女朋友写的生统资料_Part18
apply和function 差异基因的检验估计会用到function和apply。不过差异基因表达的functi...
[生信基础知识]几种常用的差异表达基因识别方法FC,T检验,SA
目前常用的几款差异表达进识别方法有：FC、T检验、SAM等一、FC FC(Fold Change)算法是最早用于...