suppressPackageStartupMessages(library(CLL))
data(sCLLex)#data函数是加载数据集,类似于library（）
exprSet=exprs(sCLLex)#提取表达矩阵
samples=sampleNames(sCLLex)
pdata=pData(sCLLex)#提取样本的表型数据,pData用来获得每个样本的描述信息，用来获得group_list
#即取出每个病人的id以及每个人的疾病进展情况，
#第一列为病人ID，第二列为病人进展情况
group_list=as.character(pdata[,2])#取出每个人的疾病进展情况，并将其转化为字符
##group_list是分组信息，有了分组信息和表达矩阵就可以做下面的差异分析了。
##用分组信息给表达矩阵的列改名
dim(exprSet)
exprSet[1:5,1:5]
#自动提取注释信息注释探针
show(sCLLex)#show函数查看芯片种类，查看芯片信息。
#使用这个函数查看芯片信息的时候，如果没有下载对应的注释包的话会
#自动下载，免掉了自己上bioconductor查找下载注释包的麻烦
#annotation：hgu95av2这一行就是显示芯片种类为hgu95av2
#对应的注释包为hgu95av2.db
###运行以下代码得到探针注释信息probe2symbol和表达矩阵exprSet
library(stats4)
library(S4Vectors)
library(IRanges)
library(AnnotationDbi)
library(org.Hs.eg.db)
library(hgu95av2.db)#加载芯片注释包
columns(hgu95av2.db)#查看芯片注释包提供的注释信息
ids=toTable(hgu95av2SYMBOL)#探针id和基因symbol，
#toTable提取注释信息，这里选择GENESYMBOL进行注释，也可以选择其他ID进行注释，详情请见注释包的支持文档
#toTable提取对象的列，hgu95av2SYMBOL列就是一个两列多行的数据框，是探针对应symbol文件
length(unique(ids$symbol))#查看有多少个symbol：[1] 8585
length(unique(ids$probe_id))#查看有多少个probe_id：[1] 11460
####注意，一般情况下，多个探针对应一个基因，上面的统计也验证了这事实
tail(sort(table(ids$symbol)))#查看每个symbol的频次，输出出现频率最高的几个基因
#YME1L1 GAPDH INPP4A MYB PTGER3 STAT1
#7 8 8 8 8 8
#可以发现 STAT1该symbol对应了8个探针
table(sort(table(ids$symbol)))#统计各频率出现的次数
# 1 2 3 4 5 6 7 8
#6555 1428 451 102 22 16 6 5
#可以发现，有6555个symbol他们只对应一个探针
# 相反，有5个symbol他们对应了8个探针
plot(table(sort(table(ids$symbol))))
table(rownames(exprSet) %in% ids$probe_id)
#查有表达谱内多少个基因再探针库里，匹配不上要把这些基因删除
#FALSE TRUE
#1165 11460
#有1165个基因不在探针库里
exprSet=exprSet[rownames(exprSet) %in% ids$probe_id,]
dim(exprSet)
#[1] 11460 22
ids=ids[match(rownames(exprSet),ids$probe_id),]
#match返回a在b中的位置
#返回表达谱内每个探针在探针库里的位置，并将其从探针库里连同对应的symbol取出来
head(ids)
exprSet[1:5,1:5]
tmp = by(exprSet,ids$symbol,function(x) rownames(x)[which.max(rowMeans(x))] )
#取均值最大的探针作为该基因对应的探针
probes = as.character(tmp)
exprSet=exprSet[rownames(exprSet) %in% probes ,]
dim(exprSet)
rownames(exprSet)=ids[match(rownames(exprSet),ids$probe_id),2]
exprSet[1:5,1:5]
library(reshape2)
exprSet_L=melt(exprSet)
#拆分数据，把整体数据打碎，让其回到一个个数据点的状态，然后再根据观测id和某变量名称，
#即可交叉定位到某数据点
colnames(exprSet_L)=c('probe','sample','value')
#打碎后的矩阵没有列明，这一步是赋予其列名
exprSet_L$group=rep(group_list,each=nrow(exprSet))
#给exprSet_L加group列，对每一个样本，都对应一种疾病进展状态，
#而每一个样本对应有nrow(exprSet)个基因，所以一个病人的疾病进展状态重复nrow(exprSet)次
head(exprSet_L)
### ggplot2
library(ggplot2)
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
print(p)
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_violin()
print(p)
p=ggplot(exprSet_L,aes(value,fill=group))+geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4)
#sample组在图中并排出现，不是重叠为单一的图形，这样关系更清晰。
#sample组按行展开，一共四行
print(p)
p=ggplot(exprSet_L,aes(value,col=group))+geom_density()+facet_wrap(~sample, nrow = 4)
print(p)
p=ggplot(exprSet_L,aes(value,col=group))+geom_density()
#col=color,边框颜色按group属性来分
print(p)
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
print(p)
p=p+stat_summary(fun.y="mean",geom="point",shape=23,size=3,fill="red")
print(p)
p=p+theme_set(theme_set(theme_bw(base_size=20)))
print(p)
p=p+theme(text=element_text(face='bold'),axis.text.x=element_text(angle=30,hjust=1),axis.title=element_blank())
print(p)
## mean,median,max,min,sd,var,mad
g_mean <- tail(sort(apply(exprSet,1,mean)),50)
g_median <- tail(sort(apply(exprSet,1,median)),50)
g_max <- tail(sort(apply(exprSet,1,max)),50)
g_min <- tail(sort(apply(exprSet,1,min)),50)
g_sd <- tail(sort(apply(exprSet,1,sd)),50)#标准差
g_var <- tail(sort(apply(exprSet,1,var)),50)#方差
g_mad <- tail(sort(apply(exprSet,1,mad)),50)#绝对中位差
g_mad
names(g_mad)
## heatmap热图的绘制，
#首先要对表达谱矩阵进行转置，按列（行为样本，列为基因）进行标准化,
#再进行转置
#最后对其使用pheatmap函数
library(pheatmap)
choose_gene=names(tail(sort(apply(exprSet,1,mad)),50))#绝对中位差
choose_matrix=exprSet[choose_gene,]
choose_matrix=t(scale(t(choose_matrix)))
#scale()为数据对象x按列进行中心化(cener=TRUE)
#或者标准化（center=TRUE，scale=TRUE）
pheatmap(choose_matrix)
## UpSetR
# https://cran.r-project.org/web/packages/UpSetR/README.html
#取不同统计学指标mean,median,max,mean,sd,var,mad的各top50基因列表，
#使用UpSetR包来看他们之间的overlap情况。
library(UpSetR)
g_all <- unique(c(names(g_mean),names(g_median),names(g_max),names(g_min),
names(g_sd),names(g_var),names(g_mad)))
dat=data.frame(g_all=g_all,
g_mean=ifelse(g_all %in% names(g_mean) ,1,0),
g_median=ifelse(g_all %in% names(g_median) ,1,0),
g_max=ifelse(g_all %in% names(g_max) ,1,0),
g_min=ifelse(g_all %in% names(g_min) ,1,0),
g_sd=ifelse(g_all %in% names(g_sd) ,1,0),
g_var=ifelse(g_all %in% names(g_var) ,1,0),
g_mad=ifelse(g_all %in% names(g_mad) ,1,0)
)
upset(dat,nsets = 7)
pdata=pData(sCLLex)#提取样本的表型数据
group_list=as.character(pdata[,2])
group_list
dim(exprSet)
exprSet[1:5,1:5]
## hclust实现层次聚类
hc=hclust(dist(t(exprSet)))
plot(hc)
##是对病人样本进行聚类，所以先要进行转置
#stats包中的hclust（）函数进行 行聚类。
#系统聚类一般首先使用dist（）函数计算欧式距离，得到的是一个下三角矩阵
#再使用hclust（）函数展开系统聚类。
#上面得到的病人聚类图是对病人id进行聚类，我们可以将每个病人的病情进行层次聚类
colnames(exprSet)=paste(group_list,1:22,sep='')
#更改原始表达矩阵exprSet的列名，改过之后我们就可以知道第1-22位病人每一个病人的疾病进展情况
#将group_list里的每个元素和1：22里的每个元素结合起来，
#以‘’为分隔符，即不添加分隔符
#这一步是添加样本的临床表型数据信息(更改样本名)
# Define nodePar
nodePar <- list(lab.cex = 0.6, pch = c(NA, 19),
cex = 0.7, col = "blue")
#list(A=a,B=b...)不再是常规的建立列表[[1]]、[[2]]...
#而是给每个列表都命名，如A，B...，以后可以通过这些名称来访问
hc=hclust(dist(t(exprSet)))
par(mar=c(5,5,5,10))
#par（）调整图形设置，我们可以通过设定par函数的各个参数来调整我们的图形
#mar 设置图形空白边界行数，
#mar = c(bottom, left, top, right)。缺省为mar = c(5.1,4.1,4.1,2.1)。
plot(as.dendrogram(hc), nodePar = nodePar, horiz = TRUE)
#plot(x1,y1,main = "我是标题",xlab = "我是x轴",ylab = "我是y轴",xlim = c(0,100),ylim = c(0,100),col = "red",pch = 19)
#as.dendrogram()将hclust生成的对象转换为另类的聚类图。
## PCA
library(ggfortify)
df=as.data.frame(t(exprSet))#df行为疾病情况，列为基因symbol
df$group=group_list#给df加group列，group列为每一个病人的疾病进展情况
View(df)
autoplot(prcomp(df[,1:(ncol(df)-1)] ), data=df,colour = 'group')
#prcomp:在给定数据矩阵上进行主成分分析，返回一个类prcomp对象
## t.test
dat = exprSet
group_list=as.factor(group_list)
group1 = which(group_list == levels(group_list)[1])
#levels(group_list)[1]取因子元素1
#which（A==a）找出A中为a的元素，将他们的位置序号作为输出
#找出疾病进展情况列为Progress的
group2 = which(group_list == levels(group_list)[2])
#找出疾病进展情况列为stable的
dat1 = dat[, group1]#取出progress情况的病人的基因表达
dat2 = dat[, group2]#取出stable情况的病人的基因表达
dat = cbind(dat1, dat2)#横向合并，行名为基因symbol
#第一列为progress情况的病人的基因表达谱
#第二列为stable情况的病人的基因表达谱
#######差异分析：分为两种
#######一种是具体步骤
#######另一种是用limma包做差异分析
####第一种差异表达分析：两类样本之间的差异表达
#对于每个基因，比较其在progres.和stable两种样本有无差异
pvals = apply(exprSet, 1, function(x){
t.test(as.numeric(x)~group_list)$p.value
})
#取P值，一行一行求P值
#对基因表达谱里的某一个基因在各样本里的表达差异，
#group_list是各样本的疾病情况种类，一共有两个
#将每一行按照group——list里的factor分成两类：group progres.和group stable，将这两个类进行T检验
#查看每个基因在两个类里表达差异对应的P值，即类的P值
#as.numeric()将字符型的数字转换成数值型的数字
#P值校正：在我们对数据（尤其是数据量较大的时候）进行
#假设检验以后（也就是我们常称的“多重检验”），得到了对应的P值，
#这时候不要高兴太早，你还没有分析完，接下来，你其实还要对
#得到的P值进行校正，然后使用校正后的P值进行后续的分析。
#那么，问题就来了：我们为什么要对多重检验得到的P值进行校正呢？
#其实关键点还是在这个“多”字上！中国老话都说了，“常在河边走哪有不湿鞋”，“言多必失”
p.adj = p.adjust(pvals, method = "BH")
#p.adjust对得到的每一个P值进行校正，采用BH方法，BH是最常用的校正方法
avg_1 = rowMeans(dat1)#对progress所有病人的每个基因的表达情况求均值
avg_2 = rowMeans(dat2)#对stable所有病人的每个基因的表达情况求均值
log2FC = avg_2-avg_1#foldchange计算，差异表达分析
DEG_t.test = cbind(avg_1, avg_2, log2FC, pvals, p.adj)
DEG_t.test=DEG_t.test[order(DEG_t.test[,4]),]
#两种样本基因表达情况很有差异的基因（p值小的）的排在最前面
#p值越小，差异越显著
DEG_t.test=as.data.frame(DEG_t.test)
head(DEG_t.test)#t检验结果
##第二种差异表达分析：使用limma包对表达矩阵及样本分组信息进行差异分析，
#得到差异分析表格，重点看logFC和P值，画个火山图(就是logFC和-log10(P值)的散点图。)。
# DEG by limma
########差异分析——————需要对比矩阵、设计矩阵、表达矩阵
###一、构建实验设计矩阵design
suppressMessages(library(limma))
design <- model.matrix(~0+factor(group_list))
#model.matrix将里面的因子转化为虚拟变量，这里有两个虚拟变量
#在progress列将group_list的progress.因子转换成1，不是progress则为0，
#在stable列将group_list的stable因子转换成1，不是stable则为0
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design
# progres. stable
#CLL11.CEL 1 0
#CLL12.CEL 0 1
#CLL13.CEL 1 0
#CLL14.CEL 1 0
#CLL15.CEL 1 0
#CLL16.CEL 1 0
#CLL17.CEL 0 1
#CLL18.CEL 0 1
#CLL19.CEL 1 0
#CLL20.CEL 0 1
#CLL21.CEL 1 0
#CLL22.CEL 0 1
#CLL23.CEL 1 0
#CLL24.CEL 0 1
#CLL2.CEL 0 1
#CLL3.CEL 1 0
#CLL4.CEL 1 0
#CLL5.CEL 1 0
#CLL6.CEL 1 0
#CLL7.CEL 1 0
#CLL8.CEL 1 0
#CLL9.CEL 0 1
###二、构建对比矩阵，比较两个实验条件下表达数据
contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
#构建对比矩阵，这个对比矩阵会声明我们要把progress.组跟stable进行差异分析比较
contrast.matrix
# Contrasts
#Levels progres.-stable#说明stable是用来比的，control是用来被比的
#progres. 1（1-0）
#stable -1（0-1）
##这个矩阵声明，我们要把progress.组跟stable进行差异分析比较
###三、差异分析（对比矩阵、设计矩阵、表达矩阵已经全部准备好）
##step1 线性模型拟合：设计矩阵+表达矩阵
fit <- lmFit(exprSet,design)
##step2 根据对比模型进行差值计算：线性拟合矩阵+对比矩阵
fit2 <- contrasts.fit(fit, contrast.matrix)##这一步很重要，大家可以自行看看效果
#fit的coficient矩阵的progress列和stable列相减，
#就是对比矩阵里的progress列对应1，stable列对应-1，即完成了coficient矩阵的progress列和stable列相减
##step3 贝叶斯检验
fit2 <- eBayes(fit2) ## default no trend !!!
##eBayes() with trend=TRUE
##step4 生成所有基因的检验结果报告
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG=na.omit(tempOutput)#删除缺失值，删除有空值的行
#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(nrDEG)
## volcano plot
DEG=nrDEG
logFC_cutoff <- with(DEG,mean(abs( logFC)) + 2*sd(abs( logFC)) )
#with把所有操作都限制在数据框上，调用数据框DEG里的logFC
DEG$change = as.factor(ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,
ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
)
#对DEG的每一行判定为UP/DOWN/NOT，判定结果添加为change列，并转化为因子格式
#P.Value < 0.05且|logFC|>logFC_cutoff且logFC > logFC_cutoff，输出'UP'
#P.Value < 0.05且|logFC|>logFC_cutoff且logFC </=logFC_cutoff，输出'DOWN'
#P.Value >/= 0.05或|logFC|</=logFC_cutoff，输出'NOT'
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
'\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
'\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',])
)
#round四舍五入，round（x,n)对x四舍五入，并保留3个小数点
###差异分析，得到差异分析表格，重点看logFC和P值，画个火山图(就是logFC和-log10(P值)的散点图）
g = ggplot(data=DEG, aes(x=logFC, y=-log10(P.Value), color=change)) +
geom_point(alpha=0.4, size=1.75) +
theme_set(theme_set(theme_bw(base_size=20)))+
xlab("log2 fold change") + ylab("-log10 p-value") +
ggtitle( this_tile ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
scale_colour_manual(values = c('blue','black','red')) ## corresponding to the levels(res$change)
#scale_colour_manual手动设置颜色
#alpha设置透明度
###每个点代表一个基因
#对P值取-log10，是由于P值越小表示越显著，所以我们进行-log10（P value）转化后，转化值越大表示差异约显著
#对foldchange取log2，缩小差异表达倍数的差距
###所以呈现出来的火山图是纵坐标越大，P值越显著，越好；横坐标大于1，差异表达越显著
###热图解读：在X=-1左侧的点是下调2倍以上的基因，在X=1右侧的点是上调2倍以上的基因。同时，平行于X轴有一条虚线Y=1.30，即-log10(0.05），在虚线以上的点表示显著性<0.05的基因。
#所以热土重点看x=1右侧和y=1.3上上侧，都是差异表达显著的基因
print(g)
####对T检验结果的P值和limma包差异分析的P值画散点图，
#看看哪些基因相差很大。
### different P values
head(nrDEG)
head(DEG_t.test)
DEG_t.test=DEG_t.test[rownames(nrDEG),]
head(DEG_t.test)
plot(DEG_t.test[,3],nrDEG[,1])
#看DEG_t.test$log2FC和 nrDEG$logFC差异
plot(DEG_t.test[,4],nrDEG[,4])
#看DEG_t.test$pvals和 nrDEG$p.Value差异，结果是几乎在y=x上，所以几乎无差异
plot(-log10(DEG_t.test[,4]),-log10(nrDEG[,4]))
#看DEG_t.test$pvals和 nrDEG$p.Value的-log10差异，结果是几乎在y=x上，所以几乎无差异
exprSet['GAPDH',]
exprSet['ACTB',]
exprSet['DLEU1',]
library(ggplot2)
library(magrittr)
library(ggpubr)
my_comparisons <- list(
c("stable", "progres.")
)
dat=data.frame(group=group_list,
sampleID= names(exprSet['DLEU1',]),
values= as.numeric(exprSet['DLEU1',]))
#data.frame()函数不是转换数据框函数，而是组合数据框的函数，
#指明了每一列的组成，dat数据框有三列，group、sampleID、values三列
#画DLEU1基因的表达箱型图
ggboxplot(
dat, x = "group", y = "values",
color = "group",
add = "jitter"
)+
stat_compare_means(comparisons = my_comparisons, method = "t.test")
## heatmap
library(pheatmap)
choose_gene=head(rownames(nrDEG),25)
#nrDEG是经过差异表达分析后按P值大小顺序排序的矩阵
#p值越小，差异越显著
choose_matrix=exprSet[choose_gene,]#选择差异表达后排名靠前的基因
choose_matrix=t(scale(t(choose_matrix)))
#scale函数，标准化函数，对数据列进行正态化处理，默认情况下是将一组数的每个数都减去这组数的平均值后再除以这组数的均方根。
#对每一个基因列标准化
pheatmap(choose_matrix)
