原題目鏈接：http://www.bio-info-trainee.com/3409.html
原答案：http://bio-info-trainee.com/tmp/basic_visualization_for_expression_matrix.html
參考：http://www.reibang.com/p/dd4e285665e1
http://www.reibang.com/p/788010093c90

此文為我在練習(xí)過程中的一些見解，半路出家尊流，若有理解錯(cuò)誤或者不足的地方凭戴，歡迎大家提出一起探討。

1.安裝一些R包
對(duì)于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數(shù)據(jù)包中自帶數(shù)據(jù)文件

data(sCLLex)#加載內(nèi)置測(cè)試數(shù)據(jù)
exprs(sCLLex)#提取表達(dá)矩陣
dim(sCLLex)#查看行、列數(shù)基协，共12625個(gè)探針（行）在22個(gè)樣本（列）中的表達(dá)量
sCLLex=sCLLex[,1:8]#只分析前八個(gè)樣本，也可對(duì)所有樣本做分析

3.了解 str,head,help函數(shù)

str(sCLLex)#顯示對(duì)象內(nèi)部結(jié)構(gòu)菇用，即對(duì)象中有什么
head(sCLLex)#與str類似澜驮，區(qū)別在于使用head函數(shù)時(shí)，若行惋鸥、列信息過多則只會(huì)顯示部分
help("sCLLex")#不懂的時(shí)候問R這個(gè)函數(shù)或包是什么杂穷，怎么用

4.安裝并了解 hgu95av2.db 包,看看 ls("package:hgu95av2.db") 后顯示的那些變量

BiocManager::install('hgu95av2.db') 
library(hgu95av2.db)#注意安裝一個(gè)包后要加載出來
ls("package:hgu95av2.db")#hgu95av2.db是一個(gè)注釋包，它為hgu95av2平臺(tái)的芯片提供注釋卦绣，這個(gè)包中有很多注釋文件,注意自己的測(cè)序數(shù)據(jù)用的哪個(gè)注釋文件

5.找到 TP53 基因?qū)?yīng)的探針I(yè)D

head(toTable(hgu95av2SYMBOL)) #hgu95av2SYMBOL是一種注釋耐量，即探針與基因的對(duì)應(yīng)關(guān)系，head函數(shù)查看前幾行
ids=toTable(hgu95av2SYMBOL)#將注釋信息建一個(gè)數(shù)據(jù)框滤港，方便使用
dim(ids)
library(dplyr)#使用該包中的filter函數(shù)過濾數(shù)據(jù)
filter(ids,symbol=='TP53')#找到ids中symbol為TP53所對(duì)應(yīng)的探針

6.理解探針與基因的對(duì)應(yīng)關(guān)系
根據(jù)dim()函數(shù)可查看到測(cè)序信息文件中一共測(cè)得了12625個(gè)探針的表達(dá)量廊蜒，但我們使用的注釋信息中只有11459個(gè)探針?biāo)鶎?duì)應(yīng)的基因名，因此在處理信息之前需要將這些沒有注釋信息的探針表達(dá)量去除溅漾。
我們?cè)谏弦活}也看到了一個(gè)TP53對(duì)應(yīng)了多個(gè)探針山叮，這也是正常的，因?yàn)橛行┗蜻^長(zhǎng)添履，通常會(huì)設(shè)計(jì)好幾個(gè)探針屁倔，最后取表達(dá)量最高值、平均值暮胧、中位數(shù)等作為該基因的表達(dá)量均可锐借。

table(table(ids$symbol)>1)#統(tǒng)計(jì)某個(gè)基因是否重復(fù)，發(fā)現(xiàn)2030個(gè)基因有一個(gè)以上的探針叔壤。第一個(gè)table判斷ids中的symbol的個(gè)數(shù)是否大于1瞎饲，返回T or F，第二個(gè)table統(tǒng)計(jì)F和T各有多少個(gè)炼绘。
table(table(ids$symbol)>1)#統(tǒng)計(jì)某個(gè)基因是否重復(fù)嗅战，發(fā)現(xiàn)2030個(gè)基因有一個(gè)以上的探針
table(ids$symbol)#統(tǒng)計(jì)symbol中每個(gè)元素出現(xiàn)的次數(shù)
max(table(ids$symbol))#symbol中元素出現(xiàn)的最大次數(shù)
table(table(ids$symbol)==8)#統(tǒng)計(jì)symbol中元素出現(xiàn)次數(shù)為8的有幾個(gè)
sort(table(ids$symbol),decreasing = T)[1:5]#將symbol按出現(xiàn)次數(shù)從大到小排序，最大的5個(gè)表示出來，但我覺得還有其他方法直接顯示出現(xiàn)次數(shù)為8的是哪幾個(gè)

7.8.9.10.第二步提取到的表達(dá)矩陣是12625個(gè)探針在22個(gè)樣本的表達(dá)量矩陣驮捍，找到那些不在 hgu95av2.db 包收錄的對(duì)應(yīng)著SYMBOL的探針疟呐。刪除沒有對(duì)應(yīng)基因名字的探針。整合表達(dá)矩陣东且，多個(gè)探針對(duì)應(yīng)一個(gè)基因的情況下启具，只保留在所有樣本里面平均表達(dá)量最大的那個(gè)探針。將表達(dá)矩陣的行名更換為對(duì)應(yīng)的基因名

exprSet=exprs(sCLLex)#用exprs函數(shù)提取sCLLex中的表達(dá)信息珊泳，建一個(gè)數(shù)據(jù)框
dim(exprSet)
dim(ids)

測(cè)序結(jié)果有12625個(gè)探針的表達(dá)量鲁冯，探針注釋信息只有11459個(gè)探針，接下來去除無對(duì)應(yīng)基因名稱的探針及其表達(dá)量色查。

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ǒng)計(jì)T秧了、F的個(gè)數(shù)跨扮，有1166個(gè)探針不存在對(duì)應(yīng)的基因名
n_exprSet=exprSet[!(rownames(exprSet)%in% ids$probe_id)==F,]#將包含在探針包里的探針測(cè)序數(shù)據(jù)新建一個(gè)數(shù)據(jù)框，這里很奇怪验毡，因該是滿足條件的包括進(jìn)來衡创，這里卻填F的時(shí)候會(huì)包括進(jìn)來
dim(n_exprSet)#dim.data.frame與dim有區(qū)別
View(n_exprSet)
tmp = by(n_exprSet,ids$symbol,
         function(x) 
  rownames(x)[which.max(rowMeans(x))] )#by函數(shù)，將對(duì)象按行或某種方式分為一個(gè)個(gè)小的子集晶通，對(duì)每個(gè)子集進(jìn)行操作璃氢。這里將對(duì)n_exprSet的行按與ids中的symbol的對(duì)應(yīng)關(guān)系，將同一個(gè)symbol的探針放在一起录择，然后對(duì)子集中的每一行取平均值拔莱，返回每個(gè)子集中平均值大的那一行的行名（探針名）
probes = as.character(tmp)#定義為字符型
View(probes)#這里就是樣本中相同基因平均值最大的探針
exprSet=exprSet[rownames(exprSet) %in% probes ,]#對(duì)exprSet的行進(jìn)行操作，若改行的行名在probes中出現(xiàn)了則保留下來
dim(exprSet)
View(exprSet)
rownames(exprSet)=ids[match(rownames(exprSet),ids$probe_id),2]#將exprSet的行名與ids的probe_id相互匹配隘竭，若相同，則exprSet的行名為ids的第二列
View(exprSet)

11.對(duì)第10步得到的表達(dá)矩陣進(jìn)行探索讼渊，畫表達(dá)量的boxplot动看、hist、density圖
低階版：

library(ggplot2)
boxplot(exprSet)#各個(gè)樣本表達(dá)量的箱型圖
hist(exprSet)#表達(dá)量的頻率分布直方圖

image.png

image.png

進(jìn)階版：

#畫圖進(jìn)階版
library(reshape2)#使用該包中的melt函數(shù)
exprSet_L=melt(exprSet)#處理數(shù)據(jù)爪幻，將多維數(shù)據(jù)轉(zhuǎn)變?yōu)橐痪S數(shù)據(jù)菱皆，即寬數(shù)據(jù)變?yōu)殚L(zhǎng)數(shù)據(jù)
colnames(exprSet_L)=c('probe','sample','value')#每一列依次為探針（基因）、樣本來源挨稿、表達(dá)量
group_list=sCLLex$Disease#提取分組信息
group_list#查看每個(gè)樣本所屬的組別
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()#畫圖，橫軸為樣本奶甘、縱軸為表達(dá)量篷店、填充色按組別來區(qū)分
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)#表達(dá)量分布直方圖，其中bins參數(shù)為對(duì)x軸進(jìn)行劃分的個(gè)數(shù)疲陕，geom_histogram等里面的參數(shù)若不規(guī)定方淤，則會(huì)自動(dòng)取最優(yōu)，facet_wrap(~sample, nrow = 4)對(duì)結(jié)果圖進(jìn)行排版蹄殃，將圖分為幾行排放
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() #將樣本的表達(dá)量密度圖按組區(qū)分畫在一張圖里
print(p)

image.png

12.計(jì)算每個(gè)基因在所有樣本中的mean,median,max,min,sd,var,mad携茂，最后按照mad值排序，取top 50 mad值的基因名诅岩，得到列表

#計(jì)算每個(gè)基因的統(tǒng)計(jì)學(xué)指標(biāo)
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值的基因列表來取表達(dá)矩陣的子集鸳谜，并且熱圖可視化子表達(dá)矩陣

exprSetmad=exprSet[madname,]#exprSetmad為包含madname的行
View(exprSetmad)
#各種形式的熱圖
heatmap(exprSetmad)
pheatmap::pheatmap(exprSetmad)

14.取不同統(tǒng)計(jì)學(xué)指標(biāo)mean,median,max,mean,sd,var,mad的各top50基因列表，使用UpSetR包來看他們之間的overlap情況
需要用到UpSetR包逮京。

#取各種指標(biāo)的前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)#建一個(gè)集合卿堂，該集合中包括7個(gè)子集中的共有元素
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進(jìn)行畫圖懒棉，目標(biāo)數(shù)據(jù)框是edat草描，包含7個(gè)子集，其他參數(shù)可以突出顯示顏色等

image.png

15.在第二步的基礎(chǔ)上面提取CLL包里面的data(sCLLex) 數(shù)據(jù)對(duì)象的樣本的表型數(shù)據(jù)
提取樣本表型數(shù)據(jù)我們?cè)诘?1題的時(shí)候也提取過一次策严，這里主要是想將樣本名更換為組別的信息

pdata=pData(sCLLex)#提取樣本表型數(shù)據(jù)
group_list=as.character(pdata[,2])#提取樣本表型穗慕，準(zhǔn)備做該樣本的名稱
group_list

以下幾種分析方法大家最好了解一下原理，便于理解
16.聚類并且繪圖妻导，添加樣本的臨床表型數(shù)據(jù)信息(更改樣本名)

colnames(exprSet)=paste(group_list,1:8,sep='')#樣本賦予名稱逛绵，便于識(shí)別類
out.dist=dist(t(exprSet),method='euclidean')#數(shù)值變距離，在聚類分析中有幾類求兩點(diǎn)距離的方法倔韭，這里采用平方歐氏距離來處理
out.hclust=hclust(out.dist,method='complete')#幾種聚類方法术浪，這里以complete法進(jìn)行聚類
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主成分分析,通過正交變換將一組可能存在相關(guān)性的變量轉(zhuǎn)換為一組線性不相關(guān)的變量，轉(zhuǎn)換后的這組變量叫主成分寿酌。主成分分析是對(duì)于原先提出的所有變量胰苏，將重復(fù)的變量（關(guān)系緊密的變量）刪去多余，建立盡可能少的新變量醇疼，使得這些新變量是兩兩不相關(guān)的硕并，而且這些新變量在反映課題的信息方面盡可能保持原有的信息。
library(ggplot2)
library("FactoMineR")#負(fù)責(zé)分析
library("factoextra")#負(fù)責(zé)可視化
res.pca <- PCA(exprSet, graph = FALSE)#標(biāo)準(zhǔn)化
print(res.pca)#查看分析結(jié)果
eig.val <- get_eigenvalue(res.pca)
eig.val#各個(gè)組成分對(duì)方差的解釋程度秧荆，主成分1解釋了98%的方差倔毙，通過觀察特征值來考慮保留幾個(gè)組成分。大于1的eigenvalue說明了這個(gè)主成分不止解釋了一個(gè)原始變量
fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 50))#將方差的解釋程度用圖表示出來
var <- get_pca_var(res.pca)
var#提取變量的分析結(jié)果
head(var$coord)
head(var$cos2)
head(var$contrib)
fviz_pca_var(res.pca, col.var = "black")#相關(guān)曲線作圖乙濒。變量組內(nèi)包括和主成分之間的關(guān)系陕赃，正相關(guān)的變量是彼此靠近的，負(fù)相關(guān)的變量師南轅北轍的，而從中心點(diǎn)到變量的長(zhǎng)度則代表著變量在這個(gè)維度所占的比例（也可以理解為質(zhì)量凯正，quality）毙玻。
library("corrplot")
corrplot(var$cos2, is.corr=FALSE)#各個(gè)組成分代表質(zhì)量作圖
fviz_cos2(res.pca, choice = "var", axes = 1:2)#用factoextra里面的fviz_cos()函數(shù)作代表質(zhì)量圖，沒上圖好看
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) #各個(gè)變量對(duì)主成分的貢獻(xiàn)程度
fviz_contrib(res.pca, choice = "var", axes = 1, top = 10)
fviz_contrib(res.pca, choice = "var", axes = 2, top = 10)#同上廊散，用factoextra里的函數(shù)作圖桑滩，每個(gè)變量對(duì)某個(gè)主成分的貢獻(xiàn)程度

上面主要是查看各個(gè)組成分對(duì)整體的貢獻(xiàn)程度，沒有放圖允睹。

#也可以用下面方法做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)) +#能最大反映樣本差異性的兩個(gè)成分（PC1、PC2）
  geom_text(aes(label=samples),hjust=-0.1, vjust=-0.3)#label=samples可以加上樣本名稱
print(p)

image.png

18.T檢驗(yàn)

#T test
library(readxl)
library(car)#調(diào)用car判斷方差齊性
exprSet<-as.data.frame(exprSet)
gl=as.factor(group_list) #將最初得到的的表型數(shù)據(jù)因子化
group1 = which(group_list == levels(gl)[1]) #levels(group_list)[1]返回第一個(gè)因子progres,從group_list中選出progres的元素缭受，用which來獲取他們所在的位置【目的是為下面分別得到兩種表型的樣本作準(zhǔn)備】
group2 = which(group_list == levels(gl)[2]) #返回第二個(gè)因子stable
et1 = exprSet[, group1] #將表型為progres的樣本選出來胁澳，因?yàn)檫@里是要求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檢驗(yàn)
})
p.adj = p.adjust(pvals, method = "BH") #多重比較時(shí)校正p值
eavg_1 = rowMeans(et1) #progres是對(duì)照組
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)#檢驗(yàn)結(jié)果
DEG_t.test
head(DEG_t.test)

19.用limma包對(duì)表達(dá)矩陣及樣本分組信息進(jìn)行差異分析韭畸，得到差異分析表格，畫火山圖

#火山圖
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

#帶標(biāo)題蔓搞、上下調(diào)基因的火山圖
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')
)#設(shè)置閾值
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',])
)#圖片標(biāo)題
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.對(duì)T檢驗(yàn)結(jié)果的P值和limma包差異分析的P值畫散點(diǎn)圖胰丁，看看哪些基因相差很大

DEG_t.test#T檢驗(yàn)的差異分析結(jié)果
DEG#limma包分析的差異結(jié)果

拿LSM2這個(gè)基因來說，兩種分析結(jié)果里面的p adjust差距還是很大的喂分。
可參考：http://www.reibang.com/p/84bfdf91118c
（寫到這兒有些累了锦庸，以后理解了再更新）