Learn R GEO

>rm(list = ls())
>library(GEOquery)
#先去网页确定是否是表达芯片数据，不是的话不能用本流程。
>gse_number = "GSE56649"
>eSet <- getGEO(gse_number, destdir = '.', getGPL = F) #getGEO下载并读取文件
> class(eSet)
[1] "list"
> length(eSet)
[1] 1
> eSet = eSet[[1]]
> class(eSet)
[1] "ExpressionSet"  #哪种数据类型
attr(,"package")      # 属性之包的名字
[1] "Biobase"          #Biobase包规定的格式"ExpressionSet"

> exp <- exprs(eSet)
> dim(exp)
[1] 54675    22
> exp[1:4,1:4]  #检查矩阵是否正常，如果是空的就会报错，空的和有负值的、有异常值的矩阵需要处理原始数据。是否取过loglog 是否有负值
          GSM1366348 GSM1366349 GSM1366350 GSM1366351
1007_s_at    279.156    202.866    406.818    232.668
1053_at      487.700    409.018    393.810    397.127
117_at       666.869    388.589    704.633    953.481
121_at       240.646    361.198    305.229    374.757
> #如果表达矩阵为空，大多数是转录组数据，不能用这个流程(后面另讲)。
> #自行判断是否需要log 可以画boxplot()
> exp = log2(exp+1) 正常log在0-20之间
> boxplot(exp)

### (2)提取临床信息
>pd <- pData(eSet)
#(3)让exp列名与pd的行名顺序完全一致
>p = identical(rownames(pd),colnames(exp));p#判断信息是否以一对应
>if(!p) exp = exp[,match(rownames(pd),colnames(exp))]
#分组信息来自临床信息，分组信息需要与表达矩阵列名一一对应
#临床信息需要与表达矩阵一一对应

>gpl_number <- eSet@annotation;gpl_number #提取子集 注意什么时候用@，什么时候用$，看图1
[1] "GPL570"
>save(gse_number,pd,exp,gpl_number,file = "step1output.Rdata")
#保存 gse_number(原本的编号),pd(临床信息),exp(表达矩阵),gpl_number(芯片编号)

# 从临床样本中获得实验分组（在表格中慢慢找，代码如何实现看下）
rm(list = ls())  
load(file = "step1output.Rdata")
library(stringr)
# 标准流程代码是二分组，多分组数据的分析后面另讲
# 生成Group向量的三种常规方法，三选一，选谁就把第几个逻辑值写成T，另外两个为F。如果三种办法都不适用，可以继续往后写else if
# 1.Group----
# 第一种方法，有现成的可以用来分组的列
  Group = pd$`disease state:ch1`  #pd$`cell type:ch1`注意这个``不能删
}else if(F){
# 第二种方法，自己生成，但需要看清哪是RA哪是control
  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")
#felse(str_detect()) 两个函数连用是用来分组的神器
#str_detect(pd$source_name_ch1,"control"）-source_name_ch1这一列是否含有control
#ifelse(str_detect(pd$source_name_ch1,"control"),"control","RA")
是control就填control不是填RA
#三种代码完整方法
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，对照组在前，处理组在后 ；因子正文与levels不对应时会产生NA
Group = factor(Group,levels = c("control","RA"))
Group
#Group是一个有重复值的向量 是分类型数据，适合用因子的形式
#factor直接转换并自动生成levels (control和RA)，顺序以字母排序为准
#levels顺序有意义，在第一个位置的水平是参考水平
#参考水平将在做差异分析时，被设为对照组
#所以需要控制levels的顺序
#levels = c("control","RA") 写了按照写的顺序，control位参考水平

>library(tinyarray) #tinyarry这个包是生信技能树写的，find_anno()是其中的一个函数
>find_anno(gpl_number) #打出找注释的代码
[1] "`library(hgu133plus2.db);ids <- toTable(hgu133plus2SYMBOL)` (方法一，R包）and `ids <- AnnoProbe::idmap('GPL570')` are both avaliable"（方法二，idmap，idmap包含这个平台）
> ids <- AnnoProbe::idmap('GPL570') 
trying URL 'http://49.235.27.111/GEOmirror/GPL/GPL570_bioc.rda'
Content type 'application/octet-stream' length 264292 bytes (258 KB)
==================================================
downloaded 258 KB
file downloaded in /Users/zhuo/learn /pipeline

>gpl_number 
#http://WWW.bio-info-trainee.com/1399.html 网站（看图）上找GPL对应的R包编号，
>if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db") #下载R包
>library(hgu133plus2.db)
>ls("package:hgu133plus2.db") #只是展示包里有什么，但只要SYMBOL这个
>ids <- toTable(hgu133plus2SYMBOL) #ids <- toTable(）提取 不同的包进行替换（看图）
>head(ids) #看到所需要的结果

if(F){
  #注：表格读取参数、文件列名不统一(如Gene Symbol 变成了Gene_Symbol)，活学活用，有的表格里没有symbol列，也有的GPL平台没有提供注释表格
  b = read.delim("GPL570-55999.txt",
                 check.names = F,
                 comment.char = "#")
  colnames(b)  #看图colnames(b)
  ids2 = b[,c("ID","Gene Symbol")] #提取出来想要的两列 
  colnames(ids2) = c("probe_id","symbol") #统一重新命名，后面就不用改代码了 看图ids2 ；symbol中有的名称有/// 前后有两个名称，代表是两个基因-非特异性探针，直接去除
  k1 = ids2$symbol!="";table(k1) #检测 空的 看图K1
  k2 = !str_detect(ids2$symbol,"///");table(k2) #检测非特异性探针
  ids2 = ids2[ k1 & k2,] #取既不是空的也不是非特异性的 & 同时满足两个条件
  # ids = ids2 ids2是为了教课时好区分 后面用的时候把前面的#去掉，省改后面的代码
}
#若用第二种方法把最前面的if(F)变成T

#仿制的前四列
dat=as.data.frame(t(exp)) #t() 转置 as.data.frame()作为数据框
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") #保存pca的数据

cg=names(tail(sort(apply(exp,1,sd)),1000)) # 挑选出1000个方差最大的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() -3，3（不同的结果设置的色带分配值不一样）是设置色带分布范围 分配颜色色带分配100种颜色
         ) 

library(limma) #limma 包
design=model.matrix(~Group) #构建模型矩阵
fit=lmFit(exp,design) #做线性拟合
fit=eBayes(fit) #贝叶斯检验
deg=topTable(fit,coef=2,number = Inf) #topTable 提取结果

#简化思维 把上面的函数写成一个自己新的函数如egp()
egp=function(exp,Group){
  design=model.matrix(~Group)
  fit=lmFit(exp,design)
  fit=eBayes(fit)
  deg=topTable(fit,coef=2,number = Inf)
  return(deg)
}
a = egp(exp,Group)
## 鉴定两个函数结果是否一致
identical(a,deg)
#这一步结束得到的是deg（六列数据，还需4列，看图差异分析后的数据整理）
#为deg数据框添加几列
  #1.加probe_id列，把行名变成一列
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))
  #2.加上探针注释
ids = ids[!duplicated(ids$symbol),]  # 这个代码是随机去重的方式
ids =distinct(ids,symbol,.keep_all = T)#这个代码也是随机去重的方式
###出现多个探针对应一个基因的情况，所以需对基因进行去重
####方法1：随机去重 
####方法2：保留行和/行平均值最大的探针 
####方法3：取多个探针的平均值
#其他去重方式在 “zz.去重方式.R”这个文件里
deg <- inner_join(deg,ids,by="probe_id") #inner_join 取交集
nrow(deg)
  #3.加change列,标记上下调基因
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"))) #mutate加上updown...
table(deg$change)
  #4.加ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join）
library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(deg$symbol,  #bitr() 数据类型转换
            fromType = "SYMBOL",  #从..ID转换为
            toType = "ENTREZID",   # ...ID
            OrgDb = org.Hs.eg.db)  #转换一句的数据库-人类
#转换完结果看图b
#其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL")) #两个表的样本名称不一样，用inner_join() 后面要写by=c(""="")
save(Group,deg,logFC_t,p_t,gse_number,file = "step4output.Rdata") #保存
#到此有了10列 用于后续做图

library(dplyr)
library(ggplot2)
dat  = distinct(deg,symbol,.keep_all = T)

p <- ggplot(data = dat, 
            aes(x = logFC,       #横坐标
                y = -log10(P.Value))) + #纵坐标
  geom_point(alpha=0.4, size=3.5,  #画点图 alpha透明度 size点的大小
             aes(color=change)) + #这步之前是做了普通的点图
  scale_color_manual(values=c("blue", "grey","red"))+ #设定颜色
  geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",linewidth=0.8) +.  #geom_vline()添加参考线（竖线）
  geom_hline(yintercept = -log10(p_t),lty=4,col="black",linewidth=0.8) +  #geom_hline()添加参考线（横线）
  theme_bw() #去掉灰色的背景
p
#在火山图上标记感兴趣的基因
for_label <- dat%>% 
  filter(symbol %in% c("HADHA","LRRFIP1")) #挑选出想要的基因所在的行

volcano_plot <- p +     #叠加涂层
  geom_point(size = 3, shape = 1, data = for_label) + #画了空心的圆
  ggrepel::geom_label_repel(   #添加标签 （基因的名+边框）
    aes(label = symbol),
    data = for_label,  #数据是来自新建的（挑选的基因所在的行，for_label)
    color="black"
  )
volcano_plot #结果图看图c

load(file = 'step2output.Rdata')
# 表达矩阵行名替换
exp = exp[dat$probe_id,] #转变为以基因为行名的
rownames(exp) = dat$symbol
if(T){
  #取前10上调和前10下调 （可按logFC取也可按P value取）
  library(dplyr)
  dat2 = dat %>%
    filter(change!="stable") %>%  #筛选基因change不是stable的-找差异的
    arrange(logFC)                     #按照logFC从小到大排序
  cg = c(head(dat2$symbol,10),      #取前十后十
         tail(dat2$symbol,10))
}else{
#全部差异基因
  cg = dat$symbol[dat$change !="stable"]
  length(cg)
}
n=exp[cg,] #统计就会出现 20个差异表达基因和22个样本
dim(n)
#差异基因热图
library(pheatmap)
annotation_col=data.frame(group=Group)
rownames(annotation_col)=colnames(n) 
heatmap_plot <- pheatmap(n,show_colnames =F,
                         scale = "row",
                         #cluster_cols = F,  #如果基因数太多，就可以将基因名删掉
                         annotation_col=annotation_col,
                         breaks = seq(-3,3,length.out = 100)
) 
heatmap_plot
# 拼图
library(patchwork)
volcano_plot +heatmap_plot$gtable #拼火山图和热图

library(corrplot)
g = sample(deg$symbol[1:500],10) # 这里是随机取样，注意换成自己感兴趣的基因
g

M = cor(t(exp[g,])) #cor()用于计算基因的相关性，提供矩阵数据，计算列于列之间的相关性，看图
pheatmap(M)

# 配色R包
library(paletteer)
my_color = rev(paletteer_d("RColorBrewer::RdYlBu")) #paletteer_d("RColorBrewer::RdYlBu") 在控制台上展示颜色，原本是蓝色-正相关，红色-负相关，rev将它转换一下 看图c
my_color = colorRampPalette(my_color)(10) #只取10种颜色
corrplot(M, type="upper",
         method="pie",  #以饼图代表相关系数大小
         order="hclust", 
         col=my_color,   #以颜色代表正负相关
         tl.col="black", 
         tl.srt=45)
library(cowplot)
cor_plot <- recordPlot() #将画板上的图记录下来，为了后续可以拼图
# 拼图
load("pca_plot.Rdata")
plot_grid(pca_plot,cor_plot,
          volcano_plot,heatmap_plot$gtable)
dev.off()
#保存
pdf("deg.pdf")
plot_grid(pca_plot,cor_plot,
          volcano_plot,heatmap_plot$gtable)
dev.off()

a  = deg$symbol[1]
boxplot(exp[a,]~Group)
deg$logFC[1]