In [ ]:
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
if (!require("airway", quietly = TRUE)) #使用数据集
  BiocManager::install("airway")
if (!require("org.Hs.eg.db", quietly = TRUE))
  BiocManager::install("org.Hs.eg.db")
if (!require("clusterProfiler", quietly = TRUE)) 
  BiocManager::install("clusterProfiler")
if (!require("DESeq2", quietly = TRUE)) #差异分析
  BiocManager::install("DESeq2")
In [ ]:
library(airway)
data(airway)#数据
head(airway)
expr <- assay(airway)
expr[1:5,1:5]

0eabd3882f181fd237239d84aa161fb.png

In [ ]:
dim(expr)

677ea5eb1b3405a420e6fc2d9d76795.png

In [ ]:
index1=c(rowSums(expr)>0)
class(index1)
exp2=expr[index1,]
dim(exp2)
#此处操作的目的在于将没有read比对上的基因行删去，简洁数据~

9f446ae8c25dd3ee0e66afc432e9d5d.png

In [ ]:
#获取分组信息
group_list <- colData(airway)$dex
group_list

ab7fa40bbad5b9aeaa5ded55a4e0b86.png

In [ ]:
library(DESeq2)
options(stringsAsFactors = F)
#构建dds对象,从表达矩阵countData和样品信息colData构建DESeqDataSet对象
colData <- data.frame(row.names=colnames(expr), group_list=group_list)
colData
dds <- DESeqDataSetFromMatrix(countData = round(exp2,0),colData = colData,design = ~ group_list)
dds
#数据过滤
dds <- dds[rowSums(counts(dds))>1,]#过滤低丰度数据
nrow(dds)

d020bdfa48197603da9dd77a0663c85.png

样本聚类

In [ ]:
#用vst来标化数据
vsd <- vst(dds, blind = F)
#再用dist这个函数计算样本间的距离
sampleDists <- dist(t(assay(vsd)))
library(ggplot2)
library(factoextra)
res <- hcut(sampleDists, k = 2, stand = TRUE)
fviz_dend(res, 
          # 加边框
          rect = TRUE,
          # 边框颜色
          rect_border="cluster",
          # 边框线条类型
          rect_lty=2,
          # 边框线条粗细
          lwd=1.2,
          # 边框填充
          rect_fill = T,
          # 字体大小
          cex = 1,
          # 字体颜色
          color_labels_by_k=T,
          # 平行放置
          horiz=T)

32b3e89939c3b18064ae7800c9190ad.png

DESeq2的标准化方法

（1）计算归一化系数sizeFactor

In [ ]:
dds <- estimateSizeFactors(dds)

（2）标准化之后的数据

In [ ]:
normalized_counts <- counts(dds,normalized=T)

差异表达分析

第一步：

In [ ]:
dds <- DESeq(dds)

2827809824cbb8a278f35a5b88721e4.png

In [ ]:
#提取差异分析结果
res <- results(dds)
resultsNames(dds) 
res=results(dds,contrast=c("group_list","trt","untrt"))
plotMA(res,ylim=c(-2,2))

11e7c931c1837235b2c70c48801a921.png

In [ ]:
names（res)

c4d9ce1f0f934ca667586133d0fbdf6.png

用summary这个函数，可以看到差异分析的结果，高表达和低表达的比例。低丰度基因所占的比例。

In [ ]:
plotMA(res,ylim=c(-2,2))

585cf740f421802b9ccef7d750be7ba.png

In [ ]:
DEG1=as.data.frame(res)
write.csv(x=DEG1,file="res.csv")

MA-plot上的点是每一个基因。横坐标是标准化后的counts平均值，纵坐标是log后的变化值。蓝色的点是p.adjust的值小于0.1的基因，由counts函数中的alpha 参数指定。 我们发现在左侧，有很多counts很小的基因，发生了很大的变化，但是没有明显意义。类似于从(1,3,9)变成了（20,12,3）他们的counts很小，波动性很大，对logFC产生了很大的影响。 GSEA分析中，排序就是按照logFC来进行的。按照这个结果往下做，GSEA那里，富集不到任何条目。

那就需要矫正。用的函数是lfcShrink

LFC校正

这一步，对于依赖logFC变化值的分析，很重要，比如GSEA分析。我们画一个MAplot图，

In [ ]:
resultsNames(dds)
dd1 <- lfcShrink(dds, coef="group_list_untrt_vs_trt", type="apeglm")
In [ ]:
plotMA(dd1,ylim=c(-2,2))

f27378b30c2ca7c690a679ba7e5ef38.png

这样，原先那些波动很大的基因，就被矫正了。而此时有LogFC以蓝色的点为主。

In [ ]:
names(dd1)

3f4bceb32b33c21e115b554d9ebb226.png

与res相比少了"stat"一列。这一步只会对LFC的值产生影响，p值是没有改变的，不会改变显著差异的基因总数。

排序和筛选

In [ ]:
res<-dd1
resOrdered=res[order(res$pvalue),] #order函数是给出从小到大排序后的位置(默认升序)
sum(res$padj<0.1, na.rm=TRUE) # na.rm 为remove na值，否则会影响统计量
sum(res$pvalue<0.05, na.rm=TRUE)  
sum(!is.na(res$pvalue))  
resSig=subset(resOrdered, padj < 0.1 & abs(log2FoldChange) > 1)
dim(resSig)
save(resSig,file="resSig.Rdata")
load(file="resSig.Rdata")
DEG2=as.data.frame(resSig)
write.csv(x=resSig, file = "results.csv")

5d1f569983d935f4002077068e44c19.png

热图

根据基因表达矩阵比对上的read数目绘制

这里选取最有意义的50个基因比对read数目为例

In [ ]:
library(pheatmap)
choose_gene=head(rownames(DEG2),50)    
#根据DEG2结果，提取目标基因名
choose_matrix=exp2[choose_gene,]
choose_matrix=t(scale(t(choose_matrix)))

一般绘制热图，都需要对数据进行归一化处理。这里用到了t()转置函数的原因是因为基因表达矩阵表中观测为基因，变量为样本。而scale()归一化函数是针对每一列单独处理的。因为每一个基因的本身情况各不相同，而且我们的目的是基于实验组与对照组的比较来研究基因表达变化，所以需要转置下表达矩阵（观测为样本，变量为基因），然后再进行对每一列的归一化处理；最后再进行一次转置，即得到下图结果。

9323c94668e3ef2c0f45799edeb1393.png

In [ ]:
pheatmap(choose_matrix,filename = 'DEG_top50_heatmap.png')

DEG_top501_heatmap.png

火山图

本质上就是根据DEG数据里基因的log2FoldChange与pvalue绘制点图

把点分成三类：上调基因、下调基因与无显著改变基因。

In [ ]:
DEG1=na.omit(DEG1)

（1）设置分类判断阈值并增加分类列

In [ ]:
logFC_cutoff=with(DEG1,mean(abs(log2FoldChange))+2*sd(abs(log2FoldChange)))
logFC_cutoff
2^logFC_cutoff   
DEG1$change=as.factor(ifelse(DEG1$pvalue<0.05 &  abs(DEG1$log2FoldChange)>logFC_cutoff,
                             ifelse(DEG1$log2FoldChange>logFC_cutoff,'UP','DOWN'),'NOT'))

（2）编写图中的注释信息

In [ ]:
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绘制火山图

In [ ]:
library(ggplot2)
g=ggplot(data=DEG1,
         aes(x=log2FoldChange,y=-log10(pvalue),   #这里将pvalue取负对数
             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'))   #设定颜色
ggsave(g,filename='volcanoo.png')

1b65d7a987f5d1d2b6f49bddaa758c8.png

把差异分析的结果转化成data.frame的格式

In [ ]:
#再把差异分析的结果转化成data.frame的格式
library(dplyr)
library(tibble)
res <- DEG2 %>% 
  data.frame() %>% 
  rownames_to_column("gene_id")
res

246f21251b712d98864a7f8a77cba1f.png

基因ID转换

In [ ]:
#基因ID转换
library(AnnotationDbi)
library(org.Hs.eg.db)
res$symbol <- mapIds(org.Hs.eg.db,
                     keys=res$gene_id,
                     column="SYMBOL",
                     keytype="ENSEMBL",
                     multiVals="first")
res$entrez <- mapIds(org.Hs.eg.db,
                     keys=res$gene_id,
                     column="ENTREZID",
                     keytype="ENSEMBL",
                     multiVals="first")
head(res)

ef4b5a1ce7017430e84c2300a08c060.png

制作geneLis

这里GSEA分析是用clusterprofier包. 使用这个包做GSEA，要制作一个genelist，这个是一个向量，他的内容是排序后的logFC值，他的名称是ENTREZID，而这两个我们都是不缺的，在上一步得到的差异结果中。

In [ ]:
library(dplyr)
gene_df <- res %>% 
  dplyr::select(gene_id,log2FoldChange,symbol,entrez) %>% 
  ## 去掉NA
  filter(entrez!="NA") %>% 
  ## 去掉重复
  distinct(entrez,.keep_all = T)

制作genelist的三部曲

In [ ]:
## 1.获取基因logFC
geneList <- gene_df$log2FoldChange
## 2.命名
names(geneList) = gene_df$entrez
## 3.排序很重要
geneList = sort(geneList, decreasing = TRUE)
head(geneList)

56bd01b1f8ff0d51554ba317208d54e.png

GSEA分析

In [ ]:
#完成KEGG的GSEA分析
library(clusterProfiler)
gseaKEGG <- gseKEGG(geneList     = geneList,
                  
                  organism     = 'hsa',
                  
                  nPerm        = 1000,
                  
                  minGSSize    = 10,
                  
                  pvalueCutoff = 0.9,
                  
                  verbose      = FALSE)
head(gseaKEGG)[,1:6]
gseaKEGG$Description

df176d4700050a1ebec1dbcd00f4800.png

In [ ]:
library(ggplot2)
dotplot(gseaKEGG,showCategory=12,split=".sign")+facet_grid(~.sign)

1246fbc31de8103eea477f017b10405.png

In [ ]:
gseaKEGG_results <- gseaKEGG@result
gseaKEGG_results

a43e393e07145640c21078cfe50a307.png

这时候，我们看到有一些通路是被激活的，有一些通路是被抑制的。比如Pathways in cancer是被抑制的，我们可以选取单个通路来作图。 把富集的结果转换成data.frame,找到Pathways in cancer的通路ID是hsa05200.

In [ ]:
library(enrichplot)
pathway.id = "hsa05200"
gseaplot2(gseaKEGG, 
          color = "red",
          geneSetID = pathway.id,
          pvalue_table = T)

320f0ee4477bba2f3d4c9d748306583.png

比如Calcium signaling pathway不是被抑制的，我们可以选取单个通路来作图。 把富集的结果转换成data.frame,找到Calcium signaling pathway的通路ID是hsa04020.

In [ ]:
library(enrichplot)
pathway.id = "hsa04020"
gseaplot2(gseaKEGG, 
          color = "red",
          geneSetID = pathway.id,
          pvalue_table = T)

0f8a337d628c45634bab6d9fa40319f.png