原题目链接：http://www.bio-info-trainee.com/3409.html
原答案：http://bio-info-trainee.com/tmp/basic_visualization_for_expression_matrix.html
参考：https://www.jianshu.com/p/dd4e285665e1
https://www.jianshu.com/p/788010093c90

此文为我在练习过程中的一些见解，半路出家，若有理解错误或者不足的地方，欢迎大家提出一起探讨。

1.安装一些R包
对于not available for R ...或下载很慢等问题可更换镜像。
其他新手常见问题可参考Jimmy老师解答：http://www.bio-info-trainee.com/3925.html

rm(list = ls()) #一键清空
options()$repos 
options()$BioC_mirror
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options()$repos 
options()$BioC_mirror
# https://bioconductor.org/packages/release/bioc/html/GEOquery.html
if (!requireNamespace("BiocManager", quietly = TRUE))
 install.packages("BiocManager")
BiocManager::install("KEGG.db",ask = F,update = F)

各种包的下载，并用library载入：

BiocInstaller::biocLite('CLL')
install.packages('corrplot')
install.packages('gpairs')
install.packages('vioplot')
BiocManager::install(c("GSEABase","GSVA","clusterProfiler" ),ask = F,update = F)
BiocManager::install(c("GEOquery","limma","impute" ),ask = F,update = F)
BiocManager::install(c("org.Hs.eg.db","hgu133plus2.db" ),ask = F,update = F)
options()$repos
install.packages('WGCNA')
install.packages(c("FactoMineR", "factoextra"))
install.packages(c("ggplot2", "pheatmap","ggpubr"))
library("FactoMineR")
library("factoextra")
library(GSEABase)
library(GSVA)
library(clusterProfiler)
library(ggplot2)
library(ggpubr)
library(hgu133plus2.db)
library(limma)
library(org.Hs.eg.db)
library(pheatmap)
BiocManager::install('cLL')
library(CLL)
library(ggplot2)
library(reshape2)
library(gpairs)
library(corrplot)

2.提取CLL数据包中自带数据文件

data(sCLLex)#加载内置测试数据
exprs(sCLLex)#提取表达矩阵
dim(sCLLex)#查看行、列数，共12625个探针（行）在22个样本（列）中的表达量
sCLLex=sCLLex[,1:8]#只分析前八个样本，也可对所有样本做分析

3.了解 str,head,help函数

str(sCLLex)#显示对象内部结构，即对象中有什么
head(sCLLex)#与str类似，区别在于使用head函数时，若行、列信息过多则只会显示部分
help("sCLLex")#不懂的时候问R这个函数或包是什么，怎么用

4.安装并了解 hgu95av2.db 包,看看 ls("package:hgu95av2.db") 后显示的那些变量

BiocManager::install('hgu95av2.db') 
library(hgu95av2.db)#注意安装一个包后要加载出来
ls("package:hgu95av2.db")#hgu95av2.db是一个注释包，它为hgu95av2平台的芯片提供注释，这个包中有很多注释文件,注意自己的测序数据用的哪个注释文件

5.找到 TP53 基因对应的探针ID

head(toTable(hgu95av2SYMBOL)) #hgu95av2SYMBOL是一种注释，即探针与基因的对应关系，head函数查看前几行
ids=toTable(hgu95av2SYMBOL)#将注释信息建一个数据框，方便使用
dim(ids)
library(dplyr)#使用该包中的filter函数过滤数据
filter(ids,symbol=='TP53')#找到ids中symbol为TP53所对应的探针

6.理解探针与基因的对应关系
根据dim()函数可查看到测序信息文件中一共测得了12625个探针的表达量，但我们使用的注释信息中只有11459个探针所对应的基因名，因此在处理信息之前需要将这些没有注释信息的探针表达量去除。
我们在上一题也看到了一个TP53对应了多个探针，这也是正常的，因为有些基因过长，通常会设计好几个探针，最后取表达量最高值、平均值、中位数等作为该基因的表达量均可。

table(table(ids$symbol)>1)#统计某个基因是否重复，发现2030个基因有一个以上的探针。第一个table判断ids中的symbol的个数是否大于1，返回T or F，第二个table统计F和T各有多少个。
table(table(ids$symbol)>1)#统计某个基因是否重复，发现2030个基因有一个以上的探针
table(ids$symbol)#统计symbol中每个元素出现的次数
max(table(ids$symbol))#symbol中元素出现的最大次数
table(table(ids$symbol)==8)#统计symbol中元素出现次数为8的有几个
sort(table(ids$symbol),decreasing = T)[1:5]#将symbol按出现次数从大到小排序，最大的5个表示出来，但我觉得还有其他方法直接显示出现次数为8的是哪几个

7.8.9.10.第二步提取到的表达矩阵是12625个探针在22个样本的表达量矩阵，找到那些不在 hgu95av2.db 包收录的对应着SYMBOL的探针。删除没有对应基因名字的探针。整合表达矩阵，多个探针对应一个基因的情况下，只保留在所有样本里面平均表达量最大的那个探针。将表达矩阵的行名更换为对应的基因名

exprSet=exprs(sCLLex)#用exprs函数提取sCLLex中的表达信息，建一个数据框
dim(exprSet)
dim(ids)

测序结果有12625个探针的表达量，探针注释信息只有11459个探针，接下来去除无对应基因名称的探针及其表达量。

rownames(exprSet)%in% ids$probe_id#a %in% table a是否位于table中，返回T or F。这里判断exprSet的行名是否位于ids的probe_id这一列中   
table(rownames(exprSet)%in% ids$probe_id)#统计T、F的个数，有1166个探针不存在对应的基因名
n_exprSet=exprSet[!(rownames(exprSet)%in% ids$probe_id)==F,]#将包含在探针包里的探针测序数据新建一个数据框，这里很奇怪，因该是满足条件的包括进来，这里却填F的时候会包括进来
dim(n_exprSet)#dim.data.frame与dim有区别
View(n_exprSet)
tmp = by(n_exprSet,ids$symbol,
         function(x) 
  rownames(x)[which.max(rowMeans(x))] )#by函数，将对象按行或某种方式分为一个个小的子集，对每个子集进行操作。这里将对n_exprSet的行按与ids中的symbol的对应关系，将同一个symbol的探针放在一起，然后对子集中的每一行取平均值，返回每个子集中平均值大的那一行的行名（探针名）
probes = as.character(tmp)#定义为字符型
View(probes)#这里就是样本中相同基因平均值最大的探针
exprSet=exprSet[rownames(exprSet) %in% probes ,]#对exprSet的行进行操作，若改行的行名在probes中出现了则保留下来
dim(exprSet)
View(exprSet)
rownames(exprSet)=ids[match(rownames(exprSet),ids$probe_id),2]#将exprSet的行名与ids的probe_id相互匹配，若相同，则exprSet的行名为ids的第二列
View(exprSet)

11.对第10步得到的表达矩阵进行探索，画表达量的boxplot、hist、density图
低阶版：

library(ggplot2)
boxplot(exprSet)#各个样本表达量的箱型图
hist(exprSet)#表达量的频率分布直方图

image.png

image.png

进阶版：

#画图进阶版
library(reshape2)#使用该包中的melt函数
exprSet_L=melt(exprSet)#处理数据，将多维数据转变为一维数据，即宽数据变为长数据
colnames(exprSet_L)=c('probe','sample','value')#每一列依次为探针（基因）、样本来源、表达量
group_list=sCLLex$Disease#提取分组信息
group_list#查看每个样本所属的组别
exprSet_L$group=rep(group_list,each=nrow(exprSet))#将样本按stable、progres.分组
head(exprSet_L)
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()#画图，横轴为样本、纵轴为表达量、填充色按组别来区分
print(p)

image.png

p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_violin()#同理，画小提琴图
print(p)

image.png

p=ggplot(exprSet_L,aes(value,fill=group))+geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4)#表达量分布直方图，其中bins参数为对x轴进行划分的个数，geom_histogram等里面的参数若不规定，则会自动取最优，facet_wrap(~sample, nrow = 4)对结果图进行排版，将图分为几行排放
print(p)

image.png

p=ggplot(exprSet_L,aes(value,col=group))+geom_density()+facet_wrap(~sample, nrow = 4)#密度图
print(p)

image.png

p=ggplot(exprSet_L,aes(value,col=group))+geom_density() #将样本的表达量密度图按组区分画在一张图里
print(p)

image.png

12.计算每个基因在所有样本中的mean,median,max,min,sd,var,mad，最后按照mad值排序，取top 50 mad值的基因名，得到列表

#计算每个基因的统计学指标
exmean=apply(exprSet,1,mean)#取exprSet的每一行（1表示行，2表示列）的平均值，下同，不再一一赘述
View(exmean)
exmax=apply(exprSet,1,max)
View(exmax)
exmad=apply(exprSet,1,mad)
View(exmad)
madname=names(sort(apply(exprSet,1,mad),decreasing=T)[1:50])#选取mad top50的基因的名字

13.得到top 50 mad值的基因列表来取表达矩阵的子集，并且热图可视化子表达矩阵

exprSetmad=exprSet[madname,]#exprSetmad为包含madname的行
View(exprSetmad)
#各种形式的热图
heatmap(exprSetmad)
pheatmap::pheatmap(exprSetmad)

14.取不同统计学指标mean,median,max,mean,sd,var,mad的各top50基因列表，使用UpSetR包来看他们之间的overlap情况
需要用到UpSetR包。

#取各种指标的前50的基因的名字
meanname=names(sort(apply(exprSet,1,mean),decreasing=T)[1:50])
medianname=names(sort(apply(exprSet,1,median),decreasing=T)[1:50])
minname=names(sort(apply(exprSet,1,min),decreasing=T)[1:50])
maxname=names(sort(apply(exprSet,1,max),decreasing=T)[1:50])
sdname=names(sort(apply(exprSet,1,sd),decreasing=T)[1:50])
varname=names(sort(apply(exprSet,1,var),decreasing=T)[1:50])

BiocManager::install("UpSetR")
library(UpSetR)
all_name=c(madname,meanname,medianname,minname,maxname,sdname,varname)#建一个集合，该集合中包括7个子集中的共有元素
edat=data.frame(all_name,
                e_mean=ifelse(all_name %in%meanname,1,0),
                e_median=ifelse(all_name %in% medianname ,1,0),
                e_max=ifelse(all_name %in% maxname ,1,0),
                e_min=ifelse(all_name%in% minname,1,0),
                e_sd=ifelse(all_name %in% sdname,1,0),
                e_var=ifelse(all_name%in% varname ,1,0),
                e_mad=ifelse(all_name %in% madname ,1,0)
)#逐步判断共有集合中的元素是否位于子集中并返回yes（1）or no（0），这一步不是特别理解
upset(edat,nsets=7)#用upset进行画图，目标数据框是edat，包含7个子集，其他参数可以突出显示颜色等

image.png

15.在第二步的基础上面提取CLL包里面的data(sCLLex) 数据对象的样本的表型数据
提取样本表型数据我们在第11题的时候也提取过一次，这里主要是想将样本名更换为组别的信息

pdata=pData(sCLLex)#提取样本表型数据
group_list=as.character(pdata[,2])#提取样本表型，准备做该样本的名称
group_list

以下几种分析方法大家最好了解一下原理，便于理解
16.聚类并且绘图，添加样本的临床表型数据信息(更改样本名)

colnames(exprSet)=paste(group_list,1:8,sep='')#样本赋予名称，便于识别类
out.dist=dist(t(exprSet),method='euclidean')#数值变距离，在聚类分析中有几类求两点距离的方法，这里采用平方欧氏距离来处理
out.hclust=hclust(out.dist,method='complete')#几种聚类方法，这里以complete法进行聚类
plot(out.hclust)

image.png

#下面这种好看一些
nodePar <- list(lab.cex = 0.6, pch = c(NA, 19),
                cex = 0.7, col = "blue")
hc=hclust(dist(t(exprSet)))
par(mar=c(5,5,5,10))
plot(as.dendrogram(hc), nodePar = nodePar, horiz = TRUE)

image.png

17.PCA绘图

#PCA主成分分析,通过正交变换将一组可能存在相关性的变量转换为一组线性不相关的变量，转换后的这组变量叫主成分。主成分分析是对于原先提出的所有变量，将重复的变量（关系紧密的变量）删去多余，建立尽可能少的新变量，使得这些新变量是两两不相关的，而且这些新变量在反映课题的信息方面尽可能保持原有的信息。
library(ggplot2)
library("FactoMineR")#负责分析
library("factoextra")#负责可视化
res.pca <- PCA(exprSet, graph = FALSE)#标准化
print(res.pca)#查看分析结果
eig.val <- get_eigenvalue(res.pca)
eig.val#各个组成分对方差的解释程度，主成分1解释了98%的方差，通过观察特征值来考虑保留几个组成分。大于1的eigenvalue说明了这个主成分不止解释了一个原始变量
fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 50))#将方差的解释程度用图表示出来
var <- get_pca_var(res.pca)
var#提取变量的分析结果
head(var$coord)
head(var$cos2)
head(var$contrib)
fviz_pca_var(res.pca, col.var = "black")#相关曲线作图。变量组内包括和主成分之间的关系，正相关的变量是彼此靠近的，负相关的变量师南辕北辙的，而从中心点到变量的长度则代表着变量在这个维度所占的比例（也可以理解为质量，quality）。
library("corrplot")
corrplot(var$cos2, is.corr=FALSE)#各个组成分代表质量作图
fviz_cos2(res.pca, choice = "var", axes = 1:2)#用factoextra里面的fviz_cos()函数作代表质量图，没上图好看
fviz_pca_var(res.pca, col.var = "cos2",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), 
             repel = TRUE # Avoid text overlapping
)#上色
fviz_pca_var(res.pca, alpha.var = "cos2")#通过改变变量的透明度来说明其重要性
corrplot(var$contrib, is.corr=FALSE) #各个变量对主成分的贡献程度
fviz_contrib(res.pca, choice = "var", axes = 1, top = 10)
fviz_contrib(res.pca, choice = "var", axes = 2, top = 10)#同上，用factoextra里的函数作图，每个变量对某个主成分的贡献程度

上面主要是查看各个组成分对整体的贡献程度，没有放图。

#也可以用下面方法做PCA图，更好看
pc <- prcomp(t(exprSet),scale=TRUE) 
pcx=data.frame(pc$x)
pcr=cbind(samples=rownames(pcx),group_list, pcx) 
p=ggplot(pcr, aes(PC1, PC2))+geom_point(size=5, aes(color=group_list)) +#能最大反映样本差异性的两个成分（PC1、PC2）
  geom_text(aes(label=samples),hjust=-0.1, vjust=-0.3)#label=samples可以加上样本名称
print(p)

image.png

18.T检验

#T test
library(readxl)
library(car)#调用car判断方差齐性
exprSet<-as.data.frame(exprSet)
gl=as.factor(group_list) #将最初得到的的表型数据因子化
group1 = which(group_list == levels(gl)[1]) #levels(group_list)[1]返回第一个因子progres,从group_list中选出progres的元素，用which来获取他们所在的位置【目的是为下面分别得到两种表型的样本作准备】
group2 = which(group_list == levels(gl)[2]) #返回第二个因子stable
et1 = exprSet[, group1] #将表型为progres的样本选出来，因为这里是要求t值，可以命名为e矩阵的t值，即et
et2 = exprSet[, group2] #将表型为stable的样本选出来
et = cbind(et1, et2) #按列合并
pvals = apply(exprSet, 1, function(x){
  t.test(as.numeric(x)~group_list)$p.value # 多组样本的t检验
})
p.adj = p.adjust(pvals, method = "BH") #多重比较时校正p值
eavg_1 = rowMeans(et1) #progres是对照组
eavg_2 = rowMeans(et2) #stable是使用药物处理后的——处理组
log2FC = eavg_2-eavg_1
DEG_t.test = cbind(eavg_1, eavg_2, log2FC, pvals, p.adj)
DEG_t.test=DEG_t.test[order(DEG_t.test[,4]),] #从小到大排序
DEG_t.test=as.data.frame(DEG_t.test)#检验结果
DEG_t.test
head(DEG_t.test)

19.用limma包对表达矩阵及样本分组信息进行差异分析，得到差异分析表格，画火山图

#火山图
suppressMessages(library(limma))
design=model.matrix(~factor(group_list))#分组信息
fit=lmFit(exprSet,design)
fit=eBayes(fit)
DEG=topTable(fit,coef=2,n=Inf)
with(DEG, plot(logFC, -log10(P.Value), pch=20, main="Volcano plot")) 

image.png

#带标题、上下调基因的火山图
logFC_cutoff <- with(DEG,mean(abs( logFC)) + 2*sd(abs( logFC)) )
DEG$change = as.factor(ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,
                              ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'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',])
)#图片标题
library(ggplot2)
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)
print(g)

image.png

20.对T检验结果的P值和limma包差异分析的P值画散点图，看看哪些基因相差很大

DEG_t.test#T检验的差异分析结果
DEG#limma包分析的差异结果

拿LSM2这个基因来说，两种分析结果里面的p adjust差距还是很大的。
可参考：https://www.jianshu.com/p/84bfdf91118c
（写到这儿有些累了，以后理解了再更新）