一.樣本表達(dá)分布

箱式圖

##魔鬼操作一鍵清空，以防前面出現(xiàn)的變量名干擾后續(xù)的操作
rm(list = ls())
options(stringsAsFactors = F)

# 加載原始表達(dá)的數(shù)據(jù)
lname <- load(file = "Step01-airwayData.Rdata")
lname

exprSet <- express_cpm
exprSet[1:6,1:6]

# 樣本表達(dá)總體分布-箱式圖
library(ggplot2)
# 構(gòu)造繪圖數(shù)據(jù)
data <- data.frame(expression=c(exprSet),sample=rep(colnames(exprSet),each=nrow(exprSet)))
head(data)

p <- ggplot(data = data,aes(x=sample,y=expression,fill=sample))
p1 <- p + geom_boxplot() + theme(axis.text.x = element_text(angle = 90)) + xlab(NULL) + ylab("log2(CPM+1)")
p1

# 保存圖片
pdf(file = "../Analysis/sample_cor/sample_boxplot.pdf",width = 6,height = 8)
print(p1)
dev.off()

為何即使組與組之間即使有生物學(xué)上的差異牲芋，但是箱式圖中卻沒(méi)有表現(xiàn)出來(lái)呢?這是 基于前提假設(shè)：大多數(shù)情況下發(fā)生差異表達(dá)的基因只是很少的一部分询兴，并且發(fā)生差異表達(dá)的基因不足以改變整體所有基因表達(dá)水平的分布；如果看樣本的總體表達(dá)分布差異非常大的話(huà)熬甫，不是由于生物水平差異導(dǎo)致的，是誤差導(dǎo)致的【比如不同人做實(shí)驗(yàn)蔓罚、機(jī)子椿肩、protocol】

小提琴圖

image.png

# 樣本表達(dá)總體分布-小提琴圖
p2 <- p + geom_violin() + theme(axis.text = element_text(size = 12),axis.text.x = element_text(angle = 90)) + xlab(NULL) + ylab("log2(CPM+1)")
p2

# 保存圖片
pdf(file = "../Analysis/sample_cor/sample_violin.pdf",width = 6,height = 8)
print(p2)
dev.off()  #打開(kāi)畫(huà)板要養(yǎng)成關(guān)掉的習(xí)慣

與箱式圖結(jié)合發(fā)現(xiàn)，其實(shí)位于中位數(shù)的數(shù)量不是最多的

概率密度分布圖

# 樣本表達(dá)總體分布-概率密度分布圖
m <- ggplot(data=data, aes(x=expression))
p3 <- m +  geom_density(aes(fill=sample, colour=sample),alpha = 0.2) + xlab("log2(CPM+1)")
p3

# 保存圖片
pdf(file = "../Analysis/sample_cor/sample_density.pdf",width = 7,height = 8)
print(p3)
dev.off()

不僅可以看單個(gè)樣本的表達(dá)趨勢(shì)豺谈，還可以看多個(gè)樣本的表達(dá)趨勢(shì)郑象；可以判斷是否有異常的樣本；圖上不同樣本的表達(dá)性就非常一致

以上三個(gè)圖結(jié)合起來(lái)可以看總體樣本表達(dá)分布圖

二.樣本的相關(guān)性

• 層次聚類(lèi)樹(shù)————通過(guò)計(jì)算點(diǎn)與點(diǎn)之間的空間距離對(duì)樣本進(jìn)行類(lèi)別劃分

主要看縱軸Height茬末，也叫刻畫(huà)樣本之間的相似性厂榛；先從下到上看，每個(gè)樣本先是單獨(dú)的一類(lèi)团南，距離相近的一類(lèi)劃分在一起

# 魔幻操作噪沙，一鍵清空
rm(list = ls())  
options(stringsAsFactors = F)

# 加載數(shù)據(jù)并檢查
lname <- load(file = '../Analysis/data/Step01-airwayData.Rdata')
lname

dat <- express_cpm
dat[1:6,1:6]
dim(dat)


## 1.樣本之間的相關(guān)性-層次聚類(lèi)樹(shù)----
sampleTree <- hclust(dist(t(dat)), method = "average")  ##dist是計(jì)算距離【按照行來(lái)算的】，所以要用t來(lái)轉(zhuǎn)置 吐根，把距離矩陣算完之后正歼，開(kāi)始聚類(lèi)了使用的是：hclust函數(shù)
plot(sampleTree)

• PCA主成分分析————提取樣本的綜合特征，即主成分（第一主成分拷橘、第二主成分）對(duì)樣本進(jìn)行分類(lèi)

image.png

## 2.樣本之間的相關(guān)性-PCA----
# 第一步局义，數(shù)據(jù)預(yù)處理
dat <- as.data.frame(t(dat))  ##轉(zhuǎn)置
dat$group_list <- group_list  ##添加group信息

# 第二步喜爷，繪制PCA圖
library(FactoMineR)
library(factoextra)

# 畫(huà)圖僅需要數(shù)值型數(shù)據(jù)，去掉最后一列的分組信息
dat_pca <- PCA(dat[,-ncol(dat)], graph = FALSE)
class(dat_pca)

p <- fviz_pca_ind(dat_pca,
                  geom.ind = "text", # 只顯示點(diǎn)萄唇，不顯示文字
                  col.ind = dat$group_list, # 用不同顏色表示分組
                  palette = c("#00AFBB", "#E7B800"),
                  addEllipses = T, # 是否圈起來(lái)
                  legend.title = "Groups")
p

• 相關(guān)性分析

image.png

# 使用500個(gè)基因的表達(dá)量來(lái)做相關(guān)性圖
library(corrplot)
dim(exprSet)

# 計(jì)算相關(guān)性
M <- cor(exprSet)  #M是一個(gè)對(duì)稱(chēng)矩陣
g <- corrplot(M,order = "AOE",addCoef.col = "white")
##可視化corrplot
corrplot(M,order = "AOE",type="upper",tl.pos = "d")
corrplot(M,add=TRUE, type="lower", method="number",order="AOE",diag=FALSE,tl.pos="n", cl.pos="n")

# 繪制樣本相關(guān)性的熱圖
anno <- data.frame(sampleType=group_list)  ##最重要的地方就是矩陣的行名是數(shù)據(jù)框的列名
rownames(anno) <- colnames(exprSet)
anno
p <- pheatmap::pheatmap(M,display_numbers = T,annotation_col = anno,fontsize = 11,cellheight = 28,cellwidth = 28)
p

pdf(file = "../Analysis/sample_cor/cor.pdf")
print(p)
dev.off()

image.png

image.png

三檩帐、差異分析

三大分析方法：limma、edgeR另萤、DESeq2【原始輸入數(shù)據(jù)都是count值】————用于高通量測(cè)序湃密，不能用于基因芯片

基本理論

image.png

很多指標(biāo)是對(duì)FC取了log，那么logFC>0四敞，B相對(duì)于A就是上調(diào)泛源；反之亦然【對(duì)表達(dá)水平較低的基因敏感，對(duì)表達(dá)水平高的不敏感】

image.png

可以理解為對(duì)P值做的矯正

image.png

limma差異分析

# 清空當(dāng)前對(duì)象
rm(list = ls())
options(stringsAsFactors = F)

# 讀取基因表達(dá)矩陣
lname <- load(file = "Step01-airwayData.Rdata")
lname

exprSet <- filter_count
# 檢查表達(dá)譜
dim(exprSet)
exprSet[1:6,1:6]
table(group_list)

# 加載包
library(limma)
library(edgeR)

## 第一步忿危，創(chuàng)建設(shè)計(jì)矩陣和對(duì)比：假設(shè)數(shù)據(jù)符合正態(tài)分布达箍，構(gòu)建線(xiàn)性模型
# 0代表x線(xiàn)性模型的截距為0
design <- model.matrix(~0+factor(group_list))  ##0是線(xiàn)性模型的截距
colnames(design) <- levels(factor(group_list))
rownames(design) <- colnames(exprSet)
design

兩列八行

# 設(shè)置需要進(jìn)行對(duì)比的分組，需要修改
comp <- 'trt-untrt'   ##不要設(shè)置反了
cont.matrix <- makeContrasts(contrasts=c(comp),levels = d esign)



## 第二步铺厨，進(jìn)行差異表達(dá)分析
# 將表達(dá)矩陣轉(zhuǎn)換為edgeR的DGEList對(duì)象
dge <- DGEList(counts=exprSet)

# 進(jìn)行標(biāo)準(zhǔn)化
dge <- calcNormFactors(dge)   

#Use voom() [15] to convert the read counts to log2-cpm, with associated weights, ready for linear modelling:
v <- voom(dge,design,plot=TRUE, normalize="quantile") 
fit <- lmFit(v, design)
fit2 <- contrasts.fit(fit,cont.matrix)
fit2 <- eBayes(fit2)


## 第三步缎玫，提取過(guò)濾差異分析結(jié)果
tmp <- topTable(fit2, coef=comp, n=Inf,adjust.method="BH")
DEG_limma_voom <- na.omit(tmp)
head(DEG_limma_voom)

# 篩選上下調(diào)，設(shè)定閾值
fc_cutoff <- 2
fdr <- 0.05

DEG_limma_voom$regulated <- "normal"

loc_up <- intersect(which(DEG_limma_voom$logFC>log2(fc_cutoff)),which(DEG_limma_voom$adj.P.Val<fdr))
loc_down <- intersect(which(DEG_limma_voom$logFC< (-log2(fc_cutoff))),which(DEG_limma_voom$adj.P.Val<fdr))

DEG_limma_voom$regulated[loc_up] <- "up"
DEG_limma_voom$regulated[loc_down] <- "down"
  
table(DEG_limma_voom$regulated)

# 添加一列g(shù)ene symbol
library(org.Hs.eg.db)
keytypes(org.Hs.eg.db)

library(clusterProfiler)
id2symbol <- bitr(rownames(DEG_limma_voom), fromType = "ENSEMBL", toType = "SYMBOL", OrgDb = org.Hs.eg.db )
head(id2symbol)

symbol <- rep("NA",time=nrow(DEG_limma_voom))
symbol[match(id2symbol[,1],rownames(DEG_limma_voom))] <- id2symbol[,2]
DEG_limma_voom <- cbind(rownames(DEG_limma_voom),symbol,DEG_limma_voom)
colnames(DEG_limma_voom)[1] <- "GeneID"


# 保存
write.table(DEG_limma_voom,"DEG_limma_voom_all-1.xls",row.names = F,sep="\t",quote = F)

## 取表達(dá)差異倍數(shù)和p值,矯正后的pvalue
DEG_limma_voom <- DEG_limma_voom[,c(1,2,3,6,7,9)]
save(DEG_limma_voom, file = "Step03-limma_voom_nrDEG.Rdata")


## 檢查是否上下調(diào)設(shè)置錯(cuò)了
# 挑選一個(gè)差異表達(dá)基因
head(DEG_limma_voom)

exp <- c(t(express_cpm[match("ENSG00000178695",rownames(express_cpm)),]))
test <- data.frame(value=exp,group=group_list)
library(ggplot2)
ggplot(data=test,aes(x=group,y=value,fill=group)) + geom_boxplot()

edgeR差異分析

rm(list = ls())
options(stringsAsFactors = F)

# 讀取基因表達(dá)矩陣信息
lname <- load(file = "../Analysis/data/Step01-airwayData.Rdata")
lname 

# 查看分組信息和表達(dá)矩陣數(shù)據(jù)
exprSet <- filter_count
dim(exprSet)
exprSet[1:6,1:6]
table(group_list)

# 加載包
library(DESeq2)

# 第一步解滓，構(gòu)建DESeq2的DESeq對(duì)象
colData <- data.frame(row.names=colnames(exprSet),group_list=group_list)
dds <- DESeqDataSetFromMatrix(countData = exprSet,colData = colData,design = ~ group_list)

# 第二步赃磨，進(jìn)行差異表達(dá)分析
dds2 <- DESeq(dds)

# 提取差異分析結(jié)果，trt組對(duì)untrt組的差異分析結(jié)果
tmp <- results(dds2,contrast=c("group_list","trt","untrt"))
DEG_DESeq2 <- as.data.frame(tmp[order(tmp$padj),])
head(DEG_DESeq2)

# 去除差異分析結(jié)果中包含NA值的行
DEG_DESeq2 = na.omit(DEG_DESeq2)

# 篩選上下調(diào)伐蒂，設(shè)定閾值
fc_cutoff <- 2
fdr <- 0.05

DEG_DESeq2$regulated <- "normal"

loc_up <- intersect(which(DEG_DESeq2$log2FoldChange>log2(fc_cutoff)),which(DEG_DESeq2$padj<fdr))
loc_down <- intersect(which(DEG_DESeq2$log2FoldChange< (-log2(fc_cutoff))),which(DEG_DESeq2$padj<fdr))

DEG_DESeq2$regulated[loc_up] <- "up"
DEG_DESeq2$regulated[loc_down] <- "down"

table(DEG_DESeq2$regulated)

# 添加一列g(shù)ene symbol
library(org.Hs.eg.db)
keytypes(org.Hs.eg.db)

library(clusterProfiler)
id2symbol <- bitr(rownames(DEG_DESeq2), fromType = "ENSEMBL", toType = "SYMBOL", OrgDb = org.Hs.eg.db )
head(id2symbol)

symbol <- rep("NA",time=nrow(DEG_DESeq2))
symbol[match(id2symbol[,1],rownames(DEG_DESeq2))] <- id2symbol[,2]
DEG_DESeq2 <- cbind(rownames(DEG_DESeq2),symbol,DEG_DESeq2)
colnames(DEG_DESeq2)[1] <- "GeneID"
head(DEG_DESeq2)

# 保存
write.table(DEG_DESeq2,"../Analysis/deg_analysis/DEG_DESeq2_all.xls",row.names = F,sep="\t",quote = F)

## 取表達(dá)差異倍數(shù)和p值,矯正后的pvalue并保存
colnames(DEG_DESeq2)
DEG_DESeq2 <- DEG_DESeq2[,c(1,2,4,7,8,9)]
save(DEG_DESeq2, file = "../Analysis/deg_analysis/Step03-DESeq2_nrDEG.Rdata")


## 檢查是否上下調(diào)設(shè)置錯(cuò)了
# 挑選一個(gè)差異表達(dá)基因
head(DEG_DESeq2)

exp <- c(t(express_cpm[match("ENSG00000152583",rownames(express_cpm)),]))
test <- data.frame(value=exp,group=group_list)
library(ggplot2)
ggplot(data=test,aes(x=group,y=value,fill=group)) + geom_boxplot()

image.png

差異結(jié)果可視化

image.png

image.png

rm(list = ls())
options(stringsAsFactors = F)

# 加載原始表達(dá)矩陣
load(file = "Step01-airwayData.Rdata")

# 讀取3個(gè)軟件的差異分析結(jié)果
load(file = "Step03-limma_voom_nrDEG.Rdata")
load(file = "Step03-DESeq2_nrDEG.Rdata")
load(file = "Step03-edgeR_nrDEG.Rdata")
ls()

# 根據(jù)需要修改DEG的值
data <- DEG_limma_voom
colnames(data)


# 繪制火山圖
library(ggplot2)
colnames(data)
p <- ggplot(data=data, aes(x=logFC, y=-log10(adj.P.Val),color=regulated)) + 
     geom_point(alpha=0.5, size=1.8) + theme_set(theme_set(theme_bw(base_size=20))) + 
     xlab("log2FC") + ylab("-log10(FDR)") +scale_colour_manual(values = c('blue','black','red'))
p

image.png