目录
(3)提取差异表达基因
一、数据准备
1.数据加载
#数据表达数据加载
exp1=read.csv('F:\\bioinformatics\\research program\\materialXXXXXXXXXX\\20221101XXXXXXXXXX\\datafile\\GSE5281\\GSE5281traned.csv')
dim(exp1) #查看数据维度
[1] 20857 162
#加载样本信息数据
sample_info=read.csv('F:\\bioinformatics\\research program\\materialXXXXXXXXXX\\20221101XXXXXXXXXX\\datafile\\GSE5281\\GSE5281 sample.csv')
dim(sample_info) #查看数据维度
[1] 161 122.做分组信息数据
提取样本数据关键信息（分组信息）
sample_data=sample_info[,c(1,2)] #提取样本信息数据的关键数据
head(sample_data) #查看数据
Accession Title
1 GSM238763 EC_affected_1
2 GSM238790 EC_affected_2
3 GSM238791 EC_affected_3
4 GSM238792 EC_affected_4
5 GSM238793 EC_affected_5
6 GSM238794 EC_affected_6
新建修改过的样本分组列
sample_data1=sample_data
sample_data1$group_list=c(rep('AD',87),rep('Control',74))
head(sample_data1) #查看数据
Accession Title group_list
1 GSM238763 EC_affected_1 AD
2 GSM238790 EC_affected_2 AD
3 GSM238791 EC_affected_3 AD
4 GSM238792 EC_affected_4 AD
5 GSM238793 EC_affected_5 AD
6 GSM238794 EC_affected_6 AD删除第二列
sample_info1=sample_data1[,-2] #删除第二列
head(sample_info1) #得到分组信息
Accession group_list
1 GSM238763 AD
2 GSM238790 AD
3 GSM238791 AD
4 GSM238792 AD
5 GSM238793 AD
6 GSM238794 AD查看两个数据长度
length(colnames(exp1)) #表达数据的列数，多了一列gene symbol ID
length(sample_info1[,1]) #样本分组类型的长度
[1] 162
[1] 1613.表达数据样本ID顺序与样本信息数据匹配
表达数据中的列名样本ID号是长成这样的：
colnames(exp1)[2] #随便看一一个
[1] "X.GSM119615."为了匹配，需要先修改一下
list1<-list() #新建一个空列表
cols=colnames(exp1[,2:ncol(exp1)]) #提取表达数据的所有样本ID
cols[2] #确认字符串样式
[1] "X.GSM119616."循环语句进行修改
for (i in c(1:length(cols))) {
list1[i]=substr(cols[i],3,11)
}
head(list1) #查看数据
[[1]]
[1] "GSM119615"
[[2]]
[1] "GSM119616"
[[3]]
[1] "GSM119617"
[[4]]
[1] "GSM119618"
[[5]]
[1] "GSM119619"
[[6]]
[1] "GSM119620"替换修改后的样本ID
exp2=exp1
names(exp2)[2:ncol(exp1)]<-list1 #修改列名查看数据
head(colnames(exp2))
head(sample_info1[,1])
[1] "symbol" "GSM119615" "GSM119616" "GSM119617" "GSM119618"
[6] "GSM119619"
[1] "GSM238763" "GSM238790" "GSM238791" "GSM238792" "GSM238793"
[6] "GSM238794"确认两个集合相等
#样本ID比对
FT=colnames(exp2[,2:ncol(exp2)]) %in% list1
a=0
for (j in FT) {
if (j==TRUE) {
a=a+0
}
else if (j==FALSE) {
a=a+1
}
}
a #a=0说明样本对应没有问题解决顺序的问题
###修改样本对应顺序，让两个数据的样本数据一致
sample_info2=sample_info1[match(list1,sample_info1$Accession),]
head(sample_info2 ) #现在两个数据的样本信息数据就一致啦
Accession group_list
88 GSM119615 Control
89 GSM119616 Control
90 GSM119617 Control
91 GSM119618 Control
92 GSM119619 Control
93 GSM119620 Control二、数据探索
(1)查看数据是否经过了log2转换
head(exp2) #发现数据并未经过对数转换发现数据很大，并未经过log2对数转换
#做对数转换
exp3=mutate_if(exp2,is.numeric,funs(log2)) #对数据做对数转换
head(exp3)
dim(exp2)
dim(exp3)
[1] 20857 162
[1] 20857 162(2)查看管家基因的表达量
exp3[exp3[,1] %in% 'GAPDH',] #挺高的：GAPDH，没有问题
exp3[exp3[,1] %in% 'ACTB',] #也挺高的：ACTB，没有问题
(3)检查数据是否经过了分位数归一化
画图数据准备
library(reshape2)
exp3_L=melt(exp3) #宽数据转为长数据
head(exp3_L)
dim(exp3_L)
symbol sample value group
1 DDR1 GSM119615 9.438736 Control
2 PAX8 GSM119615 7.435399 Control
3 GUCA1A GSM119615 5.717456 Control
4 UBA7 GSM119615 6.082226 Control
5 THRA GSM119615 7.710202 Control
6 ESRRA GSM119615 6.864596 Control
>
[1] 3357977 4
names(exp3_L)[2]='sample' #修改列名
head(exp3_L) #查看数据
symbol sample value group
1 DDR1 GSM119615 9.438736 Control
2 PAX8 GSM119615 7.435399 Control
3 GUCA1A GSM119615 5.717456 Control
4 UBA7 GSM119615 6.082226 Control
5 THRA GSM119615 7.710202 Control
6 ESRRA GSM119615 6.864596 Control
####将分组信息加入长数据
library(stringr)
str_detect(sample_info2$group_list,'Control')
g_list=ifelse(str_detect(sample_info2$group_list,'Control')==TRUE,'Control','AD')
length(g_list)
head(g_list) #看看g_list是个什么东西？
[1] 161
[1] "Control" "Control" "Control" "Control" "Control" "Control"
exp3_L$group=rep(g_list,each=nrow(exp3))
head(exp3_L) #查看数据
dim(exp3_L) #查看数据维度
symbol sample value group
1 DDR1 GSM119615 9.438736 Control
2 PAX8 GSM119615 7.435399 Control
3 GUCA1A GSM119615 5.717456 Control
4 UBA7 GSM119615 6.082226 Control
5 THRA GSM119615 7.710202 Control
6 ESRRA GSM119615 6.864596 Control
[1] 3357977 4画箱线图检查分位数归一化的情况
library(ggplot2)
p=ggplot(exp3_L,
aes(x=sample,y=value,fill=group))+geom_boxplot() #fill参数：用分组进行颜色映射
print(p)
接下来，用lmma包的normalizeBetweenArrays函数做批次校正进行归一化
library(limma)
#先将第一列设为行索引
exp4=exp3
rownames(exp4)=exp4[,1]
exp5=exp4[,-1]
exp6=normalizeBetweenArrays(exp5) #进行分位数归一化检查校正后的效果
###查看分位数归一化后的效果
exp6_L=melt(exp6) #宽转长
head(exp6_L) #查看数据
exp6_L$groupx=rep(g_list,each=nrow(exp3))
head(exp6_L)
p0=ggplot(exp6_L,
aes(x=Var2,y=value,fill=groupx))+geom_boxplot() #fill参数：用分组进行颜色映射
print(p0)###查看分位数归一化后的效果
exp6_L=melt(exp6) #宽转长
head(exp6_L) #查看数据
exp6_L$groupx=rep(g_list,each=nrow(exp3))
head(exp6_L)
p0=ggplot(exp6_L,
aes(x=Var2,y=value,fill=groupx))+geom_boxplot() #fill参数：用分组进行颜色映射
print(p0)
上面的图还可以进一步精修
##箱线图精修版
p00=ggplot(exp6_L,
aes(x=Var2,y=value,fill=groupx))+geom_boxplot() #fill参数：用分组进行颜色映射
#去除网格线和背景
p00=p00+theme_bw()+theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = 'black'))
p00 #查看图形
(4)PCA图、层次聚类图
层次聚类图
exp7=exp6
colnames(exp7)=paste(sample_info2$group_list,1:ncol(exp6),sep = '')
head(exp7)
#定义nodePar
nodePar=list(lab.cex=0.6,pch=c(NA,19),cex=0.7,col='blue')
#聚类
hc=hclust(dist(t(exp7))) #t()的意思是转置
#绘图
plot(as.dendrogram(hc),nodePar = nodePar,horiz = TRUE)
好像看不出来什么。。。。。。。。。。。。。。。。。。。。
PCA图
##PCA图
library(ggfortify)
df=as.data.frame(t(exp6)) #转置后就变成了矩阵
dim(df) #查看数据维度
dim(exp6)
df$groupp=sample_info2$group_list #加入样本分组信息
autoplot(prcomp(df[,1:ncol(df)-1]),data=df,colour='groupp') #PCA散点图
大致是分开了
三、差异表达分析
(1)数据准备
需要三个数据：表达矩阵已经有啦（exp6）；还有：分组矩阵，差异比较矩阵
分组矩阵
design=model.matrix(~0+factor(sample_info2$group_list))
colnames(design)=levels(factor(sample_info2$group_list))
head(exp6)
row.names(design)=colnames(exp6)
design #得到分组矩阵：0代表不是，1代表是
AD Control
GSM119615 0 1
GSM119616 0 1
GSM119617 0 1
GSM119618 0 1
GSM119619 0 1
GSM119620 0 1差异比较矩阵
##差异比较矩阵
contrast_matrix=makeContrasts(paste0(c('AD','Control'),collapse = '-'),levels = design)
contrast_matrix #-1和1的意思是Control是用来被比的，AD是用来比的
Contrasts
Levels AD-Control
AD 1
Control -1
> (2)差异分析及可视化
分三步走：lmFit;eBayes;topTable
step1
#step:lmFit
fit=lmFit(exp6,design)
fit2=contrasts.fit(fit,contrast_matrix)step2
#step:eBayes
fit3=eBayes(fit2)step3
step3:topTable
tempoutput=topTable(fit3,coef = 1,n=Inf)
DEG_M=na.omit(tempoutput) #得到差异分析矩阵，重点看logFC和P值
head(DEG_M) #查看数据
logFC AveExpr t P.Value adj.P.Val B
NEAT1 3.872485 9.512774 15.17968 4.053906e-33 8.455232e-29 64.36975
MSI2 1.621265 9.321856 12.00464 2.962532e-24 2.779383e-20 44.43212
EPC1 1.476923 9.952469 11.91114 5.412282e-24 2.779383e-20 43.84258
DTNA 1.873689 12.034226 11.90096 5.779352e-24 2.779383e-20 43.77838
AMER2 1.406314 12.314261 11.85761 7.641583e-24 2.779383e-20 43.50512
NAV2 1.329965 9.248507 11.85058 7.995541e-24 2.779383e-20 43.46082接下来就是将上述结果可视化
热图：选取前40个最显著的基因做热图，查看差异是否真的很显著
library(pheatmap)
f40_gene=head(rownames(DEG_M),40)
f40_subset_matrix=exp6[f40_gene,]
head(f40_subset_matrix)
f40_subset_matrixx=t(scale(t(f40_subset_matrix))) #数据标准化。。。数据标准化和归一化的区别：平移和压缩
pheatmap(f40_subset_matrixx) #出图
火山图
火山图1，没有显示差异基因
#火山图
colnames(DEG_M)
plot(DEG_M$logFC,-log10(DEG_M$P.Value)) #查看图形
火山图2，显示差异基因
DEG=DEG_M
logFC_cutoff=with(DEG,mean(abs(logFC))+2*sd(abs(logFC))) #选差异倍数的阈值
logFC_cutoff
####给差异分析矩阵数据中的基因加上标签
DEG$change=as.factor(ifelse(DEG$P.Value<0.05 & abs(DEG$logFC)>logFC_cutoff,
ifelse(DEG$logFC>logFC_cutoff,'UP','DOWN'),'NOT'))
head(DEG) #查看数据
nrow(DEG[DEG$change=='UP',]) #上调基因数量
nrow(DEG[DEG$change=='DOWN',]) #下调基因的数量
dim(DEG)
dim(exp6)
thttile=paste0('GSE',5281,
'\nCutoff for logFC is ',round(logFC_cutoff,4),
'\nThe number of up gene is ',nrow(DEG[DEG$change=='UP',]),
'\nThe number of down gene is ',nrow(DEG[DEG$change=='UP',]))
thttile #查看准备的火山图标题
g=ggplot(data = DEG,aes(x=logFC,y=-log10(P.Value),color=change))+
geom_point(alpha=0.4,size=1.75)+
theme_set(theme_set(theme_bw(base_size = 20)))+
xlab('log2 fold change')+ylab('-log10 P.Value')+
theme_bw()+theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = 'black'))+
ggtitle(thttile)+theme(plot.title = element_text(size = 15,hjust = 0.5))+
scale_color_manual(values = c('blue','black','red'))
print(g) #输出火山图
(3)提取差异表达基因
head(DEG)
logFC AveExpr t P.Value adj.P.Val B
NEAT1 3.872485 9.512774 15.17968 4.053906e-33 8.455232e-29 64.36975
MSI2 1.621265 9.321856 12.00464 2.962532e-24 2.779383e-20 44.43212
EPC1 1.476923 9.952469 11.91114 5.412282e-24 2.779383e-20 43.84258
DTNA 1.873689 12.034226 11.90096 5.779352e-24 2.779383e-20 43.77838
AMER2 1.406314 12.314261 11.85761 7.641583e-24 2.779383e-20 43.50512
NAV2 1.329965 9.248507 11.85058 7.995541e-24 2.779383e-20 43.46082
change
NEAT1 UP
MSI2 UP
EPC1 UP
DTNA UP
AMER2 UP
NAV2 UP
DEG_MATRIX=exp6[row.names(exp6) %in% rownames(DEG[DEG$change=='UP' | DEG$change=='DOWN',]),]
dim(DEG_MATRIX) #差异表达基因表达数据
[1] 1057 161
UP_DEG=exp6[row.names(exp6) %in% rownames(DEG[DEG$change=='UP',]),]
dim(UP_DEG) #上调基因表达数据
[1] 596 161
DOWN_DEG=exp6[row.names(exp6) %in% rownames(DEG[DEG$change=='DOWN',]),]
dim(DOWN_DEG) #下调基因表达数据
[1] 461 161
好啦，就到这里，欢迎点进来的同学交流和学习！
本文含有隐藏内容，请 开通VIP 后查看