基于上游分析獲得表達矩陣后弹澎，就可以進行差異分析了欢搜，最基礎(chǔ)的莫過于利用DEseq2包尋找差異基因恕曲。學習自b站視頻生信技能樹轉(zhuǎn)錄組視頻~
0、數(shù)據(jù)準備及前處理
- 由于目前在家不方便獲得數(shù)據(jù)池户，因此參照教學視頻使用airway包的現(xiàn)成數(shù)據(jù)。
如鏈接中官方介紹，airway包背景是使用地塞米松處理的四個人氣道（支氣管校焦？）平滑肌細胞（即共8個數(shù)據(jù)）進行rna-seq實驗赊抖，分析得到的基因表達矩陣。
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("airway")
library(airway)
data(airway)
airway
class(airway)
- 如上返回結(jié)果寨典，airway為RangedSummarizedExperiment類型的對象氛雪，是將數(shù)據(jù)歸檔R包的一個好方法，使用時主要會涉及到兩個函數(shù)操作
assay(data)表示取觀測為基因名凝赛，變量為樣本信息的表達矩陣注暗；
colData(data)則是取記錄詳細的樣本信息；
exp=assay(airway) #取原始表達矩陣數(shù)據(jù)
tmp=colData(airway)
tmp=as.data.frame(airway)
#由于后續(xù)分析需要分組信息墓猎，因此這里我們?nèi)∠碌谌行畔?group=colData(airway)[,3]
class(group)
值得注意的是
colData(airway)得到的結(jié)果并非直接為真正的表格捆昏，還需要用as.data.frame()函數(shù)轉(zhuǎn)換下，后面分析還會遇到類似的問題毙沾。
index1=c(rowSums(exp)>0)
class(index1)
exp2=exp[index1,]
- 此處操作的目的在于將沒有read比對上的基因行刪去骗卜，簡潔數(shù)據(jù)~
1、DEseq2包差異分析
（1）設(shè)置單列分組樣本信息表
group
sample=c(colnames(exp2)) #取樣本名信息
f=as.data.frame(group)
rownames(f)=c(sample) #設(shè)置行名為樣本名信息
colnames(f)='condition' #列名
（2）DESeq2分析
library("DESeq2")
dds=DESeqDataSetFromMatrix(countData = exp2,
colData = f,
design = ~ condition)
dep=DESeq(dds)
res=results(dep,contrast=c("condition","trt","untrt"))
這里要注意一下左胞，所謂差異分析寇仓，即比較實驗組與對照組相比，基因表達變化情況烤宙。通過
contrast=參數(shù)分別設(shè)置樣本分組信息表的列名遍烦，實驗組名（numerator，分子）躺枕，對照組名（denominator服猪，分母）；即研究分子/分母的結(jié)果
DEG1=as.data.frame(res)
write.csv(x=DEG1,file="res.csv")
- 返回的六列值中重點關(guān)注 log2FoldChange,pvalue,padj三列內(nèi)容拐云。
log2FoldChange中的FoldChange即倍數(shù)變化的意思罢猪。某一基因表達 實驗組/對照組，可分為兩種結(jié)果：大于1叉瘩，即上調(diào)膳帕；小于1，即下調(diào)薇缅。利用log函數(shù)將其分別轉(zhuǎn)化成正數(shù)危彩、負數(shù)，更便于理解捅暴。（通常為了防止取log2時產(chǎn)生NA恬砂，我們會給表達值加1（或者一個極小的數(shù)），也就是log2(B+1) - log2(A+1)）蓬痒；
pvalue是假設(shè)檢驗的p值泻骤，值越衅岣帷（一般標準0.05）即說明log2FoldChange數(shù)據(jù)的可靠性越強；padj為adjusted p-value值狱掂，也稱Q-value演痒、FDR（false discovery rate）為效驗過后的p-value值，更有信服力趋惨。
resOrdered=res[order(res$pvalue),] #order函數(shù)是給出從小到大排序后的位置(默認升序)
sum(res$padj<0.1, na.rm=TRUE) # na.rm 為remove na值鸟顺，否則會影響統(tǒng)計量
sum(res$pvalue<0.05, na.rm=TRUE)
sum(!is.na(res$pvalue))
resSig=subset(resOrdered, padj<0.1)
DEG2=as.data.frame(resSig)
write.csv(x=resSig, file = "results.csv")
2、繪制熱圖
- 根據(jù)基因表達矩陣比對上的read數(shù)目繪制
- 這里選取最有意義的50個基因比對read數(shù)目為例
library(pheatmap)
choose_gene=head(rownames(DEG2),50)
#根據(jù)DEG2結(jié)果器虾，提取目標基因名
choose_matrix=exp2[choose_gene,]
choose_matrix=t(scale(t(choose_matrix)))
一般繪制熱圖讯嫂，都需要對數(shù)據(jù)進行歸一化處理。這里用到了
t()轉(zhuǎn)置函數(shù)的原因是因為基因表達矩陣表中觀測為基因兆沙，變量為樣本欧芽。而scale()歸一化函數(shù)是針對每一列單獨處理的。因為每一個基因的本身情況各不相同葛圃，而且我們的目的是基于實驗組與對照組的比較來研究基因表達變化千扔，所以需要轉(zhuǎn)置下表達矩陣（觀測為樣本，變量為基因）库正，然后再進行對每一列的歸一化處理曲楚；最后再進行一次轉(zhuǎn)置，即得到下圖結(jié)果褥符。
pheatmap(choose_matrix,filename = 'DEG_top50_heatmap.png')
3龙誊、繪制火山圖
- 本質(zhì)上就是根據(jù)DEG數(shù)據(jù)里基因的log2FoldChange與pvalue繪制點圖
- 把點分成三類：上調(diào)基因、下調(diào)基因與無顯著改變基因喷楣。
DEG1=na.omit(DEG1)
(1) 設(shè)定分類判斷閾值并增加分類列
logFC_cutoff=with(DEG1,mean(abs(log2FoldChange))+2*sd(abs(log2FoldChange)))
logFC_cutoff
2^logFC_cutoff
# 返回結(jié)果為 3.67
DEG1$change=as.factor(ifelse(DEG1$pvalue<0.05 & abs(DEG1$log2FoldChange)>logFC_cutoff,
ifelse(DEG1$log2FoldChange>logFC_cutoff,'UP','DOWN'),'NOT'))
(2)編寫圖中的注釋信息
this_tile=paste('Cutoff for logFC is ',round(logFC_cutoff,3),
'\nThe number of up gene is ',nrow(DEG1[DEG1$change=='UP',]),
'\nThe number of down gene is ',nrow(DEG1[DEG1$change=='DOWN',]))
- 其中
\n表示換行符
(3)利用ggplot2繪制火山圖
library(ggplot2)
g=ggplot(data=DEG1,
aes(x=log2FoldChange,y=-log10(pvalue), #這里將pvalue取負對數(shù)
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 pvalue") + #軸標簽
ggtitle(this_tile)+theme(plot.title=element_text(size=15,hjust=0.5)) +
scale_color_manual(values=c('blue','black','red')) #設(shè)定顏色
ggsave(g,filename='volcano.png')
# 如下圖载迄，一般會重點關(guān)注兩邊的，偏上方的基因點
以上是基于基因表達矩陣進行差異分析以及繪制熱圖抡蛙、火山圖的大致流程。這應(yīng)該下有分析簡單的開始魂迄，后期還有富集分析等粗截，再慢慢學習吧~
