跟小洁老师学习GEO的第二天

#数据下载
rm(list = ls())
library(GEOquery)
#先去网页确定是否是表达芯片数据，不是的话不能用本流程。
gse_number = "GSE28345"
library(AnnoProbe)
eSet <- geoChina(gse_number, destdir = '.')
class(eSet)
length(eSet)
eSet = eSet[[1]]

options("repos"="https://mirrors.ustc.edu.cn/CRAN/")
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")

cran_packages <- c('tidyr',
                   'tibble',
                   'dplyr',
                   'stringr',
                   'ggplot2',
                   'ggpubr',
                   'factoextra',
                   'FactoMineR',
                   'devtools',
                   'cowplot',
                   'patchwork',
                   'basetheme',
                   'paletteer',
                   'AnnoProbe',
                   'ggthemes',
                   'VennDiagram',
                   'tinyarray') 
Biocductor_packages <- c('GEOquery',
                         'hgu133plus2.db',
                         'ggnewscale',
                         "limma",
                         "impute",
                         "GSEABase",
                         "GSVA",
                         "clusterProfiler",
                         "org.Hs.eg.db",
                         "preprocessCore",
                         "enrichplot")

for (pkg in cran_packages){
  if (! require(pkg,character.only=T) ) {
    install.packages(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}


for (pkg in Biocductor_packages){
  if (! require(pkg,character.only=T) ) {
    BiocManager::install(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}

#前面的所有提示和报错都先不要管。主要看这里
for (pkg in c(Biocductor_packages,cran_packages)){
  require(pkg,character.only=T) 
}
#没有任何提示就是成功了，如果有warning xx包不存在，用library检查一下。

#library报错，就单独安装。

rm(list = ls())  
load(file = "step1output.Rdata")
library(stringr)
# 标准流程代码是二分组，多分组数据的分析后面另讲
# 生成Group向量的三种常规方法，三选一，选谁就把第几个逻辑值写成T，另外两个为F。如果三种办法都不适用，可以继续往后写else if
if(F){
  # 1.Group----
  # 第一种方法，有现成的可以用来分组的列
  Group = pd$`disease state:ch1` 
}else if(F){
  # 第二种方法，自己生成
  Group = c(rep("RA",times=13),
            rep("control",times=9))
  Group = rep(c("RA","control"),times = c(13,9))
}else if(T){
  # 第三种方法，使用字符串处理的函数获取分组
  Group=ifelse(str_detect(pd$source_name_ch1,"control"),
               "control",
               "RA")
}

# 需要把Group转换成因子，并设置参考水平，指定levels，对照组在前，处理组在后
Group = factor(Group,levels = c("control","RA"))
Group

#捷径
library(tinyarray)
find_anno(gpl_number) #打出找注释的代码
ids <- AnnoProbe::idmap('GPL570')

gpl_number 
#http://www.bio-info-trainee.com/1399.html
if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db")
library(hgu133plus2.db)
ls("package:hgu133plus2.db")
ids <- toTable(hgu133plus2SYMBOL)
head(ids)

if(F){
  #注：表格读取参数、文件列名不统一，活学活用，有的表格里没有symbol列，也有的GPL平台没有提供注释表格
  b = read.delim("GPL570-55999.txt",
                 check.names = F,
                 comment.char = "#")
  colnames(b)
  ids2 = b[,c("ID","Gene Symbol")]
  colnames(ids2) = c("probe_id","symbol")
  k1 = ids2$symbol!="";table(k1)
  k2 = !str_detect(ids2$symbol,"///");table(k2)
  ids2 = ids2[ k1 & k2,]
  # ids = ids2
}

save(exp,Group,ids,gse_number,file = "step2output.Rdata")

rm(list = ls())  
load(file = "step1output.Rdata")
load(file = "step2output.Rdata")
#输入数据：exp和Group
#Principal Component Analysis
#http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials

dat=as.data.frame(t(exp))
library(FactoMineR)
library(factoextra) 
dat.pca <- PCA(dat, graph = FALSE)
pca_plot <- fviz_pca_ind(dat.pca,
                         geom.ind = "point", # show points only (nbut not "text")
                         col.ind = Group, # color by groups
                         palette = c("#00AFBB", "#E7B800"),
                         addEllipses = TRUE, # Concentration ellipses
                         legend.title = "Groups"
)
pca_plot
save(pca_plot,file = "pca_plot.Rdata")

cg=names(tail(sort(apply(exp,1,sd)),1000))
n=exp[cg,]#把方差最大的基因挑选出来

# 直接画热图，对比不鲜明
library(pheatmap)
annotation_col=data.frame(group=Group)
rownames(annotation_col)=colnames(n) 
pheatmap(n,
         show_colnames =F,
         show_rownames = F,
         annotation_col=annotation_col
)

# 按行标准化，放大了行内部的差别
pheatmap(n,
         show_colnames =F,
         show_rownames = F,
         annotation_col=annotation_col,
         scale = "row",
         breaks = seq(-3,3,length.out = 100)
         ) #breaks参数：设置色带分配范围，100种数字就是100种颜色
dev.off()

rm(list = ls()) 
load(file = "step2output.Rdata")
#差异分析，用limma包来做
#需要表达矩阵和Group，不需要改
library(limma)
design=model.matrix(~Group)
fit=lmFit(exp,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)#提取结果

library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))

ids = ids[!duplicated(ids$symbol),]

#其他去重方式
rm(list = ls())
load("step2output.Rdata")
#1.保留最大值
exp2 = exp[ids$probe_id,]
identical(ids$probe_id,rownames(exp2))
ids = ids[order(rowSums(exp2),decreasing = T),]
ids = ids[!duplicated(ids$symbol),];nrow(ids)
# 拿这个ids去inner_join
#2.求平均值
rm(list = ls())
load("step2output.Rdata")
exp3 = exp[ids$probe_id,]
rownames(exp3) = ids$symbol
exp3[1:4,1:4]
exp4 = limma::avereps(exp3)
# 此时拿到的exp4已经是一个基因为行名的表达矩阵，直接差异分析，不再需要inner_join

deg <- inner_join(deg,ids,by="probe_id")
nrow(deg)

logFC_t=1
p_t = 0.05
k1 = (deg$P.Value < p_t)&(deg$logFC < -logFC_t)
k2 = (deg$P.Value < p_t)&(deg$logFC > logFC_t)
deg <- mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable")))
table(deg$change)

library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(deg$symbol, 
            fromType = "SYMBOL",
            toType = "ENTREZID",
            OrgDb = org.Hs.eg.db)#人类
#其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
save(Group,deg,logFC_t,p_t,gse_number,file = "step4output.Rdata")