寫在前面：

分析之前需要兩個文件
文件1：countMatrix.txt數(shù)據(jù)，其中存儲著實(shí)驗(yàn)組與對照組的基因表達(dá)的count值（行是基因，列是樣本名）同時注意如果數(shù)據(jù)里面有缺失值需要進(jìn)行補(bǔ)缺失值http://m.itdecent.cn/p/ed14687738f6
。示例數(shù)據(jù)如下：

image.png

image.png

【注意】

如果是無重復(fù)樣本數(shù)據(jù)請參考：(27條消息) 使用edgeR進(jìn)行無重復(fù)差異表達(dá)分析_xuzhougeng blog-CSDN博客

1. 安裝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")}

2. 載入R包

library(DESeq2) # 差異表達(dá)分析
library(ggplot2) # 繪圖
library(pheatmap) # 熱圖
library(clusterProfiler) # 功能分析
library(org.Mm.eg.db)

3. 設(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.差異表達(dá)分析

 # 構(gòu)建DESeqDataSet對象
dds <- DESeqDataSetFromMatrix(countData = countMat, colData = clinical,design = ~Type)
 # 差異表達(dá)分析
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))) + # 散點(diǎn)顏色
  labs(x = expression(paste(log[2], "(Fold Change)")),
       y = expression(paste(-log[10], "(FDR)"))) + # 坐標(biāo)軸標(biāo)簽
  scale_colour_manual(values = c("#5175A4","grey","#D94E48")) + # 散點(diǎn)顏色
  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), # 坐標(biāo)軸標(biāo)題及字體大小
        legend.text = element_text(size=15),
        legend.title = element_text(size=15)) # 圖例標(biāo)題及字體大小
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],] # 為了實(shí)驗(yàn)方便，差異結(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計(jì)算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
# 計(jì)算基因的有效長度
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進(jìn)行富集分析

(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ù)庫使用

(https://string-db.org)

image.png

9.2 結(jié)果展示

![image.png](https://upload-images.jianshu.io/upload_images/20015369-d0287c9817f40370.png?imageMogr2/auto-
節(jié)點(diǎn)間的邊表示蛋白的相互作用，不同顏色邊表示不同的相互作用類型

image.png

9.3 其他網(wǎng)絡(luò)相關(guān)分析

1. 點(diǎn)擊analysis后，對于蛋白互作網(wǎng)絡(luò)的基因，進(jìn)行GO和KEGG富集分析的結(jié)果，頁面底部可以進(jìn)行下載

   
   
   image.png

2. 在Exports頁面，可以導(dǎo)出相互作用網(wǎng)絡(luò)的圖片，支持PNG，SVG格式，也可以導(dǎo)出對應(yīng)的相互作用表格和蛋白序列，注釋等信息

   
   
   image.png

3. 通過Cluster頁面來挖掘其中的子網(wǎng)sub network/module, 本質(zhì)上是對基因進(jìn)行聚類，屬于同一類的基因所構(gòu)成的相互作用網(wǎng)絡(luò)就是一個module，支持kmeans和MCL聚類(馬爾可夫聚類算法)，聚類的結(jié)果為TSV格式，從中可以看出哪些基因?qū)儆谕活?/p>
   

   image.png