在什么情况下癌症样本与正常样本是分不开的呢？

## 魔幻操作，一键清空~
rm(list = ls())
library(AnnoProbe)
library(GEOquery)
library(ggplot2)
library(ggstatsplot)
library(patchwork)
library(reshape2)
library(stringr)
library(limma)
library(tidyverse)
getOption('timeout')
options(timeout=10000)

## 获取并且检查表达量矩阵
## ～～～gse编号需修改～～～
gse_number <- "GSE54236"
dir.create(gse_number)
setwd(gse_number)
getwd()
list.files() 

#gset <- geoChina(gse_number)
gset <- getGEO(gse_number, destdir = '.', getGPL = T)
gset[[1]]

a <- gset[[1]]
dat <- exprs(a) # a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
dim(dat)        # 看一下dat这个矩阵的维度
dat[1:4,1:4]    # 查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列
range(dat)
table(is.na(dat))
#dat[is.na(dat)] <- 0
dat <- na.omit(dat)
range(dat)
dim(dat)

## 根据生物学背景及研究目的人为分组
## 通过查看说明书知道取对象a里的临床信息用pData
## 挑选一些感兴趣的临床表型。
pd <- pData(a)
colnames(pd)
pd$title
pd$characteristics_ch1
table(pd$characteristics_ch1)

## ～～～分组信息编号需修改～～～
group_list <- ifelse(grepl('Non tumor',pd$title ,ignore.case = T ), "control","case")
group_list <- factor(group_list, levels = c("control","case"))
table(group_list)

## 查看注释平台gpl 获取芯片注释信息
a
a@annotation
gpl <- getGEO(filename=paste0(a@annotation,".soft.gz"))
class(gpl)
gpl_anno <- gpl@dataTable@table
colnames(gpl_anno)

id2name <- gpl_anno[,c("ID" ,"GENE_SYMBOL")]
colnames(id2name) <- c("ID","GENE_SYMBOL")
# 1.过滤掉空的探针
id2name <- na.omit(id2name)
id2name <- id2name[which(id2name$GENE_SYMBOL!=""), ]
# 2.过滤探针一对多
id2name <- id2name[!grepl("\\///",id2name$GENE_SYMBOL), ]
head(id2name)

# 3.多对一取均值
# 合并探针ID 与基因，表达谱对应关系
# 提取表达矩阵
dat <- dat %>% 
  as.data.frame() %>% 
  rownames_to_column("ID")

exp <- merge(id2name, dat, by.x="ID", by.y="ID")

# 多对一取均值
exp <- avereps(exp[,-c(1,2)],ID = exp$GENE_SYMBOL) %>% 
  as.data.frame()

dat <- as.matrix(exp[,pd$geo_accession])
dim(dat)
fivenum(dat['CRP',])
fivenum(dat['GAPDH',])
dat[1:5, 1:6]
save(gse_number, dat, group_list, pd, file = 'step1_output.Rdata')

## 数据检查
load( file = 'step1_output.Rdata')
source('../scripts/run_check.R')
run_check(dat, group_list, target_gene='AFP', pro=gse_number,path='./')

## 人类看家基因列表如下（部分基因，最常用的是CRP和CRP）
source('../scripts/run_DEG.R')
group_list
deg <- run_DEG(dat, group_list, target_gene=c('AFP'),pro=gse_number,path='./')
head(deg)
table(deg$g)

exp <- dat[deg$g!="stable",]
## 下面是画PCA的必须操作，需要看说明书。
exp=t(exp)#画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转换
exp=as.data.frame(exp)#将matrix转换为data.frame 

library("FactoMineR")#画主成分分析图需要加载这两个包
library("factoextra")  
#～～～主成分分析图p2～～～
dat.pca <- PCA(exp , graph = FALSE)#现在exp最后一列是group_list，需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的
pro <- gse_number
this_title <- paste0(pro,'_PCA')
p2 <- fviz_pca_ind(dat.pca,
                 geom.ind = "point", # show points only (nbut not "text")
                 col.ind = group_list, # color by groups
                 palette = "Dark2",
                 addEllipses = TRUE, # Concentration ellipses
                 legend.title = "Groups")+
ggtitle(this_title) +
theme_ggstatsplot() +
theme(plot.title = element_text(size=12,hjust = 0.5))

p2

cg <- names(tail(sort(apply(dat,1,sd)),1000))#apply按行（'1'是按行取，'2'是按列取）取每一行的方差，从小到大排序，取最大的1000个

exp <- dat[cg,]
## 下面是画PCA的必须操作，需要看说明书。
exp=t(exp)#画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转换
exp=as.data.frame(exp)#将matrix转换为data.frame 

library("FactoMineR")#画主成分分析图需要加载这两个包
library("factoextra")  
#～～～主成分分析图p2～～～
dat.pca <- PCA(exp , graph = FALSE)#现在exp最后一列是group_list，需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的
pro <- gse_number
this_title <- paste0(pro,'_PCA')
p2 <- fviz_pca_ind(dat.pca,
                 geom.ind = "point", # show points only (nbut not "text")
                 col.ind = group_list, # color by groups
                 palette = "Dark2",
                 addEllipses = TRUE, # Concentration ellipses
                 legend.title = "Groups")+
ggtitle(this_title) +
theme_ggstatsplot() +
theme(plot.title = element_text(size=12,hjust = 0.5))

p2