7. 差異分析

• 將基因計(jì)數(shù)導(dǎo)入 R/RStudio

工作流程完成后拯钻，您現(xiàn)在可以使用基因計(jì)數(shù)表作為 DESeq2 的輸入，使用 R 語(yǔ)言進(jìn)行統(tǒng)計(jì)分析。

7.1. 安裝R包

source("https://bioconductor.org/biocLite.R")
biocLite("DESeq2") ; library(DESeq2)
biocLite("ggplot2") ; library(ggplot2)
biocLite("clusterProfiler") ; library(clusterProfiler)
biocLite("biomaRt") ; library(biomaRt)
biocLite("ReactomePA") ; library(ReactomePA)
biocLite("DOSE") ; library(DOSE)
biocLite("KEGG.db") ; library(KEGG.db)
biocLite("org.Mm.eg.db") ; library(org.Mm.eg.db)
biocLite("org.Hs.eg.db") ; library(org.Hs.eg.db)
biocLite("pheatmap") ; library(pheatmap)
biocLite("genefilter") ; library(genefilter)
biocLite("RColorBrewer") ; library(RColorBrewer)
biocLite("GO.db") ; library(GO.db)
biocLite("topGO") ; library(topGO)
biocLite("dplyr") ; library(dplyr)
biocLite("gage") ; library(gage)
biocLite("ggsci") ; library(ggsci)

7.2. 導(dǎo)入表達(dá)矩陣

• 開(kāi)始導(dǎo)入文件夾中的 featureCounts 表。本教程將使用 DESeq2 對(duì)樣本組之間進(jìn)行歸一化和執(zhí)行統(tǒng)計(jì)分析。

# 導(dǎo)入基因計(jì)數(shù)表
# 使行名成為基因標(biāo)識(shí)符
countdata <- read.table("example/final_counts.txt", header = TRUE, skip = 1, row.names = 1)

# 從列標(biāo)識(shí)符中刪除 .bam 和 '..'
colnames(countdata) <- gsub(".bam", "", colnames(countdata), fixed = T)
colnames(countdata) <- gsub(".bam", "", colnames(countdata), fixed = T)
colnames(countdata) <- gsub("..", "", colnames(countdata), fixed = T)

# 刪除長(zhǎng)度字符列
countdata <- countdata[ ,c(-1:-5)]

# 查看 ID
head(countdata)  # 如下圖

countdata

7.3. 導(dǎo)入metadata

• 導(dǎo)入元數(shù)據(jù)文本文件续镇。 SampleID 必須是第一列。

# 導(dǎo)入元數(shù)據(jù)文件
# 使行名稱與 countdata 中的 sampleID 相匹配
metadata <- read.delim("example/metadata.txt", row.names = 1)

# 將 sampleID 添加到映射文件
metadata$sampleid <- row.names(metadata)

# 重新排序 sampleID 以匹配 featureCounts 列順序销部。
metadata <- metadata[match(colnames(countdata), metadata$sampleid), ]

# 查看 ID
head(metadata)  # 如下圖

metadata

7.4. DESeq2對(duì)象

• 根據(jù)計(jì)數(shù)和元數(shù)據(jù)創(chuàng)建 DESeq2 對(duì)象

# - countData : 基于表達(dá)矩陣
# - colData : 見(jiàn)上圖
# - design : 比較
ddsMat <- DESeqDataSetFromMatrix(countData = countdata,
                                 colData = metadata,
                                 design = ~Group)


# 查找差異表達(dá)基因
ddsMat <- DESeq(ddsMat)

7.5. 統(tǒng)計(jì)

• 獲取基因數(shù)量的基本統(tǒng)計(jì)數(shù)據(jù)

# 使用 FDR 調(diào)整 p-values 從檢測(cè)中獲取結(jié)果
results <- results(ddsMat, pAdjustMethod = "fdr", alpha = 0.05)

# 結(jié)果查看
summary(results)  # 如下圖

results

# 檢查 log2 fold change
## Log2 fold change is set as (LoGlu / HiGlu)
## Postive fold changes = Increased in LoGlu
## Negative fold changes = Decreased in LoGlu
mcols(results, use.names = T)  # 結(jié)果如下

mcols_result

8. 注釋基因symbol

經(jīng)過(guò)比對(duì)和總結(jié)摸航，我們只有帶注釋的基因符號(hào)。要獲得有關(guān)基因的更多信息舅桩，我們可以使用帶注釋的數(shù)據(jù)庫(kù)將基因符號(hào)轉(zhuǎn)換為完整的基因名稱和 entrez ID 以進(jìn)行進(jìn)一步分析酱虎。

• 收集基因注釋信息

# 小鼠基因組數(shù)據(jù)庫(kù)
library(org.Mm.eg.db) 

# 添加基因全名
results$description <- mapIds(x = org.Mm.eg.db,
                              keys = row.names(results),
                              column = "GENENAME",
                              keytype = "SYMBOL",
                              multiVals = "first")

# 添加基因 symbol
results$symbol <- row.names(results)

# 添加 ENTREZ ID
results$entrez <- mapIds(x = org.Mm.eg.db,
                         keys = row.names(results),
                         column = "ENTREZID",
                         keytype = "SYMBOL",
                         multiVals = "first")

# 添加 ENSEMBL
results$ensembl <- mapIds(x = org.Mm.eg.db,
                          keys = row.names(results),
                          column = "ENSEMBL",
                          keytype = "SYMBOL",
                          multiVals = "first")

# 取 (q < 0.05) 的基因
results_sig <- subset(results, padj < 0.05)

# 查看結(jié)果
head(results_sig)  # 如下圖

• 將所有重要結(jié)果寫(xiě)入 .txt 文件

# 將歸一化基因計(jì)數(shù)寫(xiě)入 .txt 文件
write.table(x = as.data.frame(counts(ddsMat), normalized = T), 
            file = 'normalized_counts.txt', 
            sep = '\t', 
            quote = F,
            col.names = NA)

# 將標(biāo)準(zhǔn)化基因計(jì)數(shù)寫(xiě)入 .txt 文件
write.table(x = counts(ddsMat[row.names(results_sig)], normalized = T), 
            file = 'normalized_counts_significant.txt', 
            sep = '\t', 
            quote = F, 
            col.names = NA)

# 將帶注釋的結(jié)果表寫(xiě)入 .txt 文件
write.table(x = as.data.frame(results), 
            file = "results_gene_annotated.txt", 
            sep = '\t', 
            quote = F,
            col.names = NA)

# 將重要的注釋結(jié)果表寫(xiě)入 .txt 文件
write.table(x = as.data.frame(results_sig), 
            file = "results_gene_annotated_significant.txt", 
            sep = '\t', 
            quote = F,
            col.names = NA)

9. 繪圖

有多種方法可以繪制基因表達(dá)數(shù)據(jù)。下面只列出了一些流行的方法擂涛。

9.1. PCA

# 將所有樣本轉(zhuǎn)換為 rlog
ddsMat_rlog <- rlog(ddsMat, blind = FALSE)

# 按列變量繪制 PCA
plotPCA(ddsMat_rlog, intgroup = "Group", ntop = 500) +
  theme_bw() +
  geom_point(size = 5) + 
  scale_y_continuous(limits = c(-5, 5)) +
  ggtitle(label = "Principal Component Analysis (PCA)", 
          subtitle = "Top 500 most variable genes") 

plotPCA

9.2. Heatmap

# 將所有樣本轉(zhuǎn)換為 rlog
ddsMat_rlog <- rlog(ddsMat, blind = FALSE)

# 收集30個(gè)顯著基因读串，制作矩陣
mat <- assay(ddsMat_rlog[row.names(results_sig)])[1:40, ]

# 選擇您要用來(lái)注釋列的列變量。
annotation_col = data.frame(
  Group = factor(colData(ddsMat_rlog)$Group), 
  Replicate = factor(colData(ddsMat_rlog)$Replicate),
  row.names = colData(ddsMat_rlog)$sampleid
)

# 指定要用來(lái)注釋列的顏色撒妈。
ann_colors = list(
  Group = c(LoGlu = "lightblue", HiGlu = "darkorange"),
  Replicate = c(Rep1 = "darkred", Rep2 = "forestgreen")
)

# 使用 pheatmap 功能制作熱圖恢暖。
pheatmap(mat = mat, 
         color = colorRampPalette(brewer.pal(9, "YlOrBr"))(255), 
         scale = "row",
         annotation_col = annotation_col, 
         annotation_colors = ann_colors,
         fontsize = 6.5, 
         cellwidth = 55,
         show_colnames = F)

pheatmap

9.3. Volcano

# 從 DESeq2 結(jié)果中收集倍數(shù)變化和 FDR 校正的 pvalue
## - 將 pvalues 更改為 -log10 (1.3 = 0.05)
data <- data.frame(gene = row.names(results),
                   pval = -log10(results$padj), 
                   lfc = results$log2FoldChange)

# 刪除任何以 NA 的行
data <- na.omit(data)

## If fold-change > 0 and pvalue > 1.3 (Increased significant)
## If fold-change < 0 and pvalue > 1.3 (Decreased significant)
data <- mutate(data, color = case_when(data$lfc > 0 & data$pval > 1.3 ~ "Increased",
                                       data$lfc < 0 & data$pval > 1.3 ~ "Decreased",
                                       data$pval < 1.3 ~ "nonsignificant"))

# 用 x-y 值制作一個(gè)基本的 ggplot2 對(duì)象
vol <- ggplot(data, aes(x = lfc, y = pval, color = color))

# 添加 ggplot2 圖層
vol +   
  ggtitle(label = "Volcano Plot", subtitle = "Colored by fold-change direction") +
  geom_point(size = 2.5, alpha = 0.8, na.rm = T) +
  scale_color_manual(name = "Directionality",
                     values = c(Increased = "#008B00", Decreased = "#CD4F39", nonsignificant = "darkgray")) +
  theme_bw(base_size = 14) + 
  theme(legend.position = "right") + 
  xlab(expression(log[2]("LoGlu" / "HiGlu"))) + 
  ylab(expression(-log[10]("adjusted p-value"))) + 
  geom_hline(yintercept = 1.3, colour = "darkgrey") + 
  scale_y_continuous(trans = "log1p") 

Volcano

9.4. MA

plotMA(results, ylim = c(-5, 5))

MA

9.5. Dispersions

plotDispEsts(ddsMat)

plotDispEsts

9.6. 單基因圖

# 將所有樣本轉(zhuǎn)換為 rlog
ddsMat_rlog <- rlog(ddsMat, blind = FALSE)

# 獲得最高表達(dá)的基因
top_gene <- rownames(results)[which.min(results$log2FoldChange)]

# 畫(huà)單基因圖
plotCounts(dds = ddsMat, 
           gene = top_gene, 
           intgroup = "Group", 
           normalized = T, 
           transform = T)

單基因圖

10. 通路富集

• 從差異表達(dá)基因中尋找通路

通路富集分析是基于單個(gè)基因變化生成結(jié)論的好方法。有時(shí)個(gè)體基因的變化是難以解釋狰右。但是通過(guò)分析基因的通路杰捂，我們可以收集基因反應(yīng)的視圖。

設(shè)置矩陣以考慮每個(gè)基因的 EntrezID 和倍數(shù)變化

# 刪除沒(méi)有任何 entrez 標(biāo)識(shí)符的基因
results_sig_entrez <- subset(results_sig, is.na(entrez) == FALSE)

# 創(chuàng)建一個(gè)log2倍數(shù)變化的基因矩陣
gene_matrix <- results_sig_entrez$log2FoldChange

# 添加 entrezID 作為每個(gè) logFC 條目的名稱
names(gene_matrix) <- results_sig_entrez$entrez

# 查看基因矩陣的格式
##- Names = ENTREZ ID
##- Values = Log2 Fold changes
head(gene_matrix)  # 如下圖

gene_matrix

10.1. KEGG

• 使用 KEGG 數(shù)據(jù)庫(kù)豐富基因

kegg_enrich <- enrichKEGG(gene = names(gene_matrix),
                          organism = 'mouse',
                          pvalueCutoff = 0.05, 
                          qvalueCutoff = 0.10)

# 結(jié)果可視化
barplot(kegg_enrich, 
        drop = TRUE, 
        showCategory = 10, 
        title = "KEGG Enrichment Pathways",
        font.size = 8)

KEGG

10.2. GO

• 使用 Gene Onotology 豐富基因

go_enrich <- enrichGO(gene = names(gene_matrix),
                      OrgDb = 'org.Mm.eg.db', 
                      readable = T,
                      ont = "BP",
                      pvalueCutoff = 0.05, 
                      qvalueCutoff = 0.10)

# 結(jié)果可視化
barplot(go_enrich, 
        drop = TRUE, 
        showCategory = 10, 
        title = "GO Biological Pathways",
        font.size = 8)

GO

11. 通路可視化

Pathview 是一個(gè)包棋蚌，它可以獲取顯著差異表達(dá)基因的 KEGG 標(biāo)識(shí)符嫁佳，還可以與 KEGG 數(shù)據(jù)庫(kù)中發(fā)現(xiàn)的其他生物一起使用，并且可以繪制特定生物的任何 KEGG 途徑谷暮。

# 加載包
biocLite("pathview") ; library(pathview)

# 可視化通路 (用 fold change) 
## pathway.id : KEGG pathway identifier
pathview(gene.data = gene_matrix, 
         pathway.id = "04070", 
         species = "mouse")

pathview

歡迎Star -> 學(xué)習(xí)目錄

國(guó)內(nèi)鏈接 -> 學(xué)習(xí)目錄