生信云平台再升级！Galaxy 新增工具 PyDESeq2，RNA-seq 差异分析一键搞定-平芜编程栈

差异基因分析是转录组下游分析的第一步，我们前面写了目前最流行的 3 个差异分析软件的介绍：

差异表达分析三剑客：DESeq2、edgeR、limma 到底该怎么选？

目前可以说，DESeq2 是差异分析的金标准，因此我之前专门写了一篇文章详细介绍了这个软件在平台的用法：

任何新手小白都可以使用这个差异基因分析神器：DESeq2 | 零代码，无需编程

不过，目前平台 R 语言版本的 DESeq2 工具的使用个人觉得还是略显复杂。它主要以每一个样本的原始计数文件作为输入，这有一个问题，就是当我们要分析的样本很多时，这就有些不方便了。

转录组上游分析完成后，通常会将所有样本的表达量数据合并到一个表格，得到表达量矩阵。那么，为什么不直接以这个表达量矩阵为输入，再指定要比较的分组情况，事情不更简单吗？

鉴于平台已经有 R 版本的 DESeq2 了，并且目前它也运行良好，似乎没必要对其进行升级。我们干脆新增一个 Python 版的 DESeq2 - PyDESeq2，这样平台 R 和 Python 版本的差异分析工具都有了。

进入网站

https://usegalaxy.cn

测试数据

测试数据位置：库 / Libraryies / RNA-seq / pasilla，找到两个文件，并将它们添加到历史记录中：

• pasilla_gene_counts.tsv，表达量矩阵
• pasilla_sample_annotation.csv，样本注释信息

找到工具

在平台搜索：pydeseq2，找到工具：

开始分析

需要输入两类信息：

• 表达量矩阵：行为基因，列为样本，Tab 键隔开。普通转录组表达量矩阵默认就是这个格式。
• 实验设计信息：可以手工输入，或者通过元数据表输入。

简单模式 - 手动输入样本名

这种模式针对样本量不多的情况，直接将样本名称粘贴到对应的文本框中就好了，名称之间可以用空格、制表符或换行进行分隔。

当样本比较多、设计比较复杂时，可以使用下面介绍的高级模式。

高级模式 - 使用元数据表

•元数据表 (Metadata File)：这种模式你先得制作一个表，将样本的相关信息填写在其中，记得保存为 csv 文件（英文逗号分隔）。

sample,condition,type treated1,treated,single treated2,treated,paired treated3,treated,paired untreated1,untreated,single untreated2,untreated,single untreated3,untreated,paired untreated4,untreated,paired

•设计公式 (Design Formula)：condition
•比较因子 (Contrast Factor)：condition
•测试组水平 (Test Level)：treated
•基准组水平 (Base Level) ：untreated

上述演示了单因素设计的情形，如果是多因素设计，可以这样书写设计公式：

注意你最关心的因子要写在最后，比如最关注 condition，就将它放到最后。

不管是单因素设计，还是多因素设计，都需要指明最终用于计算 foldChange 以及 pvalue 等值的分组的情况。简单来说，你希望计算哪一列中，哪两个组的比较情况。

比如要计算 condition 列中：

treated vs untreated

这样最终算出来的变化倍数（foldChange）是 treated 组的均值，除以 untreated 组的均值。

好了，是不是超简单？

结果解读

差异表达结果表 (CSV)包含详细统计结果和标准化表达量的表格。关键列说明如下：

•baseMean: 所有样本的标准化平均表达量。
•baseMean(Group): 某一组 (测试组或基准组) 的标准化平均表达量。
•FoldChange: 线性表达倍数变化 (计算公式为2 ^ log2FoldChange)。
•log2FoldChange: Log2 转换后的差异倍数。注意: 该值已经过LFC Shrinkage (收缩校正)处理，能够更准确地估算低表达量或高离散度基因的变化，更适合用于排序和绘图。
•pvalue: Wald Test 统计检验的 P 值。
•padj: 经 Benjamini-Hochberg 方法校正后的 P 值 (即 FDR)。
•[样本名]: 表格的最后部分包含每个样本经过 DESeq2 标准化后的具体表达量。

与 R 语言比较

单因素设计 - Galaxy 平台

gene_id baseMean log2FoldChange lfcSE pvalue padj FBgn0039155 730.6 -4.6 0.17 4.3e-163 0 FBgn0025111 1501.4 2.9 0.13 5.8e-113 0 FBgn0029167 3706.1 -2.2 0.1 2.2e-108 0 FBgn0003360 4343 -3.2 0.15 1.5e-106 0 FBgn0035085 638.2 -2.5 0.14 1.3e-76 0 FBgn0039827 261.9 -4.1 0.23 3.4e-72 0 FBgn0034736 225.9 -3.5 0.21 1.1e-60 0 FBgn0029896 489.9 -2.4 0.15 5.8e-58 0 FBgn0000071 342.2 2.7 0.18 4.1e-49 0

单因素设计 - R 语言

library("DESeq2") pasCts <- system.file("extdata", "pasilla_gene_counts.tsv.gz", package="DESeq2", mustWork=TRUE) pasAnno <- system.file("extdata", "pasilla_sample_annotation.csv", package="DESeq2", mustWork=TRUE) cts <- as.matrix(read.csv(pasCts,sep="\t",row.names="gene_id")) coldata <- read.csv(pasAnno, row.names=1) coldata <- coldata[,c("condition","type")] coldata$condition <- factor(coldata$condition) coldata$type <- factor(coldata$type) rownames(coldata) <- sub("fb", "", rownames(coldata)) all(rownames(coldata) %in% colnames(cts)) cts <- cts[, rownames(coldata)] all(rownames(coldata) == colnames(cts)) dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~ condition) smallestGroupSize <- 3 keep <- rowSums(counts(dds) >= 10) >= smallestGroupSize dds <- dds[keep,] # ------------------------------------------------------------------------- # 单因素设计 # ------------------------------------------------------------------------- dds$condition <- factor(dds$condition, levels = c("untreated","treated")) dds$condition dds$condition <- droplevels(dds$condition) dds <- DESeq(dds) res <- results(dds) resLFC <- lfcShrink(dds, coef="condition_treated_vs_untreated", type="apeglm") results = resLFC[order(resLFC$padj), ] results$baseMean = round(results$baseMean, 1) results$log2FoldChange = round(results$log2FoldChange, 1) results$lfcSE = round(results$lfcSE, 2) results$pvalue = formatC(results$pvalue, format = "e", digits = 1) results$padj = round(results$padj, 4) results write.csv(results, 'pasilla_deg_results.tsv')

可以看到，不说分毫不差吧，也几乎完全一样。

多因素设计 - Galaxy 平台

gene_id baseMean log2FoldChange lfcSE pvalue padj FBgn0003360 4343 -3.1 0.11 1.4e-175 0 FBgn0026562 43909.3 -2.5 0.09 1.6e-171 0 FBgn0039155 730.6 -4.6 0.17 3.6e-170 0 FBgn0025111 1501.4 2.8 0.1 4.1e-166 0 FBgn0029167 3706.1 -2.2 0.1 9.2e-114 0 FBgn0035085 638.2 -2.6 0.13 3.8e-91 0 FBgn0039827 261.9 -4.2 0.21 4.6e-86 0 FBgn0034736 225.9 -3.5 0.2 3.2e-69 0 FBgn0000071 342.2 2.6 0.15 8.6e-64 0

多因素设计 - R 语言

# ------------------------------------------------------------------------- # 多因素设计 # ------------------------------------------------------------------------- ddsMF <- dds ddsMF levels(ddsMF$type) levels(ddsMF$type) <- sub("-.*", "", levels(ddsMF$type)) levels(ddsMF$type) design(ddsMF) <- formula(~ type + condition) ddsMF <- DESeq(ddsMF) resMF <- results(ddsMF) results = resMF[order(resMF$padj), ] results$baseMean = round(results$baseMean, 1) results$log2FoldChange = round(results$log2FoldChange, 1) results$lfcSE = round(results$lfcSE, 2) results$pvalue = formatC(results$pvalue, format = "e", digits = 1) results$padj = round(results$padj, 4) results write.csv(results, 'pasilla_deg_results_multi-factor_designs.tsv')