轉(zhuǎn)錄組數(shù)據(jù)繞不過(guò)差異分析。
為什么選擇這兩個(gè)包呢？
DEseq2針對(duì)有生物學(xué)重復(fù)的樣本。（一般情況下應(yīng)該是都需要生物學(xué)重復(fù)的）
edgeR對(duì)于單個(gè)樣本是比較好的。（但是細(xì)胞材料真的沒(méi)辦法，細(xì)胞是又貴又難養(yǎng)。下一篇介紹。）
一、DEseq2差異表達(dá)分析
1、安裝
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("DESeq2")
library(DESeq2)

2、準(zhǔn)備數(shù)據(jù)
featureCounts定量后的數(shù)據(jù)，或者FPKM數(shù)據(jù)（下一遍講如何獲取FPKM）
定量后的數(shù)據(jù)的數(shù)據(jù)格式如下：

image.png

colData，其實(shí)就是表型數(shù)據(jù)。
格式如下

image.png

3、計(jì)算差異基因

library(tidyverse)
library(DESeq2)
#載入定量后的數(shù)據(jù)
setwd("E:/4/")
mycounts <- read.table("07counts_matrix.txt", header=TRUE)

#如果重復(fù)之間需要就平均值，那就計(jì)算下均值。
mycounts$Control <- round((mycounts$control2.sorted.bam + mycounts$control3.sorted.bam)/2)
mycounts$A <- round((mycounts$A1.sorted.bam + mycounts$A2.sorted.bam + mycounts$A3.sorted.bam)/3)
mycounts$B <- round((mycounts$B1.sorted.bam + mycounts$B2.sorted.bam + mycounts$B3.sorted.bam)/3)
mycounts <- subset(mycounts, select=c(Geneid,Control,A,B))

#前面我輸入的時(shí)候沒(méi)有把ensemble名字修改成symble，所以這里就不需要。
#如果處理成名字之后的數(shù)據(jù)會(huì)存在重復(fù)，所以需要把重復(fù)的刪掉。
rows <- rownames(unique(mycounts['Geneid']))
mycounts <- mycounts[rows,]

#這里有個(gè)Geneid，需要去除，先把第一列當(dāng)作行名來(lái)處理
rownames(mycounts)<-mycounts[,1]
#把帶Geneid的列刪除
mycounts <- mycounts[,-1]
#head(mycounts)

#載入colData文件
colData <-read.csv("08counts_matrix_condition.csv")
#colData <- colData[1:5,]
#colData <- colData[c(1,2,6,7,8),]
#構(gòu)建數(shù)據(jù)矩陣
dds <- DESeqDataSetFromMatrix(mycounts, colData, design= ~ condition)
dds <- DESeq(dds)

#接下來(lái)，我們要查看treatmentt VS control的總體結(jié)果，并根據(jù)p-value進(jìn)行重新排序。
#利用summary命令統(tǒng)計(jì)顯示一共多少個(gè)genes上調(diào)和下調(diào)（FDR0.1）
res <- results(dds, contrast=c("condition", "control", "treatment"))##或者res= results(dds)
res <- res[order(res$pvalue),]
summary(res)
#所有結(jié)果先進(jìn)行輸出
write.csv(res,file="DESeq2_results.csv")

#獲取padj（p值經(jīng)過(guò)多重校驗(yàn)校正后的值）小于0.05，表達(dá)倍數(shù)取以2為對(duì)數(shù)后大于1.5或者小于-1.5的差異表達(dá)基因。
#DS和PD都是合理的，所以選擇1.5。
diff_gene_deseq2 <-subset(res, padj < 0.05 & abs(log2FoldChange) > 1.5) 
##或者diff_gene_deseq2 <-subset(res,padj < 0.05 & (log2FoldChange > 1.5 | log2FoldChange < -1.5))
dim(diff_gene_deseq2)
write.csv(diff_gene_deseq2,file= "DESeq2_diffExpression.csv")

####提取出你所需要的log2FoldChange和pvalue列。
diff_gene_deseq2 <- na.omit(diff_gene_deseq2)
###這里把需要的兩行提取出來(lái)先，是log2FoldChange和padj
nrDEG <- diff_gene_deseq2[,c(2,6)]
View(nrDEG)
###然后換個(gè)名字即可。
colnames(nrDEG) <- c('log2FoldChange','pvalue')
View( nrDEG)
write.csv(nrDEG,file= "DESeq2_diffExpression_logFC.csv")

4、數(shù)據(jù)可視化
4.1、火山圖
加載包

library(ggplot2)
library(ggrepel)
library(export)

上步驟得到了差異基因，賦值給一個(gè)新的參數(shù)。
之所以這么干，是因?yàn)槲疫@邊是兩個(gè)分開(kāi)寫(xiě)的腳本。

DEG <- res
dim(DEG)
#這里用的是DESeq2的DEG，刪掉NA。
DEG <- na.omit(DEG)
dim(DEG)
# 使用基礎(chǔ)函數(shù)plot繪圖
plot(DEG$log2FoldChange,-log2(DEG$padj))
# 確定差異表達(dá)倍數(shù),abs表示絕對(duì)值
logFC_cutoff <- with(DEG,mean(abs(log2FoldChange)) + 2*sd(abs(log2FoldChange)))
# 取前兩位小數(shù)
logFC_cutoff <- round(logFC_cutoff, 2)
logFC_cutoff
查看一下logFC_cutoff的值。

# 確定上下調(diào)表達(dá)基因。
#方法1：按照l(shuí)ogFC_cutoff繪圖
#方法2：按照差異基因的log2FoldChange大于1.5或者其他值繪圖。
DEG$change = as.factor(ifelse(DEG$padj < 0.05 & abs(DEG$log2FoldChange) > logFC_cutoff,
                        ifelse(DEG$log2FoldChange > logFC_cutoff ,'UP','DOWN'),'STABLE'))
                        
DEG$change = as.factor(ifelse(DEG$padj < 0.05 & abs(DEG$log2FoldChange) > 1.5,
                        ifelse(DEG$log2FoldChange > 1.5 ,'UP','DOWN'),'STABLE'))
 
##確定需要在火山圖上展示的字。
#this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,2),
                    '\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
                    '\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',]))
                    
#this_tile <- paste0('Cutoff for logFC is ',round(2.0,1),
                    '\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
                    '\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',]))                    

###繪圖的時(shí)候必須是dataframe，所以需要轉(zhuǎn)換一下。
DEG <- data.frame(DEG)
## 顏色與分組一一對(duì)應(yīng)
g <- ggplot(data=DEG, aes(x=log2FoldChange, y=-log10(padj),color=change)) + 
  geom_point(shape = 16, size=2) + 
  theme_set(theme_set(theme_bw(base_size=20))) + 
  xlab("log2 fold change") + 
  ylab("-log10 p-value") +
  ggtitle( this_tile ) + ##這個(gè)可以不要
  theme(plot.title = element_text(size=15,hjust = 0.5)) + 
  theme_classic()+
  scale_colour_manual(values = c('blue','black','red'))
graph2ppt(file="Volcano_Plot.ppt", width=10, aspectr=1)

4.2、熱圖
載入包

library(pheatmap)

載入數(shù)據(jù)

dat <- mycounts
###載入表型數(shù)據(jù)
View(colData$condition)
# 提取差異倍數(shù)
# 根據(jù)需要修改DEG的值
FC <- DEG$log2FoldChange
###根據(jù)FC篩選。
#FC <- subset(FC, abs(FC)>1)
#View(FC)
names(FC) <- rownames(DEG)
#View(FC)

# 排序差異倍數(shù)，提取前100和后100的基因名
#或者全部的基因名字
DEG_200 <- c(names(head(sort(FC),100)),names(tail(sort(FC),100)))
DEG_all <- c(names(sort(FC)))
###將篩選得到的基因，提取并進(jìn)行標(biāo)準(zhǔn)化。
dat <- t(scale(t(dat[DEG_200,])))
dat <- t(scale(t(dat[DEG_all,])))
###針對(duì)不同數(shù)據(jù)集，求平均值
dat.fram <- data.frame(dat)
dat.fram$Control <- (dat.fram$control2.sorted.bam + dat.fram$control3.sorted.bam)/2
dat.fram$A <- (dat.fram$A1.sorted.bam + dat.fram$A2.sorted.bam + dat.fram$A3.sorted.bam)/3
dat.fram$B <- (dat.fram$B1.sorted.bam + dat.fram$B2.sorted.bam + dat.fram$B3.sorted.bam)/3
dat2 <- subset(dat.fram, select=c(Control,A,B))
#這里需要重新導(dǎo)入表型數(shù)據(jù)的名字。
colData <-read.csv("09counts_matrix_condition_mean.csv")
colData <- colData[1:2,]
colData <- colData[c(1,3),]
group <- data.frame(colData$condition)
rownames(group)=colnames(dat2)
#繪圖
#show_colnames =T，show_rownames = F，cluster_cols = T這三個(gè)參數(shù)根據(jù)需要設(shè)置
#顏色根據(jù)需要設(shè)置colorRampPalette(c("green", "balck", "red"))
g2 <- pheatmap(dat2,show_colnames =T,show_rownames = F, cluster_cols = T,
         annotation_col=group,border=FALSE,
         color = colorRampPalette(c("navy", "white", "firebrick3"))(50))     
graph2ppt(file="pheatmap_cluster_mean.ppt", width=10, aspectr=1)