写在前面：
分析之前需要两个文件
文件1：countMatrix.txt数据，其中存储着实验组与对照组的基因表达的count值（行是基因，列是样本名）同时注意如果数据里面有缺失值需要进行补缺失值https://www.jianshu.com/p/ed14687738f6
。示例数据如下：
image.png
文件2：samplelInfo.txt数据，其中存储着对样本的介绍。如下：
image.png
【注意】
如果是无重复样本数据请参考：(27条消息) 使用edgeR进行无重复差异表达分析_xuzhougeng blog-CSDN博客
- 安装R包
if(!requireNamespace("BiocManager", quietly = TRUE)){
install.packages("BiocManager")}
if(!requireNamespace("DESeq2", quietly = TRUE)){
BiocManager::install("DESeq2")}
if(!requireNamespace("clusterProfiler", quietly = TRUE)){
BiocManager::install("clusterProfiler")}
if(!requireNamespace("org.Mm.eg.db", quietly = TRUE)){
BiocManager::install("org.Mm.eg.db")}
if(!requireNamespace("ggplot2", quietly = TRUE)){
install.packages("ggplot2")}
if(!requireNamespace("pheatmap", quietly = TRUE)){
install.packages("pheatmap")}
- 载入R包
library(DESeq2) # 差异表达分析
library(ggplot2) # 绘图
library(pheatmap) # 热图
library(clusterProfiler) # 功能分析
library(org.Mm.eg.db)
- 设置工作路径; 读取数据
##设置路径方便文件读取输出
setwd("C:\\Users\\experiment-RNA_seq")#注意需要双反斜线
# 读入count矩阵
countMat <- read.table("countMatrix.txt", sep = "\t", header = T, row.names = 1, check.names = F)
#调整gene ID ENSG00000005513.9->ENSG00000005513 并设为countMat 行名
rownames(countMat) <- rownames(countMat) %>% strsplit(., split = "[.]") %>% lapply(., function(x) x[1]) %>% unlist
samplelinfo <- read.table("samplelInfo.txt", sep = "\t", header = T, row.names = 1)
4.差异表达分析
# 构建DESeqDataSet对象
dds <- DESeqDataSetFromMatrix(countData = countMat, colData = clinical,design = ~Type)
# 差异表达分析
dds2 <- DESeq(dds)
# 提取分析结果
res <- results(dds2)
5.绘制火山图
# 火山图数据构建
volcano_data <- data.frame(log2FC = res$log2FoldChange,
padj = res$padj,
sig = "not",
stringsAsFactors = F) # 提取差异log2FoldChange和padj
volcano_data$sig[volcano_data$padj < 0.01 & volcano_data$log2FC > 1] <- "up"
volcano_data$sig[volcano_data$padj < 0.01 & volcano_data$log2FC < -1] <- "down"
# 火山图绘制
theme_set(theme_bw())
p <- ggplot(volcano_data, aes(log2FC, -1*log10(padj))) +
geom_point(aes(colour = factor(sig))) + # 散点颜色
labs(x = expression(paste(log[2], "(Fold Change)")),
y = expression(paste(-log[10], "(FDR)"))) + # 坐标轴标签
scale_colour_manual(values = c("#5175A4","grey","#D94E48")) + # 散点颜色
geom_hline(yintercept = -log10(0.01), linetype = 4) +
geom_vline(xintercept = c(-1,1), linetype = 4) + # 画线
theme(panel.grid = element_blank()) + # 去背景
theme(axis.text = element_text(size=15),
axis.title = element_text(size=15), # 坐标轴标题及字体大小
legend.text = element_text(size=15),
legend.title = element_text(size=15)) # 图例标题及字体大小
p
image.png
6.提取差异基因
DEG <- subset(res, padj < 0.01 & abs(log2FoldChange) > 1) # 提取 padj<0.01 和 log2FoldChange 差异结果
DEG <- DEG[order(DEG$padj)[1:300],] # 为了实验方便，差异结果排序后取前top300
paste("The number of up-regulated genes:", sum(DEG$log2FoldChange>1)) %>% print()
transID <- bitr(rownames(DEG), fromType="ENSEMBL",
toType="SYMBOL", OrgDb="org.Mm.eg.db") # gene ID转换
DEG[transID$ENSEMBL, 7] <- transID$SYMBOL
colnames(DEG)[7] <- "symbol"
DEG %>% head
7.差异基因聚类分析
#从count计算TPM
BiocManager::install("biomaRt")
library(biomaRt)
#查看基因组参数
mart = useMart('ensembl')
listDatasets(mart)
#你需要哪个基因组，就复制它在dataset列里的词，放在下面这行的`dataset = `参数里
#此处以人类为例，植物参考注一
bmart <- biomaRt::useMart(biomart = "ENSEMBL_MART_ENSEMBL",
dataset = "mmusculus_gene_ensembl",
host = "www.ensembl.org")
bmart <- biomaRt::useMart(biomart = "ENSEMBL_MART_ENSEMBL",
dataset = "hsapiens_gene_ensembl",
host = "www.ensembl.org")
# 从输入数据里提取基因名
feature_ids <- rownames(countMat)##就是你需要ENSEMBL的例如ENSG00000000003
attributes = c(
"ensembl_gene_id",
#"hgnc_symbol",
"chromosome_name",
"start_position",
"end_position"
)
filters = "ensembl_gene_id"
feature_info <- biomaRt::getBM(attributes = attributes,
filters = filters,
values = feature_ids, mart = bmart)
mm <- match(feature_ids, feature_info[[filters]])
feature_info_full <- feature_info[mm, ]
rownames(feature_info_full) <- feature_ids
# 计算基因的有效长度
eff_length <- abs(feature_info_full$end_position - feature_info_full$start_position)
feature_info_full <- cbind(feature_info_full,eff_length)
eff_length2 <- feature_info_full[,c(1,5)]
x <- countMat/ eff_length2$eff_length
tpm = t(t(x)/sum(x[,1])*1000000)
head(tpm)
#差异基因聚类分析
TPMMat_DEG <- tpm[rownames(DEG),]
ann_colors <- list(Type = c(Normal ="#5175A4",Tumor ="#D94E48")) # 列注释颜色
pheatmap(TPMMat_DEG,
cluster_rows = T, cluster_cols = T, # 行列均聚类
show_rownames=F, show_colnames = F, # 行名列名显示
clustering_method = "ward.D2", # 聚类方法
annotation_col = clinical, annotation_colors = ann_colors, # 列注释
scale="row", breaks = c(seq(-2,2, length=101)), border_color = NA)
image.png
8.差异基因功能富集分析
8.1 使用DAVID进行富集分析
(https://david.ncifcrf.gov)(https://david.ncifcrf.gov/content.jsp?file=functional_annotation.html)
image.png
8.2 上传差异基因
image.png
8.3 选择物种及注释数据库
image.png
8.4 KEGG注释结果展示
image.png
9. 差异基因PPI网络构建
9.1 STRING数据库使用
image.png
9.2 结果展示
image.png
9.3 其他网络相关分析
-
点击analysis后，对于蛋白互作网络的基因，进行GO和KEGG富集分析的结果，页面底部可以进行下载
image.png -
在Exports页面，可以导出相互作用网络的图片，支持PNG，SVG格式，也可以导出对应的相互作用表格和蛋白序列，注释等信息
image.png -
通过Cluster页面来挖掘其中的子网sub network/module, 本质上是对基因进行聚类，属于同一类的基因所构成的相互作用网络就是一个module，支持kmeans和MCL聚类(马尔可夫聚类算法)，聚类的结果为TSV格式，从中可以看出哪些基因属于同一类
image.png
