【轉錄組入門】7：差異基因分析

本文轉載自查看原文 2018-07-03 21:47 5049 生物信息學/ RNA-seq

作業要求：

使用R語言，載入表達矩陣，然后設置好分組信息，統一用DEseq2進行差異分析，當然也可以走走edgeR或者limma的voom流程。

基本任務是得到差異分析結果，進階任務是比較多個差異分析結果的異同點。

【1】安裝DESeq2

1 # 下面是在R語言中操作
2 # 載入安裝工具
3 > source("http://bioconductor.org/biocLite.R")
4 # 安裝包
5 > biocLite("DESeq2")
6 # 載入包
7 > library("DESeq2")

DESeq2對於輸入數據的要求：

1.DEseq2要求輸入數據是由整數組成的矩陣。

2.DESeq2要求矩陣是沒有標准化的。

【2】DESeq2進行差異表達分析

DESeq2分析差異表達基因簡單來說只有三步：構建dds矩陣，標准化，以及進行差異分析。

# dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design= ~ batch + condition) #~在R里面用於構建公式對象，~左邊為因變量，右邊為自變量。

# dds <- DESeq(dds) #標准化

# res <- results(dds, contrast=c("condition","treated","control")) #差異分析結果

【3】構建dds矩陣

1 > library(DESeq2)     # 加載包
2 > countData <- raw_count_filter[2:7]     # 中括號中的數量要與condition中數量一致
3 > condition <- factor(c("control","control","control","hypoxia","hypoxia","hypoxia"))
4 > colData <- data.frame(row.names=colnames(countData),condition)
5 # raw_count_filter：是所有樣品的count按照gene id融合后生成的矩陣。行為各個基因，列為各個樣品，中間為計算reads。

1 # 正式構建dds矩陣
2 > dds <- DESeqDataSetFromMatrix(countData,DataFrame(condition),design=~condition)
3 # 注意，condition前面是波浪線
4 > head(dds)     # 查看一下構建好的矩陣

【4】對原始dds進行標准化

1 > dds2 <- DESeq(dds)     # 對原始dds進行normalize

1 # 查看結果的名稱，本次實驗中是“intercept”,”condition_akap95_vs_control”
2 > resultsNames(dds2)
3 # 將結果用results()函數來獲取，賦值給res變量
4 > res <- results(dds2)
5 # summary一下，看一下結果的概要信息
6 > summary(res)

【5】提取差異分析結果

1 # 獲取padj（p值經過多重校驗校正后的值）小於0.05，log2FoldChange的絕對值大於1的差異表達基因。
2 > table(res$padj<0.05)     # 取p值小於0.05的結果
3 > res <- res[order(res$padj).]
4 > diff_gene_deseq2 <- subset(res,padj<0.05 & (log2FoldChange > 1 | log2FoldChange < -1))
5 > diff_gene_deseq2 <- row.names(diff_gene_deseq2)
6 > resdata <- merge(as.data.frame(res),as.data.frame(counts(dds2,normalize=TRUE)),by="row.names",sort=FALSE)
7 # 將差異表達分析結果輸出到csv文件
8 > write.csv(resdata,"d:/cs_file/r_file/diff_gene_control_vs_hypoxia.csv",row.names=FALSE)

補充：

用edgeR進行基因差異表達分析

 1 # 加載edgeR包
 2 library(edgeR)
 3 ##跟DESeq2一樣，導入數據，預處理（用了cpm函數）
 4 exprSet<- read.table(file = "mouse_all_count.txt", sep = "\t", header = TRUE, row.names = 1, stringsAsFactors = FALSE)
 5 group_list <- factor(c(rep("control",2),rep("Akap95",2)))
 6 exprSet <- exprSet[rowSums(cpm(exprSet) > 1) >= 2,]
 7                  
 8 ##設置分組信息，並做TMM標准化
 9 exprSet <- DGEList(counts = exprSet, group = group_list)
10 exprSet <- calcNormFactors(exprSet)
11                  
12 ##使用qCML（quantile-adjusted conditional maximum likelihood）估計離散度（只針對單因素實驗設計）
13 exprSet <- estimateCommonDisp(exprSet)
14 exprSet <- estimateTagwiseDisp(exprSet)
15                  
16 ##尋找差異gene(這里的exactTest函數還是基於qCML並且只針對單因素實驗設計)，然后按照閾值進行篩選即可
17 et <- exactTest(exprSet)
18 tTag <- topTags(et, n=nrow(exprSet))
19 diff_gene_edgeR <- subset(tTag$table, FDR < 0.05 & (logFC > 1 | logFC < -1))
20 diff_gene_edgeR <- row.names(diff_gene_edgeR)
21 write.csv(tTag$table,file = "control_vs_Akap95_edgeR.csv")