生信入门DAY7-9

rm(list = ls())
#打破下载时间的限制,改前60秒，改后10w秒
options(timeout = 100000) 
options(scipen = 20)#不要以科学计数法表示

#传统下载方式
library(GEOquery)
eSet = getGEO("GSE7305", destdir = '.', getGPL = F)
#上面的函数做两件事情，一是下载二是读取，外源放进来的文件仍然要运行这句代码，后面两个参数先不用管
#不用网络直接读取：getGEO(file="GSE7305_series_matrix.txt.gz")
#网速太慢，下不下来怎么办
#1.从网页上下载/发链接让别人帮忙下，放在工作目录里，但上面那句代码还是得运行
#2.试试geoChina,只能下载2019年前的表达芯片数据，曾老师写的包
#library(AnnoProbe)
#eSet = geoChina("GSE7305") #选择性代替第8行
#研究一下这个eSet
class(eSet)
length(eSet)

eSet = eSet[[1]] 
class(eSet)
?ExpressionSet#侠义的对象，包括哪些由GSE作者设定

#(1)提取表达矩阵exp，只是芯片能用，其它组学用不了
exp <- exprs(eSet)
#⭐第一个要检查的地方👇，表达矩阵行列数，正常是几万行，列数=样本数，
#行代表探针数，如果0行说明不是表达芯片或者是遇到特殊情况，不能用此流程分析
dim(exp)
#⭐二个要检查的地方👇
range(exp)#看数据范围决定是否需要log，是否有负值，异常值，如有负值，结合箱线图进一步判断
#⭐可能要修改的地方👇
exp = log2(exp+1) #需要log才log，不需要log要注释掉这一句，+1是为了避免负值
#取过log的数值范围在0-20之间，取两次log的话数值在0-4之间。
#⭐第三个要检查的地方👇
boxplot(exp,las = 2) #看是否有异常样本
#看样本齐性，也就是中位线和上下四分位线基本在一条线上；看纵坐标范围是否在0-20之间。
#减少样本数量的代码，样本数量很多的时候
#boxplot(exp,las=2)
#boxplot(exp[,seq(1,ncol(exp),4)],las=2)
#随机取样的方式画图，样本数量很多的时候
#boxplot(exp[,sample(1:ncol(exp),7)],las=2)
#异常样本怎么处理：①删掉异常样本；②exp=limma:normlizeBetweenArrays(exp)

#(2)提取临床信息（找分组），其它组学也可以用
pd <- pData(eSet)#所有组学都有，一级一级去取eSet@phenoData@data
#⭐多分组中提取两分组的代码示例，二分组不需要
library(stringr)
if(F){
  #因为现在这个例子不是多分组，所以编造一列做示例。
  pd$fake = paste0(rep(c("a","b","c","d"),each = 5),1:5)
  k1 = str_detect(pd$fake,"b");table(k1)
  k2 = str_detect(pd$fake,"c");table(k2)
  pd = pd[k1|k2,]
}
#(3)让exp列名与pd的行名顺序完全一致
#函数identical，判断是否一致，包括数据类型和数据结构
p = identical(rownames(pd),colnames(exp));p
if(!p) {
  s = intersect(rownames(pd),colnames(exp))
  exp = exp[,s]
  pd = pd[s,]
}

#(4)提取芯片平台编号，后面要根据它来找探针注释
gpl_number <- eSet@annotation;gpl_number
save(pd,exp,gpl_number,file = "step1output.Rdata")

# 原始数据处理的代码，按需学习
# https://mp.weixin.qq.com/s/0g8XkhXM3PndtPd-BUiVgw

#当箱线图样本太多的解决办法：
boxplot(exp,las = 2) #看是否有异常样本
boxplot(exp[,seq(1,ncol(exp),4)],las = 2)
students = c("曾颖","王小冬","武文豪","弘毅","青柠","燃烧小宇宙","叶茂艺")
sample(students,10000,replace = T) |> table() |> sort(decreasing = T)
boxplot(exp[,sample(1:ncol(exp),10)],las = 2)
#管道符号：|>

# Group(实验分组)和ids(探针注释)
rm(list = ls())  
load(file = "step1output.Rdata")
# 1.Group----
library(stringr)
# 标准流程代码是二分组，多分组数据的分析后面另讲
#⭐要修改的地方：分组信息,必须学会ifelse和str_detect
k = str_detect(pd$title,"Normal");table(k) #不在title就在pd的其他列
Group = ifelse(k,"Normal","Disease")

# 需要把Group转换成因子，并设置参考水平，指定levels
#⭐要修改的地方，对照组在前，处理组在后。
#如果不人为设定levels，直接factor(Group)，那么factor会自定义一个水平。
Group = factor(Group,levels = c("Normal","Disease"))
Group

 #⭐检查自己得到的分组是否正确
data.frame(pd$title,Group)
#2.探针注释的获取-----------------
#四种方法，方法1里找不到就从方法2找，以此类推。
#方法1 BioconductorR包(最常用)
#⭐要操作的地方
library(tinyarray)
gpl_number #首先看看编号是多少
#View(pkg_all) 
#然后在pkg_all里搜索gpl编号，找到对应的R包前缀(第二列)，没搜到就是没有R包,再看方法2。
#也可以用代码直接得到对应的R包前缀：
pkg_all[pkg_all$gpl==gpl_number,2]
#如果输出character(0)则表示pkg_all里没有这个gpl编号，表示没有结果。
#⭐要操作的地方
#用上面找到的前缀替换下面所有的hgu133plus2，共5处
if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db",ask = F,update = F)
library(hgu133plus2.db)
ls("package:hgu133plus2.db") #列出R包里都有啥
ids <- toTable(hgu133plus2SYMBOL) #把R包里的注释表格变成数据框
# 方法2 下载并读取GPL网页的表格文件，按列取子集
#⭐要操作的地方
library(tinyarray)
get_gpl_txt(gpl_number) #获取表格文件的下载链接
# 接下来是复制网址去浏览器下载、放在工作目录下、读取、提取探针id和基因symbol(没有现成的需要拆分和转换)，不同文件代码不统一，等看同学们的例子。
# 注意:最终的数据ids只能有两列，第一列列名是probe_id,第二列列名是symbol,且都是字符型，否则后面代码要报错咯。
# 方法3 官网下载注释文件并读取
# 方法4 自主注释，了解一下
#https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA
save(exp,Group,ids,file = "step2output.Rdata")

#比较复杂的探针注释参考资料
#资料1：拆分取列https://www.yuque.com/xiaojiewanglezenmofenshen/kzgwzl/sv262capcgg9o8s5
#资料2：多种id的转换https://www.yuque.com/xiaojiewanglezenmofenshen/kzgwzl/pn0s1mmsaxocfynb?singleDoc# 《又一个有点难的探针注释（多种id的转换）》

rm(list = ls())  
#⭐不用改代码，会看结果就行
load(file = "step2output.Rdata")
#输入数据：exp和Group
#Principal Component Analysis
#下面这个链接有画PCA图的链接
#http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials

# 1.PCA 图（只用到了exp和group）----
#因为别人的原始数据是以样本作为行，而我们的样本放在列，所以需要转置。
class(exp)
dat=as.data.frame(t(exp))#函数t实现了矩阵的行列互换，并转换成数据框
library(FactoMineR)
library(factoextra) 
dat.pca <- PCA(dat, graph = FALSE)
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"
)

#如果少于4个点就不会画出圈，是正常的。因为圈是置信区间，样本太少无法计算，不是必须的。
# 2.top 1000 sd 热图，是表达量最离散的1000个基因，包括了组间的---- 
#因为整个表达矩阵画不了热图，我们挑出一部分来画
g = names(tail(sort(apply(exp,1,sd)),1000)) #day7-apply的思考题
n = exp[g,]
library(pheatmap)
annotation_col = data.frame(row.names = colnames(n),
                            Group = Group)
pheatmap(n,
         show_colnames =F,
         show_rownames = F,
         #color = colorRampPalette(c("navy", "white", "firebrick3"))(100),
         annotation_col=annotation_col,#添加注释条用的
         scale = "row", #按行标准化，现在的行是GSM，只保留行内差别，不保留行间差别，会把数据范围缩放到大概-5~5之间
         breaks = seq(-3,3,length.out = 100) #设置色带分布范围为-3~3之间，使中间色带加深颜色，超出此范围的数字显示极限颜色，100个刻度。
         ) 
#画图代码不报错，但不出图用dev.off解决
#关于breaks的进一步学习
#?pheatmap::pheatmap
#library(RColorBrewer)
#colorRampPalette(rev(brewer.pal(n = 7, name =
#                                  "RdYlBu")))(100)

rm(list = ls()) 
load(file = "step2output.Rdata")
#差异分析
library(limma)
design = model.matrix(~Group)
fit = lmFit(exp,design)#线性拟合
fit = eBayes(fit)#计算P值
deg = topTable(fit,coef = 2,number = Inf)
#提取差异分析结果的函数toptable,这个表格就是我们差异分析的结果，inf是正无穷，有多少行显示多少行。

#为deg数据框添加几列
#1.加probe_id列，把行名变成一列
library(dplyr)
deg = mutate(deg,probe_id = rownames(deg))
#2.加上探针注释
ids = distinct(ids,symbol,.keep_all = T)
#多个探针对应一个基因，把ids进行去重，只保留该基因出现的第一个探针
#其他去重方式在zz.去重方式.R
deg = inner_join(deg,ids,by="probe_id")
#验证一下ids里面的20815个探针是否全部出自于deg里面
#table(ids$probe_id %in% deg$probe_id)
#⭐检查
nrow(deg) #如果行数为0就是你找的探针注释是错的。

#3.加change列,标记上下调基因，是做过的练习题
#⭐阈值，可按需修改
logFC_t = 1
p_t = 0.05
#⭐思考，如何使用padj而非p值
#colnames(deg) 哪里出现了P.Value，就把哪改成adj.P.Val
#adj.P.Val为了避免假阳性的存在，会把所有的P值变大
#table(deg$adj.P.Val>deg$P.Value)
#which(deg$adj.P.Val<=deg$P.Value)
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(ggplot2)
ggplot(data = deg, aes(x = logFC, y = -log10(P.Value))) +
  geom_point(alpha=0.4, size=3.5, 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_hline(yintercept = -log10(p_t),lty=4,col="black",linewidth=0.8) +
  theme_bw()

# 差异基因热图----
# 表达矩阵行名替换为基因名
exp = exp[deg$probe_id,]#把表达矩阵中的探针也按要求筛选
#dim(exp)
#identical(rownames(exp),deg$probe_id)
#运行顺序错乱：subscript out of bounds，超纲的意思
rownames(exp) = deg$symbol#探针表达矩阵换成了基因表达矩阵
diff_gene = deg$symbol[deg$change !="stable"]
n = exp[diff_gene,]
library(pheatmap)
annotation_col = data.frame(group = Group)
rownames(annotation_col) = colnames(n) 
pheatmap(n,show_colnames =F,
         show_rownames = F,
         scale = "row",
         #cluster_cols = F, 如果聚类树和分组情况不一致，第一种办法是调整阈值logFC，第二种办法是用这句代码可以让聚类树消失。
         annotation_col=annotation_col,
         breaks = seq(-3,3,length.out = 100)
) 

#4.加ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join）
library(clusterProfiler)
library(org.Hs.eg.db)
s2e = bitr(deg$symbol, 
           fromType = "SYMBOL",
           toType = "ENTREZID",
           OrgDb = org.Hs.eg.db)#⭐Hs是人类,注意物种
#一部分基因没匹配上是正常的。只要物种正确，<30%的失败都没事。
#其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
#构建自己研究的物种Orgdb，在刘小泽老师的简书上有教程。
#⭐检查行数，减少太多说明不正常
nrow(deg) 
deg = inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
#多了一些行少了一些行都正常，SYMBOL与ENTREZID不是一对一的。
nrow(deg)
save(exp,Group,deg,logFC_t,p_t,file = "step4output.Rdata")

rm(list = ls())  
load(file = 'step4output.Rdata')
library(clusterProfiler)
library(ggthemes)
library(org.Hs.eg.db)
library(dplyr)
library(ggplot2)
library(stringr)
library(enrichplot)

#(1)输入数据
gene_diff = deg$ENTREZID[deg$change != "stable"] 

#(2)富集
#⭐下面的三句都要注意物种
ekk <- enrichKEGG(gene = gene_diff,organism = 'hsa')
#其他物种https://www.genome.jp/kegg/catalog/org_list.html
ekk <- setReadable(ekk,OrgDb = org.Hs.eg.db,keyType = "ENTREZID")
#setReadable这个函数会显示出symbol
#如果ekk是空的，这句就会报错，因为没富集到任何通路。
ego <- enrichGO(gene = gene_diff,OrgDb= org.Hs.eg.db,
                ont = "ALL",readable = TRUE)
#"ALL"意思时BP、CC和MF都选
#setReadable和readable = TRUE都是把富集结果表格里的基因名称转为symbol
class(ekk)

#(3)可视化
#barplot可以换成dotplot
barplot(ego, split = "ONTOLOGY") + 
  facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y") 
barplot(ekk)
#如果ekk中没有padj<0.05的通路，就会报错，因为只画padj<0.05,没有参数