DESeq2進(jìn)行差異分析

### 安裝包
############################
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("DESeq2")
install.packages("corrplot")
install.packages("PerformanceAnalytics")
install.packages("factoextra")
############################

1 ###數(shù)據(jù)讀入和處理

rm(list=ls()) 
countdata<- read.table("all_fcount.matrix.txt", header=TRUE,row.names = 1) #導(dǎo)入數(shù)據(jù)
head(countdata) # 查看數(shù)據(jù)前幾行(列名 太長(zhǎng)自己展示看)
##  過濾featurecounts后，每個(gè)樣本計(jì)數(shù)小于等于10的?。。。。。?！
countdata.filter<-countdata[rowSums(countdata)>=0&apply(countdata,1,function(x){all(x>=10)}),]
head(countdata.filter) 

colnames(countdata.filter) <- c("RIBEN_1","RIBEN_2","RIBEN_3","CHINA_1","CHINA_2","CHINA_3") # 修改列名
head(countdata.filter)
dim(countdata.filter) # 查看數(shù)據(jù)維度，即幾行幾列

2 #### dds矩陣的創(chuàng)建 ####

library(DESeq2)
condition <- factor(c("RIBEN","RIBEN","RIBEN","CHINA","CHINA","CHINA")) # 賦值因子，即變量
coldata <- data.frame(row.names=colnames(countdata.filter), condition) # 創(chuàng)建一個(gè)condition數(shù)據(jù)框
dds <- DESeqDataSetFromMatrix(countData=countdata.filter, colData=coldata, design=~condition) ##構(gòu)建dds矩陣(后面很多分析都是基于這個(gè)dds矩陣)

3 ### DESeq2進(jìn)行差異分析 #####

dds <- DESeq(dds)
resdata<- results(dds,contrast = c("condition","RIBEN","CHINA"))##此為前比后
table(resdata$padj<0.05 ) # Benjamini-Hochberg矯正后的p<0.05的基因數(shù)?。。。。。?！
res_padj <- resdata[order(resdata$padj), ]  ##按照padj(矯正后的p值)列的值排序

write.table(res_padj,"diffexpr_padj_results.txt",,quote = F,sep = '\t')#### 將結(jié)果文件保存到本地，打開在第一列頭加gene

4 #### DESeq2分析中得到的resdata進(jìn)行繪制火山圖 #####

rm(list=ls())  
resdata <- read.table("diffexpr_padj_results.txt",header = T,sep = '\t',row.names = 1)### 加載DESeq2中生成的resdata文件
resdata$label <- c(rownames(resdata)[1:10] ,rep(NA,(nrow(resdata)-10)))##對(duì)前10個(gè)基因進(jìn)行標(biāo)注

library(ggplot2)
ggplot(resdata,aes(x=log2FoldChange,y=-log10(padj))) +
geom_hline(yintercept = -log10(0.05),linetype = "dashed",color = "#999999")+ ##橫向水平參考線，顯著性--p值
geom_vline(xintercept = c(-1 , 1),linetype = "dashed", color = "#999999")+ ##縱向垂直參考線，差異倍數(shù)--foldchange
geom_point(aes(size = -log10(padj),color = -log10(padj))) + ##散點(diǎn)圖
scale_color_gradientn(values = seq(0,1,0.2),
                       colors = c("#39489f","#39bbec","#f9ed36","#f38466","#b81f25"))+ ##指定顏色漸變模式
scale_size_continuous(range = c(0,3))+##指定散點(diǎn)漸變模式
  
##主題調(diào)整
###  theme_grey()為默認(rèn)主題，theme_bw()為白色背景主題，theme_classic()為經(jīng)典主題
theme_bw()+
##調(diào)整主題和圖例位置
theme(panel.grid.major = element_blank(), #刪除主網(wǎng)格線
        panel.grid.minor = element_blank(), #刪除次網(wǎng)格線
        legend.position = c(0.08,0.9),      #設(shè)置圖例位置
        legend.justification = c(0,1))+
##設(shè)置部分圖例不顯示
guides(col = guide_colourbar(title = "-Log10_q-value"),
         size = "none")+
##添加標(biāo)簽
geom_text(aes(label=c(label),color = -log10(padj)), size = 3, vjust = 1.5, hjust = 1)+
##修改坐標(biāo)軸
xlab("Log2FC")+ylab("-Log10(FDR q-value)") 
##保存圖片
ggsave("vocanol_plot.pdf", height = 9 , width = 10)

5 篩選差異基因

subset(resdata,pvalue < 0.05) -> diff  ## 先篩選pvalue < 0.05的行?。。。?！
subset(diff,log2FoldChange < -0.585) -> down  ## 篩選出下調(diào)的,此處定位1.5倍
subset(diff,log2FoldChange > 0.58) -> up  ## 篩選出上調(diào)的
dim(down) # 查看數(shù)據(jù)維度，即幾行幾列
dim(up)
up1=as.data.frame(up )
#導(dǎo)出上、下調(diào)基因的那一列
up_names<-rownames(up)
write.table(up_names,'up_gene.txt',quote = F,sep = '\t',row.names = F)
down_names <- rownames(down)
write.table(down_names,'down_gene.txt',quote = F,sep = '\t',row.names = F)

6 表達(dá)矩陣探索,選取差異表達(dá)的基因做熱圖 (注意是用矯正后p值排外序的)

library(pheatmap)
choose_gene=head(rownames(resdata),50)  
choose_matrix=countdata.filter[choose_gene,]  #抽取差異表達(dá)顯著的前50個(gè)基因
choose_matrix=t(scale(t(choose_matrix)))  #用t函數(shù)轉(zhuǎn)置，scale函數(shù)標(biāo)準(zhǔn)化

pheatmap(choose_matrix)

7 PCA分析 使用對(duì)dds矩陣處理后的vst或rld值

dds_pca <- estimateSizeFactors(dds) #計(jì)算每個(gè)樣本的歸一化系數(shù)
vsd <- vst(dds_pca,blind = FALSE)  ## 方差穩(wěn)定變換

rld <- rlog(dds_pca,blind = FALSE)  # 正則化對(duì)數(shù)變換
# rlog函數(shù)將count data轉(zhuǎn)換為log2尺度，以最小化有small counts的行的樣本間差異，并使library size標(biāo)準(zhǔn)化
## vst函數(shù)快速估計(jì)離散趨勢(shì)并應(yīng)用方差穩(wěn)定變換。
## 數(shù)據(jù)集小于30個(gè)樣品可以用rlog，數(shù)據(jù)集大于30個(gè)樣品用vst，因?yàn)閞log速度慢
# 轉(zhuǎn)換的目的是，為了確保所有基因有大致相同的貢獻(xiàn)
plotPCA(vsd, intgroup=c("condition"))
plotPCA(rld, intgroup=c("condition"))