前面在学习python的时候，在学员群里看到很多人问，能不能用python来做bulk转录组的差异表达分析 ，这就来看看！
我搜了一下相关的资料，发现有py版本的DESeq2模块：pydeseq2。网址：
https://pypi.org/project/pydeseq2/
PyDESeq2 是一个用 Python 编写的差异表达分析工具，它实现了最初在 R 语言中开发的 DESeq2 方法。这个工具主要是为了帮助 Python 用户进行 bulk RNA 测序数据的差异表达分析。尽管 PyDESeq2 是 DESeq2 的重新实现，但在结果值和功能上可能存在一些差异。
目前，PyDESeq2 提供的功能与 DESeq2 1.34.0 版本的默认设置相似，支持单因素和多因素分析，包括分类和连续因素，并使用 Wald 检验。
模块安装 我这里直接安装在conda小环境 sc 中：
## bash 命令 conda activate sc pip install pydeseq2 -i https://mirrors.tuna.tsinghua.edu.cn/pypi/web/simple 加载模块：
import pandas as pd import numpy as np import anndata as ad import math import seaborn as sns import matplotlib.pyplot as plt from pydeseq2.dds import DeseqDataSet from pydeseq2.ds import DeseqStats 数据 数据为：bulk转录组的count矩阵，以及一个样本表型文件group.txt，里面是两种样本类型：实验组DEX，对照组untreated
现在读取count值，需要对读入的数据进行转置，因为使用pydeseq2包进行分析时，count矩阵需要的是行为样本，列为基因名称 ，和R语言中的DESeq2包刚好相反。
## 读取count矩阵 df = pd.read_table( "counts.txt" ,index_col= 0 ).T # 查看前6列 df.head().iloc[:,range( 6 )] 样本分组信息：
注意：分组信息中不能有NA值，NaN，group的行名需要与count矩阵中样本的顺序一一对应
## 读取样本分组 group = pd.read_table( "group.txt" , index_col= 0 ) group 过滤低表达基因 筛选掉那些总读取计数少于10的基因：
## 数据过滤 genes_to_keep = df.columns[df.sum(axis= 0 ) >= 10 ] df = df[genes_to_keep] df.shape 过滤前后基因数分别为：21930，21600
差异分析 首先构建DeseqDataSet 对象：
# 构建DeseqDataSet 对象 inference = DefaultInference(n_cpus= 8 ) dds = DeseqDataSet( counts=df, metadata=group, design= "~sample_title" , refit_cooks= True , inference=inference, # n_cpus=8, # n_cpus can be specified here or in the inference object ) 一旦初始化了一个 DeseqDataSet 对象，我们就可以使用 deseq2() 方法来拟合离散度（dispersions）和对数折叠变化（LFCs）。
拟合离散度 ：离散度是基因表达计数变异性的度量，DESeq2 使用一个广义线性模型来估计每个基因的离散度。 计算对数折叠变化（LFCs） ：LFCs 是基因在两个比较条件下表达水平变化的对数比率，是差异表达分析中的关键度量。 dds.deseq2() 结果：
print(dds) AnnData object with n_obs × n_vars = 8 × 21930 obs: 'sample_title' , 'size_factors' , 'replaceable' var: '_normed_means' , 'non_zero' , '_MoM_dispersions' , 'genewise_dispersions' , '_genewise_converged' ,
'fitted_dispersions' , 'MAP_dispersions' , '_MAP_converged' , 'dispersions' , '_outlier_genes' , '_LFC_converged' , 'replaced' , 'refitted' , '_pvalue_cooks_outlier' uns: 'trend_coeffs' , 'disp_function_type' , '_squared_logres' , 'prior_disp_var' obsm: 'design_matrix' , '_mu_LFC' , '_hat_diagonals' varm: 'LFC' layers: 'normed_counts' , '_mu_hat' , 'cooks' # access dispersions print(dds.var[ "dispersions" ]) # LFCs (in natural log scale) print(dds.varm[ "LFC" ]) DeseqDataSet 类通过继承 AnnData 类，使得它能够利用 AnnData 类中定义的数据结构和方法，同时增加额外的功能来支持 DESeq2 分析流程。
X ：主要数据矩阵，通常存储原始的 counts 数据，其中行代表基因，列代表单个样本。 obs ：一个数据框（DataFrame），用于存储与每个样本相关的一维数据，如样本的设计因素（例如，实验条件、时间点等）和标准化所需的大小因子。 obsm ：一个或多个矩阵，用于存储每个样本的多维数据，如设计矩阵，这可能包括连续变量或其他复杂的样本特征。 var ：一个数据框，用于存储与每个基因相关的一维数据，如基因名称和每个基因的离散度估计。 varm ：一个或多个矩阵，用于存储每个基因的多维数据，如对数折叠变化（LFC），这表示基因在不同条件下表达水平的变化。 已经拟合了离散度和对数折叠变化（LFCs），我们可以继续进行统计测试，以计算差异表达的 p 值和调整后的 p 值。这就是 DeseqStats 类的作用。它有两个必需的参数：
dds ，应该是一个拟合好的 DeseqDataSet 对象【必须】， contrast ，是一个包含三个字符串的列表，格式为 ["变量名", "测试水平", "对照水平"] ，或者直接是一个对比向量【必须】。 alpha ：p 值和调整后 p 值的显著性阈值（默认为 0.05）， cooks_filter ：是否基于 Cook's 异常值来过滤 p 值（默认为 True），
independent_filter ：是否执行独立过滤以校正 p 值趋势（默认为 True）。 # 统计检验 ds = DeseqStats(dds, contrast=[ "sample_title" , "Dex" , "untreated" ], inference=inference) # PyDESeq2 使用 Wald 检验来计算 p 值。这可以通过 summary() 方法来完成 ds.summary() 结果随后存储在 results_df 属性中（ ds.results_df ）。
保存结果：
## 提取结果 res = ds.results_df res # 保存 # 保存为 CSV 文件 res.to_csv( 'diff_results.txt' ,sep= "\t" ,index= True ) 绘制火山图就用上次的帖子啦：
python版本的差异表达基因火山图绘制
结果如下：
下期分享一下多分组的分析~
