PCA图显示分组无差异，怎么办？

rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)#在调用as.data.frame的时，将stringsAsFactors设置为FALSE可以避免character类型自动转化为factor类型
# 注意查看下载文件的大小，检查数据
f='GSE22633_eSet.Rdata'
library(GEOquery)
# 这个包需要注意两个配置，一般来说自动化的配置是足够的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
  gset <- getGEO('GSE22633', destdir=".",
                 AnnotGPL = F,     ## 注释文件
                 getGPL = F)       ## 平台文件
  save(gset,file=f)   ## 保存到本地
}
load('GSE22633_eSet.Rdata')  ## 载入数据
class(gset)  #查看数据类型
length(gset)  #
class(gset[[1]])
gset
# assayData: 48701 features, 63 samples

# 因为这个GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的list
a=gset[[1]] #
dat=exprs(a) #a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
dim(dat)#看一下dat这个矩阵的维度
#
dat[1:4,1:4] #查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列
boxplot(dat,las=2)
dat=dat[apply(dat,1,sd)>0,]
dat[dat<0]=1
boxplot(dat,las=2)
#经观察，表达量比较小，推测已经log了

#提取临床数据
pd=pData(a)
#因关注小白菊内酯处理细胞系cell line的变化，所以要去除cell type
#且cell line: CK系列没有处理组，所以也需要去掉
#删除部分数据
pd=pd[grepl('cell line',pd$characteristics_ch1),]
pd=pd[-(33:38),]
table(pd$characteristics_ch1)
dat=dat[,rownames(pd)]

##分组 32个
group_list=ifelse(grepl('parthenolide',pd$characteristics_ch1.2),'treat','normal')
#设置参考水平，对照在前，处理在后
group_list = factor(group_list,
                    levels = c("normal","treat"))

#ID转换 GPL6102
#http://www.bio-info-trainee.com/1399.html
#BiocManager::install('illuminaHumanv2.db')
library(illuminaHumanv2.db)
## 获取表达矩阵的行名，就是探针名称
probeset <- rownames(dat)
## 使用lookup函数，找到探针在illuminaHumanv2.db中的对应基因名称
## 如果分析别的芯片数据，把illuminaHumanv2.db更换即可
SYMBOL <-  annotate::lookUp(probeset,"illuminaHumanv2.db", "SYMBOL")
## 转换为向量
symbol = as.vector(unlist(SYMBOL))
#制作ids转换文件
ids = data.frame(probeset,symbol,stringsAsFactors = F)

colnames(ids)=c('probe_id','symbol')  
ids=ids[ids$symbol != 'NA',]

ids=ids[ids$probe_id %in%  rownames(dat),]
dat[1:4,1:4]   
dat=dat[ids$probe_id,]

ids$median=apply(dat,1,median) #ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids
ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果s
dat=dat[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名
dat[1:4,1:4]  #保留每个基因ID第一次出现的信息

save(pd,dat,group_list,file = 'step1-output.Rdata')

rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)

pca_plot = function(dddd,ggggg){
  library("FactoMineR")
  library("factoextra")
  df.pca <- PCA(t(dddd), graph = FALSE)
  fviz_pca_ind(df.pca,
               #axes = c(2,3),
               geom.ind = "point",
               col.ind = ggggg ,
               addEllipses = TRUE,
               legend.title = "Groups"
  )
}
# 下面的 step1-output.Rdata 文件
load(file = 'step1-output.Rdata')
# 每次都要检测数据
dat[1:4,1:4]
table(group_list)
#normal  treat
#    16     16
g=factor( group_list )
g
#对象顺序
g=relevel(g,'normal')
pca_plot(dat,g)
ggsave('all_samples_PCA.png')
#热图
cg=names(tail(sort(apply(dat,1,sd)),1000))#apply按行（'1'是按行取，'2'是按列取）取每一行的方差，从小到大排序，取最大的1000个
library(pheatmap)
n=t(scale(t(dat[cg,]))) # 'scale'可以对log-ratio数值进行归一化
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
#两个分组
batch=pd$characteristics_ch1
batch=factor(batch)
annotation_col=data.frame(group=as.character(group_list),
                          pair = as.character(batch))
rownames(annotation_col)=colnames(dat)
pheatmap(n,show_colnames =F,show_rownames = F,
         annotation_col=annotation_col，filename = 'top1000_sd.png')

rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load(file = 'step1-output.Rdata')
# 每次都要检测数据
dat[1:4,1:4]
table(group_list)
library(limma)
g=factor( group_list )
g
g=relevel(g,'normal')
design=model.matrix(~g)
#影响表达量的批次
batch=pd$characteristics_ch1
batch=factor(batch)
## 使用 limma 的 removeBatchEffect 函数
dat[1:4,1:4]
ex_b_limma <- removeBatchEffect(dat,
                                batch = batch,
                                design = design)
#ex_b_limma 去除批次效应后的矩阵                               
dim(ex_b_limma)
ex_b_limma[1:4,1:4]
##差异分析
#去除批次效应前差异分析
fit=lmFit(dat,design)
fit=eBayes(fit)
options(digits = 4)
deg1=topTable(fit,coef=2,adjust='BH',number = Inf)
save(deg1,file = 'deg1.Rdata')
#后
fit=lmFit(ex_b_limma,design)  
fit=eBayes(fit)
options(digits = 4)
topTable(fit,coef=2,adjust='BH')
deg2=topTable(fit,coef=2,adjust='BH',number = Inf)
pca_plot(ex_b_limma,g)
ggsave('ex_b_limma.png')
save(deg2,ex_b_limma,group_list,file = 'deg2.Rdata')

rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load(file = 'deg1.Rdata')
head(deg1)
## for heatmap
if(T){
  load(file = 'step1-output.Rdata')
  # 每次都要检测数据
  dat[1:4,1:4]
  table(group_list)
  x=deg1$logFC #deg取logFC这列并将其重新赋值给x
  names(x)=rownames(deg1) #deg取probe_id这列，并将其作为名字给x
  cg=c(names(head(sort(x),100)),#对x进行从小到大排列，取前100及后100，并取其对应的探针名，作为向量赋值给cg
       names(tail(sort(x),100)))
  library(pheatmap)
  pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #对dat按照cg取行，所得到的矩阵来画热图
  n=t(scale(t(dat[cg,])))#通过“scale”对log-ratio数值进行归一化，现在的dat是行名为探针，列名为样本名，由于scale这个函数应用在不同组数据间存在差异时，需要行名为样本，因此需要用t(dat[cg,])来转换，最后再转换回来

  n[n>2]=2
  n[n< -2]= -2
  n[1:4,1:4]
  pheatmap(n,show_colnames =F,show_rownames = F)
  ac=data.frame(group=group_list)
  rownames(ac)=colnames(n) #将ac的行名也就分组信息（是‘no TNBC’还是‘TNBC’）给到n的列名，即热图中位于上方的分组信息
  pheatmap(n,show_colnames =F,
           show_rownames = F,
           cluster_cols = T,
           annotation_col=ac,filename = 'pheatmap_top200_DEG.png') #列名注释信息为ac即分组信息


}
rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load(file = 'deg2.Rdata')
head(deg2)
## for heatmap
if(T){
  # 每次都要检测数据
  ex_b_limma[1:4,1:4]
  table(group_list)
  x=deg2$logFC #deg取logFC这列并将其重新赋值给x
  names(x)=rownames(deg2) #deg取probe_id这列，并将其作为名字给x
  cg=c(names(head(sort(x),100)),#对x进行从小到大排列，取前100及后100，并取其对应的探针名，作为向量赋值给cg
       names(tail(sort(x),100)))
  library(pheatmap)
  pheatmap(ex_b_limma[cg,],show_colnames =F,show_rownames = F) #对dat按照cg取行，所得到的矩阵来画热图
  n=t(scale(t(ex_b_limma[cg,])))#通过“scale”对log-ratio数值进行归一化，现在的dat是行名为探针，列名为样本名，由于scale这个函数应用在不同组数据间存在差异时，需要行名为样本，因此需要用t(dat[cg,])来转换，最后再转换回来

  n[n>2]=2
  n[n< -2]= -2
  n[1:4,1:4]
  pheatmap(n,show_colnames =F,show_rownames = F)
  ac=data.frame(group=group_list)
  rownames(ac)=colnames(n) #将ac的行名也就分组信息（是‘no TNBC’还是‘TNBC’）给到n的列名，即热图中位于上方的分组信息
  pheatmap(n,show_colnames =F,
           show_rownames = F,
           cluster_cols = T,
           annotation_col=ac,filename = 'heatmap_top200_DEG2.png') #列名注释信息为ac即分组信息


}