2020-02-05

GSE95070和GSE95426进行联合差异基因分析

下载两个小鼠的数据集，想对他进行联合差异基因分析。有两种思路，一是对各个数据集进行差异分析，然后再进行差异数据取交集，利用韦恩图展现出来。第二个就是除批次效应将两个数据和在一起

GSE95070这个数据集没有所以去官网下载，然后读取。

数据处理

rm(list = ls())
.libPaths(c("C:/Users/22349/R3.6.1", .libPaths()))
.libPaths()
install.packages("devtools")
install.packages("rlang")
library(devtools)
install_github("jmzeng1314/AnnoProbe")
library(AnnoProbe)
suppressPackageStartupMessages(library(GEOquery))
gset1=AnnoProbe::geoChina('GSE95426')
#gset2=AnnoProbe::geoChina('GSE95070')
dat2= read.table(file="./GSE95070_series_matrix.txt.gz",
                header = T,sep = "\t",quote = "",fill = T,
                comment.char = "!")
rownames(dat2)= dat2[,1]
dat2 = dat2[,-1]
library(stringr)
colnames(dat2) <- str_split(colnames(dat2),'[.]',simplify = T)[,2]
eset=gset1[[1]]
dat1 <- exprs(gset1[[1]]) #（提取表达矩阵）
dat1[1:4,1:4]
pd1<- pData(gset1[[1]]) # (提取临床信息)
pd1[1:4,1:4] 
boxplot(dat1,las=2)

dat2=log2(dat2+1)
library(limma)
dat=normalizeBetweenArrays(dat)
boxplot(dat2,las=2)
save(pd1,gpl1,gpl2,dat1,dat2,file = "step1-output.Rdata"

image-20200205163219192.png

image-20200205162846323.png

分组id注释转换（这里闹出了乌龙）

pd1=pd1[,apply(pd1, 2, function(x){
  length(unique(x))>1})]  #缩小范围，从所有临床信息中筛选出含有大于2个元素的信息列
dim(pd1)
apply(pd1,2,table)

table(pd1$source_name_ch1)
control=rownames(pd1[grepl('Control hippocampal tissue',as.character(pd1$source_name_ch1)),])
POCD=rownames(pd1[grepl('POCD hippocampal tissue',as.character(pd1$source_name_ch1)),]) 
dat=dat1[,c(control,POCD)]
class(dat)
dim(dat)
group_list=c(rep('control',length(control)) ,
             rep('POCD',length(POCD))) #分组信息
table(group_list)

eset@annotation
GPL=eset@annotation
library(AnnoProbe)
ids<- idmap(eset@annotation)
ids <- idmap(gpl = "GPL22782", type = "pipe", mirror = "tercent")
exprset<- filterEM(dat1,ids)

image-20200205173643414.png

image-20200205174655111.png

colD=data.frame(group_list=group_list)
rownames(colD)=colnames(exprset)
pheatmap::pheatmap(cor(exprset),
                   annotation_col = colD,
                   show_rownames = F,
                   show_colnames =T)

差异分析

group_list2 <- factor(group_list2)
class(group_list2)

design <- model.matrix(~0+factor(group_list2))
colnames(design)=levels(factor(group_list2))
rownames(design)=colnames(dat2)
library(limma)
####构建比较矩阵,设置用来对比的组别
contrast.matrix<-makeContrasts(contrasts=
                                 c('control-POCD'),levels = design)
######limma三部曲,只需要归一化后的数据、实验矩阵、比较矩阵的输入
fit <- lmFit(dat2,design)
fit1 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit1)
####上面得到了p值等统计的结果，topTable对p值校验，对基因排序
tT <- topTable(fit2, adjust="fdr", number=nrow(fit2))
deg2 <- tT
dim(deg2)
save(deg1,deg2,file='GSE95426-deg.Rdata')
####提取出我们需要的指标
tT <- subset(tT, select=c('adj.P.Val',"P.Value","logFC"))


library(VennDiagram) 
venn.diagram(list( GSE95426=deg1$logFC ,
                   GSE95070=deg2$logFC
                   ), 
             fill=c("red","green"),
             filename="VennDiagram.tiff")

image-20200205203524209.png

这两个基因好像没有交集，这就难办了。后面想做两个数据交集的kegg分析和go分析，得再看看原因

2020-02-05

GSE95070和GSE95426进行联合差异基因分析

推荐阅读更多精彩内容