友情提示：一定要有R語言的基礎(chǔ)，建議學習生信技能樹的《生信人這樣這樣學R語言》之后（最好再把中級20題做完），再來學習本系列視頻！
總的來說，理解了R語言中級20題之后，這一塊的內(nèi)容學起來就輕松了很多。
1.通讀文獻閱讀及規(guī)律
這一塊兒主要就是在你感興趣的文章里，找到其測序信息的GSE號，在GEO數(shù)據(jù)庫中直接檢索該文件號。
2.了解GEO數(shù)據(jù)庫
比如我們對GSE42872 這個測序結(jié)果比較感興趣，我們現(xiàn)在GEO Database中搜索該文件號，點擊進入。
我們可以通過原文章了解實驗設計及摘要，也可以在GEO Database中瀏覽。原文章是比較簡單的一個實驗設計思路，直接用藥物vemurafenib或空白對照處理了黑色素瘤細胞后去做RNA seq，每組3個對照。
再向下看，就是測序信息了，這里要說明三個號：
數(shù)據(jù)集：GSE(GEO series)，也就是本次查詢的號，是一個實驗設計下的測序結(jié)果的集合，一個數(shù)據(jù)集下有多個樣本，一篇文章也可能會出現(xiàn)多個數(shù)據(jù)集。
平臺：GPL(GEO platform)，測序公司所提供的的測序系統(tǒng)，在這次實驗中運用的是 Affymetrix Human Gene 1.0 ST Array。
樣本：GSM(GEO sample)，每個樣本的測序結(jié)果都會有一個GSM號。在這個實驗中一共有6個樣本，3個實驗組，3個對照組。
3.數(shù)據(jù)下載的3種方法
①RAW data處下載，不推薦，不好分析。舊版本芯片用readaffy包來處理，illumina用iumiR.batch來處理。
②下載matrix表達矩陣，可直接導入R。（我比較喜歡這種方法）
gset=read.table('Downloads/GSE42872_series_matrix.txt.gz',sep='\t',comment.char='!',header = T)#讀入，txt文件的分隔符為\t，前面會有一些以！開頭的描述信息，刪去這些，表頭要有
class(gset)
str(gset)
rownames(gset)=gset[,1]#原本的第一列為探針名，現(xiàn)在將行名改為探針名
gset=gset[,-1]#刪去原來的第一列
③在R中直接用GSE號下載：
source("[http://www.bioconductor.org/biocLite.R)](http://www.bioconductor.org/biocLite.R))
biocLite("GEOquery")
library(GEOquery)
gset<-getGEO("GEO number",GSEMatrix=TURE,AnnotGPL=TURE，destdir=".")
需要注意，如果是第二種方法下載，導入的直接是dataframe，第三種方法下載的是對象，里面還包含了描述信息，但是getGEO會幫助提取其中的表達矩陣。
4.ID轉(zhuǎn)換技巧大全
同R中級20題中的一樣，不過需要注意的是該套測序的探針用的是哪套注釋信息，大家可以從platforms這里查看：
可以在網(wǎng)上找注釋信息對應合集，如下述鏈接，在其中找GPL6244 [HuGene-1_0-st]對應的R包是hugene10sttranscriptcluster，下載該注釋包。
https://blog.csdn.net/weixin_40739969/article/details/103186027
#探針注釋信息下載
BiocManager::install('hugene10sttranscriptcluster.db') #下載，注意在包的后面手動加上.db
library('hugene10sttranscriptcluster.db')#載入
ls("package:hugene10sttranscriptcluster.db")#查看，找到后面是SYMBOL的，復制名字
ids=toTable(hugene10sttranscriptclusterSYMBOL)#提取注釋信息
table(table(ids$symbol)>1)#統(tǒng)計有多少個基因?qū)粋€以上的探針
table(sort(table(ids$symbol)))#統(tǒng)計有1、2、3等個探針數(shù)的基因分別有多少個，大多數(shù)基因只有一個探針
下載好之后，就開始將探針名更改為對應的基因名。先判斷探針有無對應的基因名，然后對有多個探針的基因保留平均數(shù)或中位數(shù)或其他指標作為該基因的表達量。
#更改ID
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)計T、F的個數(shù)，有13483個探針不存在對應的基因名
n_exprSet=exprSet[!(rownames(exprSet)%in% ids$probe_id)==F,]#將包含在探針包里的探針測序數(shù)據(jù)新建一個數(shù)據(jù)框，這里很奇怪，因該是滿足條件的包括進來，這里卻填F的時候會包括進來
tmp = by(n_exprSet,ids$symbol,
function(x)
rownames(x)[which.max(rowMeans(x))] )#by函數(shù)，將對象按行或某種方式分為一個個小的子集，對每個子集進行操作。這里將對n_exprSet的行按與ids中的symbol的對應關(guān)系，將同一個symbol的探針放在一起，然后對子集中的每一行取平均值，返回每個子集中平均值大的那一行的行名（探針名）
probes = as.character(tmp)#定義為字符型
View(probes)#這里就是樣本中相同基因平均值最大的探針
exprSet=exprSet[rownames(exprSet) %in% probes ,]#對exprSet的行進行操作，若該行的行名在probes中出現(xiàn)了則保留下來
dim(exprSet)#現(xiàn)在為18821個探針在6個樣本中的表達量
View(exprSet)
rownames(exprSet)=ids[match(rownames(exprSet),ids$probe_id),2]#將exprSet的行名與ids的probe_id相互匹配，若相同，則exprSet的行名為ids的第二列
View(exprSet)
5.了解你的表達矩陣
通過一些簡單的分析查看自己的表達矩陣是否合乎邏輯、前面的處理有無錯誤。
①查看管家基因的表達量，和表達量的中位數(shù)（boxplot）作比較，管家基因表達量一般比較高。
GAPDH的表達量均在14以上，β-ACTIN的表達量均在13以上，而表達量的中位數(shù)是在8到9之間，關(guān)鍵基因的表達量在所有表達量中處于較高的水平，這一項是合格的。
exprSet['GAPDH',]#查看GAPDH在六個樣本中的表達量
exprSet['ACTB',]#查看β-actin在六個樣本中的表達量
boxplot(exprSet)#畫六個樣本表達量的箱型圖
②看表達矩陣分布圖，查看所有樣本表達量的中位數(shù)，應該是處于差不多水平的才可以進行比較。
六個樣本的表達量中位數(shù)基本處于同一水平線上。
③畫hclust圖，查看同一分組是否在一起。
首先處理數(shù)據(jù)：
#添加樣本分組信息
group_list=c(rep('control',3),rep('case',3))#新建list，根據(jù)GEO描述得知前三個為對照組、后3個為實驗組
group_list
#將寬數(shù)據(jù)變?yōu)殚L數(shù)據(jù)
library(reshape2)#使用該包中的melt函數(shù)
a=rownames(exprSet)#這里需要操作一波，在實踐中出現(xiàn)的，我單獨寫一篇專門講述melt函數(shù)
a=as.data.frame(a)
exprSet=cbind(a,exprSet)
exprSet=exprSet[,-1]
exprSet1=cbind(a,exprSet)
exprSet_L=melt(exprSet1)#處理數(shù)據(jù)，將多維數(shù)據(jù)轉(zhuǎn)變?yōu)橐痪S數(shù)據(jù)，即寬數(shù)據(jù)變?yōu)殚L數(shù)據(jù)
colnames(exprSet_L)=c('probe','sample','value')#每一列依次為探針（基因）、樣本來源、表達量
exprSet_L$group=rep(group_list,each=nrow(exprSet))
head(exprSet_L)
聚類分析，相同處理的樣本應該是聚在一起的：
#聚類分析
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)
④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)
6.差異分析（limma包），火山圖、熱圖的繪制及美化
limma包分析時需要提前準備好三個矩陣：表達矩陣、分組矩陣、比較矩陣。表達矩陣已經(jīng)處理好了，下面只需制作分組矩陣（用0、1表示實驗組和對照組的樣本）和比較矩陣（誰和誰比，實驗組比上對照組）。
library(limma)#就是給表達矩陣、分組矩陣、比較矩陣，limma給你分析差異結(jié)果
design=model.matrix(~0+factor(group_list))#分組矩陣
design#分組矩陣制作好了，用0和1代替
colnames(design)=levels(factor(group_list))#給分組矩陣的列名用分組信息來表示
rownames(design)=colnames(exprSet)#行名用樣本名來表示
contrast.matrix=makeContrasts(case-control,levels = design)
contrast.matrix #比較矩陣，makeContrasts用來告訴誰比誰
fit=lmFit(exprSet,design)#lmFit函數(shù)為線性模型擬合，給表達矩陣和分組矩陣，輸出
fit2 <- contrasts.fit(fit, contrast.matrix) #把上一步的線性擬合矩陣矩陣和比較矩陣給它，計算差異和標準誤
fit2 <- eBayes(fit2)#通過eBayes將標準誤調(diào)整到一個相同的水平再來計算t、logFC等值
tempOutput = topTable(fit2, coef=1, n=Inf)#用topTable提取
nrDEG = na.omit(tempOutput) #該函數(shù)是刪去缺少某些值的對象
head(nrDEG)
write.csv(nrDEG,"limma_notrend.results.csv",quote = F)#將差異分析結(jié)果保存到本地
檢查一下差異分析的結(jié)果：
exprSet['CD36',]#查看表達情況
#熱圖畫一下前25個基因，檢驗是否有差異
choose_gene=head(rownames(nrDEG),25)
choose_matrix=exprSet[choose_gene,]
choose_matrix=t(scale(t(choose_matrix)))
pheatmap::pheatmap(choose_matrix)#可以看到這些基因的表達量在試驗和對照組中差異還是很大的
畫火山圖（和R20題的操作一樣）：
#帶標題、上下調(diào)基因的火山圖
logFC_cutoff <- with(nrDEG,mean(abs( logFC)) + 2*sd(abs( logFC)) )
nrDEG$change = as.factor(ifelse(nrDEG$P.Value < 0.05 & abs(nrDEG$logFC) > logFC_cutoff,
ifelse(nrDEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
)#設置閾值
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
'\nThe number of up gene is ',nrow(nrDEG[nrDEG$change =='UP',]) ,
'\nThe number of down gene is ',nrow(nrDEG[nrDEG$change =='DOWN',])
)#圖片標題
library(ggplot2)
g = ggplot(data=nrDEG, 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)
7.基因集的富集分析（注釋）KEGG
富集分析或注釋等通俗來講就是這些差異基因處于哪些通路或者是功能中，即造成了哪些通路、功能發(fā)生改變。這里用到的是超幾何分布的原理，簡單介紹一下：假設一共有20000個基因，差異基因有1000個，也就是說每個基因被抽中的概率為1/20?，F(xiàn)在有一個通路有100個基因，若隨機抽，則差異基因中屬于該通路的有5個，但這1000個基因中有20個基因都是屬于這個通路，這就說明該通路是有問題的。
在用clusterProfiler包進行分析時，輸入的ID類型為ENTREZID，但之前我們用的是SYMBOL，需要轉(zhuǎn)換一下。ENTREZID即gene id。
#ID轉(zhuǎn)換
gene=head(rownames(nrDEG),1000)#我們在這里只取前1000個差異基因來做，具體數(shù)字看自己的需求
library(clusterProfiler)#之前的gene ID是symbol，現(xiàn)在要用ENTREZID，所以需要進行轉(zhuǎn)換
library(org.Hs.eg.db)
gene.df <- bitr(gene, fromType = "SYMBOL",
toType = c("ENSEMBL","ENTREZID"),
OrgDb = org.Hs.eg.db)#從SYMBOL轉(zhuǎn)換為ENTREZID和ENSEMBL，可根據(jù)自己的需要轉(zhuǎn)換
head(gene.df)#檢查幾個看看做得對不對
#分析
kk <- enrichKEGG(gene = gene.df$ENTREZID,
organism = 'hsa',
pvalueCutoff = 0.05)
head(kk)[,1:6]
如第一個通路是和DNA復制相關(guān)的，這條通路中一共有36個基因，有18個基因都出現(xiàn)在差異基因中了，說明該項實驗中，藥物處理后DNA的復制發(fā)生了很大很大的變化。
超幾何分布檢驗的弊端就是只知道哪些通路發(fā)生了改變，但不知道是哪些基因在其中起到了比較重要的作用，這個時候可以用GSE進行分析（后續(xù)單獨寫一篇）。
