PyDESeq2使用

参考

GitHub - owkin/PyDESeq2: A Python implementation of the DESeq2 pipeline for bulk RNA-seq DEA.

安装

使用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]

构建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