使用bioconductorlimma進行基因表達數據分析

一、安裝和載入limma包

Limma是一款R軟體的包，可用於在微陣列和RNA-Seq下處理基因表達數據。首先，我們需要安裝limma包。代碼如下：

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("limma")

載入limma包：

library(limma)

二、讀取基因表達數據

我們可以使用read.table函數來讀取數據。在下面的示例中，數據文件名為”Gene_expression.txt”，每一行代表一個基因。首列包含基因名字，接下來的列包含的是不同樣本的基因表達數據。這裡，我們假設第一列是基因名字，第二列到第六列是5個不同組織（樣本）的基因表達數據。示例代碼如下：

data=read.table("Gene_expression.txt",header=TRUE,row.names=1)

三、數據預處理

在進行微陣列或RNA-Seq數據分析之前，我們必須先進行數據預處理。這包括數據歸一化、過濾掉低表達基因、Batch效應的去除等。在這裡，我們使用normalizeBetweenArrays來進行數據標準化。相關代碼如下：

dataNorm=normalizeBetweenArrays(data,method="quantile")

在進行limma分析之前，我們需要在數據中去除所有低表達的基因。

keep=filterByExpr(dataNorm)
dataNormFilter=dataNorm[keep,]

接下來，我們進行Batch效應的去除。

batch = c(1,1,2,2,3)
design = model.matrix(~0+batch)
fit = lmFit(dataNormFilter, design)
contrMatrix = makeContrasts(batch2-batch1, levels=design)
fit2 = contrasts.fit(fit, contrasts=contrMatrix)
fit2 = eBayes(fit2)

四、配對差異分析

接下來，我們進行配對差異分析，以比較2個組織（樣本）之間的差異性。

例如，我們想要比較第一組織和第二組織的差異性：

treatment=c(1,1,2,2,1)
contrasts=makeContrasts(treatment1-treatment2,levels=design)
fit2 = contrasts.fit(fit,contrasts=contrasts)
fit2 = eBayes(fit2)

此處我們存儲了差異表達儲存

result = topTable(fit2,coef=1,number=1000)

五、GO/GSEA富集分析

使用富集分析可以幫助我們了解差異表達基因的生物學意義。

我們首先需要用DEGs獲得差異表達基因的集合。例如，我們使用默認參數在一個基因集合上執行Label的卡方富集分析，來計算差異表達基因集合的富集結果。相應代碼如下：

de <- rownames(result)  # 使用導出的差異表達基因
ego <- enrichGO(gene = de,
                OrgDb = org.Mm.eg.db,
                keyType = "ENSEMBL",
                ont = "BP",
                pvalueCutoff = 0.05,
                qvalueCutoff = 1,
                readable = TRUE)
ego <- setReadable(ego, 'org.Mm.eg.db', 'ENSEMBL')  # 用基因符號代替ENSEMBL ID
head(as.data.frame(ego), 10)

六、不同圖展示

最後，我們可以使用不同類型的圖表展示分析結果。例如，我們可以使用pheatmap包來生成熱圖。

library(pheatmap)
mat <- assay(dataNormFilter)
rownames(mat)=paste("G",1:nrow(mat),sep="")
samples <- data.frame(Tissue=colnames(mat), row.names=colnames(mat))
colnames(samples) <- c("Tissue")
colors <- colorRampPalette(rev(brewer.pal(9, "RdBu")))(75)
pheatmap(mat,
         color=colors,
         cluster_rows=TRUE,
         cluster_cols=TRUE,
         annotation_col=samples,
         fontsize=8,
         main="Heatmap of Gene Expression")

總結

在這篇文章中，我們介紹了如何使用bioconductorlimma進行基因表達數據分析。我們展示了如何讀取數據、進行數據預處理、進行配對差異分析、進行GO/GSEA富集分析以及不同類型的圖表生成。