寫在前面：
分析之前需要兩個文件
文件1：countMatrix.txt數(shù)據(jù)蝶溶，其中存儲著實驗組與對照組的基因表達的count值（行是基因砾隅，列是樣本名）同時注意如果數(shù)據(jù)里面有缺失值需要進行補缺失值http://www.reibang.com/p/ed14687738f6
。示例數(shù)據(jù)如下：
image.png
文件2：samplelInfo.txt數(shù)據(jù)嘶是，其中存儲著對樣本的介紹钙勃。如下：
image.png
【注意】
如果是無重復(fù)樣本數(shù)據(jù)請參考：(27條消息) 使用edgeR進行無重復(fù)差異表達分析_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)
- 設(shè)置工作路徑; 讀取數(shù)據(jù)
##設(shè)置路徑方便文件讀取輸出
setwd("C:\\Users\\experiment-RNA_seq")#注意需要雙反斜線
# 讀入count矩陣
countMat <- read.table("countMatrix.txt", sep = "\t", header = T, row.names = 1, check.names = F)
#調(diào)整gene ID ENSG00000005513.9->ENSG00000005513 并設(shè)為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.差異表達分析
# 構(gòu)建DESeqDataSet對象
dds <- DESeqDataSetFromMatrix(countData = countMat, colData = clinical,design = ~Type)
# 差異表達分析
dds2 <- DESeq(dds)
# 提取分析結(jié)果
res <- results(dds2)
5.繪制火山圖
# 火山圖數(shù)據(jù)構(gòu)建
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 差異結(jié)果
DEG <- DEG[order(DEG$padj)[1:300],] # 為了實驗方便，差異結(jié)果排序后取前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轉(zhuǎn)換
DEG[transID$ENSEMBL, 7] <- transID$SYMBOL
colnames(DEG)[7] <- "symbol"
DEG %>% head
7.差異基因聚類分析
#從count計算TPM
BiocManager::install("biomaRt")
library(biomaRt)
#查看基因組參數(shù)
mart = useMart('ensembl')
listDatasets(mart)
#你需要哪個基因組聂喇，就復(fù)制它在dataset列里的詞辖源，放在下面這行的`dataset = `參數(shù)里
#此處以人類為例蔚携，植物參考注一
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")
# 從輸入數(shù)據(jù)里提取基因名
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 選擇物種及注釋數(shù)據(jù)庫
image.png
8.4 KEGG注釋結(jié)果展示
image.png
9. 差異基因PPI網(wǎng)絡(luò)構(gòu)建
9.1 STRING數(shù)據(jù)庫使用
image.png
9.2 結(jié)果展示
點間的邊表示蛋白的相互作用，不同顏色邊表示不同的相互作用類型
image.png
9.3 其他網(wǎng)絡(luò)相關(guān)分析
-
點擊analysis后克饶，對于蛋白互作網(wǎng)絡(luò)的基因酝蜒，進行GO和KEGG富集分析的結(jié)果，頁面底部可以進行下載
image.png -
在Exports頁面矾湃，可以導(dǎo)出相互作用網(wǎng)絡(luò)的圖片亡脑，支持PNG，SVG格式邀跃，也可以導(dǎo)出對應(yīng)的相互作用表格和蛋白序列霉咨，注釋等信息
image.png -
通過Cluster頁面來挖掘其中的子網(wǎng)sub network/module, 本質(zhì)上是對基因進行聚類，屬于同一類的基因所構(gòu)成的相互作用網(wǎng)絡(luò)就是一個module拍屑，支持kmeans和MCL聚類(馬爾可夫聚類算法)躯护，聚類的結(jié)果為TSV格式，從中可以看出哪些基因?qū)儆谕活?/p>
image.png
