使用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富集分析以及不同类型的图表生成。