一.样本表达分布

箱式图

##魔鬼操作一键清空，以防前面出现的变量名干扰后续的操作
rm(list = ls())
options(stringsAsFactors = F)

# 加载原始表达的数据
lname <- load(file = "Step01-airwayData.Rdata")
lname

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

# 样本表达总体分布-箱式图
library(ggplot2)
# 构造绘图数据
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()

为何即使组与组之间即使有生物学上的差异，但是箱式图中却没有表现出来呢?这是 基于前提假设：大多数情况下发生差异表达的基因只是很少的一部分，并且发生差异表达的基因不足以改变整体所有基因表达水平的分布；如果看样本的总体表达分布差异非常大的话，不是由于生物水平差异导致的，是误差导致的【比如不同人做实验、机子、protocol】

小提琴图

image.png

# 样本表达总体分布-小提琴图
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()  #打开画板要养成关掉的习惯

与箱式图结合发现，其实位于中位数的数量不是最多的

概率密度分布图

# 样本表达总体分布-概率密度分布图
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()

不仅可以看单个样本的表达趋势，还可以看多个样本的表达趋势；可以判断是否有异常的样本；图上不同样本的表达性就非常一致

以上三个图结合起来可以看总体样本表达分布图

二.样本的相关性

• 层次聚类树————通过计算点与点之间的空间距离对样本进行类别划分

主要看纵轴Height，也叫刻画样本之间的相似性；先从下到上看，每个样本先是单独的一类，距离相近的一类划分在一起

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

# 加载数据并检查
lname <- load(file = '../Analysis/data/Step01-airwayData.Rdata')
lname

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


## 1.样本之间的相关性-层次聚类树----
sampleTree <- hclust(dist(t(dat)), method = "average")  ##dist是计算距离【按照行来算的】，所以要用t来转置 ，把距离矩阵算完之后，开始聚类了使用的是：hclust函数
plot(sampleTree)

• PCA主成分分析————提取样本的综合特征，即主成分（第一主成分、第二主成分）对样本进行分类

image.png

## 2.样本之间的相关性-PCA----
# 第一步，数据预处理
dat <- as.data.frame(t(dat))  ##转置
dat$group_list <- group_list  ##添加group信息

# 第二步，绘制PCA图
library(FactoMineR)
library(factoextra)

# 画图仅需要数值型数据，去掉最后一列的分组信息
dat_pca <- PCA(dat[,-ncol(dat)], graph = FALSE)
class(dat_pca)

p <- fviz_pca_ind(dat_pca,
                  geom.ind = "text", # 只显示点，不显示文字
                  col.ind = dat$group_list, # 用不同颜色表示分组
                  palette = c("#00AFBB", "#E7B800"),
                  addEllipses = T, # 是否圈起来
                  legend.title = "Groups")
p

• 相关性分析

image.png

# 使用500个基因的表达量来做相关性图
library(corrplot)
dim(exprSet)

# 计算相关性
M <- cor(exprSet)  #M是一个对称矩阵
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")

# 绘制样本相关性的热图
anno <- data.frame(sampleType=group_list)  ##最重要的地方就是矩阵的行名是数据框的列名
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【原始输入数据都是count值】————用于高通量测序，不能用于基因芯片

基本理论

image.png
很多指标是对FC取了log，那么logFC>0，B相对于A就是上调；反之亦然【对表达水平较低的基因敏感，对表达水平高的不敏感】 image.png
可以理解为对P值做的矫正
image.png

limma差异分析

# 清空当前对象
rm(list = ls())
options(stringsAsFactors = F)

# 读取基因表达矩阵
lname <- load(file = "Step01-airwayData.Rdata")
lname

exprSet <- filter_count
# 检查表达谱
dim(exprSet)
exprSet[1:6,1:6]
table(group_list)

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

## 第一步，创建设计矩阵和对比：假设数据符合正态分布，构建线性模型
# 0代表x线性模型的截距为0
design <- model.matrix(~0+factor(group_list))  ##0是线性模型的截距
colnames(design) <- levels(factor(group_list))
rownames(design) <- colnames(exprSet)
design

两列八行

# 设置需要进行对比的分组，需要修改
comp <- 'trt-untrt'   ##不要设置反了
cont.matrix <- makeContrasts(contrasts=c(comp),levels = d esign)



## 第二步，进行差异表达分析
# 将表达矩阵转换为edgeR的DGEList对象
dge <- DGEList(counts=exprSet)

# 进行标准化
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)


## 第三步，提取过滤差异分析结果
tmp <- topTable(fit2, coef=comp, n=Inf,adjust.method="BH")
DEG_limma_voom <- na.omit(tmp)
head(DEG_limma_voom)

# 筛选上下调，设定阈值
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)

# 添加一列gene 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)

## 取表达差异倍数和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")


## 检查是否上下调设置错了
# 挑选一个差异表达基因
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)

# 读取基因表达矩阵信息
lname <- load(file = "../Analysis/data/Step01-airwayData.Rdata")
lname 

# 查看分组信息和表达矩阵数据
exprSet <- filter_count
dim(exprSet)
exprSet[1:6,1:6]
table(group_list)

# 加载包
library(DESeq2)

# 第一步，构建DESeq2的DESeq对象
colData <- data.frame(row.names=colnames(exprSet),group_list=group_list)
dds <- DESeqDataSetFromMatrix(countData = exprSet,colData = colData,design = ~ group_list)

# 第二步，进行差异表达分析
dds2 <- DESeq(dds)

# 提取差异分析结果，trt组对untrt组的差异分析结果
tmp <- results(dds2,contrast=c("group_list","trt","untrt"))
DEG_DESeq2 <- as.data.frame(tmp[order(tmp$padj),])
head(DEG_DESeq2)

# 去除差异分析结果中包含NA值的行
DEG_DESeq2 = na.omit(DEG_DESeq2)

# 筛选上下调，设定阈值
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)

# 添加一列gene 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)

## 取表达差异倍数和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")


## 检查是否上下调设置错了
# 挑选一个差异表达基因
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

差异结果可视化

image.png
image.png

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

# 加载原始表达矩阵
load(file = "Step01-airwayData.Rdata")

# 读取3个软件的差异分析结果
load(file = "Step03-limma_voom_nrDEG.Rdata")
load(file = "Step03-DESeq2_nrDEG.Rdata")
load(file = "Step03-edgeR_nrDEG.Rdata")
ls()

# 根据需要修改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