参考
GitHub - owkin/PyDESeq2: A Python implementation of the DESeq2 pipeline for bulk RNA-seq DEA.
A simple PyDESeq2 workflow — PyDESeq2 0.2.1 documentation
安装
使用conda来新建一个虚拟环境，然后使用pip安装。
conda create -n pydeseq2 python=3.8
conda activate pydeseq2
pip install pydeseq2
查看是否安装成功
(pydeseq2) [zhaoyuhu@localhost ~]$ python
Python 3.8.15 | packaged by conda-forge | (default, Nov 22 2022, 08:49:35)
[GCC 10.4.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
#加载pydeseq2包,如果安装成功的话，是可以被成功加载的
from pydeseq2.DeseqDataSet import DeseqDataSet
from pydeseq2.DeseqStats import DeseqStats
import numpy as np
import pandas as pd
用法
参考链接：PyDESeq2/test_pydeseq2.py at main · owkin/PyDESeq2 · GitHub
数据来源： 数据集是一个100个样本，每个样本10个基因的小测试集。而其中50个样本属于条件A,另50个样本属于条件B。
①test_counts.csv 文件：https://github.com/owkin/PyDESeq2/blob/main/tests/data/test_counts.csv
②test_clinical.csv文件：https://github.com/owkin/PyDESeq2/blob/main/tests/data/test_clinical.csv
作为Python版的DESeq2， 用法和R里差不多。
#数据读取
counts_df = pd.read_csv("test_counts.csv", index_col=0).T
condition_df = pd.read_csv( "test_clinical.csv", index_col=0)
>>> counts_df.shape
(100, 10)
>>> counts_df.head()
gene1 gene2 gene3 gene4 gene5 gene6 gene7 gene8 gene9 gene10
sample1 12 22 2 187 15 2 13 57 56 6
sample2 10 6 20 99 55 0 35 96 43 1
sample3 0 28 3 96 38 2 9 54 27 14
sample4 7 28 10 170 16 10 17 38 18 16
sample5 2 31 5 126 23 2 19 53 31 18
>>> condition_df
condition
sample1 A
sample2 A
sample3 A
sample4 A
sample5 A
... ...
sample96 B
sample97 B
sample98 B
sample99 B
sample100 B
[100 rows x 1 columns]
数据过滤
在进行DEA之前，对数据进行预处理是一种好的做法，例如，删除注释缺失的样本，并排除表达水平非常低的基因。这在我们的合成数据的情况下是不必要的，但是如果你使用真实数据，不要忘记这一步。为此，您可以使用下面的代码。
首先移除条件为NaN的样本。如果您正在使用另一个数据集，不要忘记更改您希望在分析中用作设计因素的clinical_df列的“条件”。
samples_to_keep = ~clinical_df.condition.isna()
counts_df = counts_df.loc[samples_to_keep]
clinical_df = clinical_df.loc[samples_to_keep]
接下来，我们过滤掉总共阅读次数少于10次的基因。
genes_to_keep = counts_df.columns[counts_df.sum(axis=0) >= 10]
counts_df = counts_df[genes_to_keep]
构建DeseqDataSet 对象
# 构建DeseqDataSet 对象
dds = DeseqDataSet(counts_df, condition_df, design_factor="condition")
# 离散度和log fold-change评估.
dds.deseq2()
#Fitting size factors...
#... done in 0.02 seconds.
#Fitting dispersions...
#... done in 0.68 seconds.
#Fitting dispersion trend curve...
#... done in 0.16 seconds.
#Fitting MAP dispersions...
#... done in 0.77 seconds.
#Fitting LFCs...
#... done in 0.71 seconds.
#Refitting 0 outliers.
统计分析
差异表达统计检验分析
res = DeseqStats(dds)
# 执行统计分析并返回结果
res_df = res.summary()
#结果如下
res_df
baseMean log2FoldChange lfcSE stat pvalue padj
gene1 10.306788 1.007045 0.225231 4.471161 7.779603e-06 2.593201e-05
gene2 24.718815 -0.059670 0.165606 -0.360311 7.186146e-01 7.186146e-01
gene3 4.348135 -0.166592 0.325445 -0.511891 6.087275e-01 6.763639e-01
gene4 98.572300 -2.529204 0.136752 -18.494817 2.273125e-76 2.273125e-75
gene5 38.008562 1.236663 0.151824 8.145377 3.781028e-16 1.890514e-15
gene6 4.734285 0.212656 0.304487 0.698408 4.849222e-01 6.061527e-01
gene7 30.011855 -0.445855 0.150575 -2.961017 3.066249e-03 5.110415e-03
gene8 59.330642 0.372080 0.118911 3.129070 1.753603e-03 3.507207e-03
gene9 46.779546 0.547280 0.124922 4.380966 1.181541e-05 2.953853e-05
gene10 11.963156 0.494775 0.229494 2.155940 3.108836e-02 4.441194e-02
