参考：http://www.bio-info-trainee.com/3754.html
下载GSE76275数据地址https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi
1、下载表达矩阵
#step1
rm(list=ls())
library(GEOmirror)
gset <- geoChina('GSE76275')
gset=gset[[1]]
exp <- exprs(gset) #表达矩阵
dim(exp)#265个样本，54675个探针
exp[1:4,1:4]
pdata <- pData(gset) #分组信息
#通过观察，找到分组信息列
table(pdata$characteristics_ch1.1)
pdata$characteristics_ch1.1=='triple-negative status: not TN'
#利用ifelse语句生成二分类分组信息(是否为三阴性乳腺癌)
group_list <- ifelse(pdata$characteristics_ch1.1=='triple-negative status: not TN',
'noTNBC', 'TNBC')
save(exp, group_list, file = "exp_group.Rdata")
2、检查数据
2.1PCA分析
step2 检查
##PCA
## 下面是画PCA的必须操作，需要看说明书。
install.packages(c('FactoMineR','factoextra'))
library(FactoMineR)#画主成分分析图需要加载这两个包
library(factoextra)
exp1 <- t(exp) #画PCA图时要求是行名是样本名，列名时探针名，因此此时需要转换
exp1 <- as.data.frame(exp1)#将matrix转换为data.frame
exp1[1:4,1:4]
exp1 <- cbind(exp1,group_list) #cbind横向追加，即将分组信息追加到最后一列
exp1.pca <- PCA(exp1[,-54676],graph=FALSE)
exp1.pca <- PCA(exp1[,-col(exp1],graph=FALSE)#和上面方式相同
fviz_pca_ind(exp1.pca,
geom.ind = "point", # show points only (nbut not "text")
col.ind = exp1$group_list, # color by groups
palette = c("#00AFBB", "#E7B800"),
addEllipses = TRUE, # Concentration ellipses
legend.title = "Groups"
)
pca.png
可以看到那些不属于TNBC的样本跟TNBC组有着很清晰的界限，而且可以看到TNBC组本身也有着界限分明的两列。
2.2 heatmap
##热图
rm(list=ls())
load('exp_group.Rdata')
Sys.setenv(LANGUAGE='EN')
exp[1:4,1:4]
system.time(apply(exp,1,sd))
system.time(for (i in 1:nrow(exp)){
sd(exp[1,])
})
tail(sort(apply(exp,1,sd)),1000)
#取出方差最大的1000个基因
cg <- names(tail(sort(apply(exp,1,sd)),1000))
#apply按行（'1'是按行取，'2'是按列取）取每一行的方差，从小到大排序，取最大的1000个
#sort 排序，默认升序排序。tail从后向前显示，默认后6行
library(pheatmap)
pheatmap(exp[cg,],show_colnames = F,
show_rownames =F)
##将表达量进行归一化
n=t(scale(t(exp[cg,])))
#通过“scale”对log-ratio数值进行归一化，现在的dat是行名为探针，列名为样本名，
#由于scale这个函数应用在不同组数据间存在差异时，需要行名为样本，因此需要用t(dat[cg,])来转换，最后再转换回来
#t 对矩阵进行行列转换
n[n>2]=2 #限定上限，使表达量大于2的等于2
n[n< -2]= -2 #限定下限，使表达量小于-2的等于-2
pheatmap(n,show_colnames = F,
show_rownames =F)
#标上分组信息
group_exp <- data.frame(group_list) #把character向量转换为数据框
head(group_exp)
rownames(group_exp) <- colnames(n)
pheatmap(n,show_colnames = F,
show_rownames =F,
annotation_col = group_exp)
简化方式
rm(list=ls())
load('exp_group.Rdata')
Sys.setenv(LANGUAGE='EN')
tail(sort(apply(exp,1,sd)),1000)
#取出方差最大的1000个基因
cg <- names(tail(sort(apply(exp,1,sd)),1000))
cm <- exp[cg,] ##choose_matrix
cm[1:4,1:4]
cm <- t(scale(t(cm)))#标准化处理
cm[1:4,1:4]
cm[cm>2]=2 #限定上限，使表达量大于2的等于2
cm[cm< -2]= -2 #限定下限，使表达量小于-2的等于-2
cm[cm>2]=2;cm[cm< -2]= -2
#热图的特定分组格式
group_dat <- data.frame(group=group_list,row.names = colnames(exp))
head(group_dat)
library(pheatmap)
pheatmap(cm,show_rownames = F,
show_colnames = F,
annotation_col = group_dat)
#补充
rm(list=ls())
cg <- names(tail(sort(apply(exp,1,sd)),1000))
cm <- exp[cg,]
cm <- t(scale(t(cm)))
cm[cm>2 ]=2
cm[cm< -2]=-2
group_exp <- data.frame(group=group_list,row.names= colnames(exp))
head(group_exp)
library(pheatmap)
pheatmap(cm,show_rownames = F,
show_colnames = F,
annotation_col = group_exp)
## 单独看 TNBC 的 top 1000 基因的热图
tnbc_exp <- exp[,group_list=='TNBC']
cg1 <- names(tail(sort(apply(tnbc_exp,1,sd)),1000))
cm1 <- tnbc_exp[cg1,]
cm1 <- t(scale(t(cm1)))
cm1[cm1>2]=2
cm1[cm1< -2]= -2
cm1[1:4,1:4]
pheatmap(cm1,show_rownames = F,
show_colnames = F)
## 单独看 noTNBC TNBC 的 top 1000 基因的热图
notnbc_exp <- exp[,group_list=='noTNBC']
cg2 <- names(tail(sort(apply(notnbc_exp,1,sd)),1000))
cm2 <- notnbc_exp[cg2,]
cm2 <- t(scale(t(cm2)))
cm2[cm2>2]=2
cm2[cm2< -2]= -2
pheatmap(cm2,show_rownames = F,
show_colnames = F)
heatmap
TNBC.heatmap
noTNBC.heatmap
网友评论