參考
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來新建一個虛擬環(huán)境沥潭，然后使用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
數(shù)據(jù)來源： 數(shù)據(jù)集是一個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里差不多。
#數(shù)據(jù)讀取
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]
數(shù)據(jù)過濾
在進行DEA之前滋恬，對數(shù)據(jù)進行預(yù)處理是一種好的做法，例如抱究，刪除注釋缺失的樣本恢氯，并排除表達水平非常低的基因。這在我們的合成數(shù)據(jù)的情況下是不必要的媳维，但是如果你使用真實數(shù)據(jù)酿雪，不要忘記這一步。為此侄刽，您可以使用下面的代碼指黎。
首先移除條件為NaN的樣本。如果您正在使用另一個數(shù)據(jù)集州丹，不要忘記更改您希望在分析中用作設(shè)計因素的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]
接下來，我們過濾掉總共閱讀次數(shù)少于10次的基因墓毒。
genes_to_keep = counts_df.columns[counts_df.sum(axis=0) >= 10]
counts_df = counts_df[genes_to_keep]
構(gòu)建DeseqDataSet 對象
# 構(gòu)建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.
統(tǒng)計分析
差異表達統(tǒng)計檢驗分析
res = DeseqStats(dds)
# 執(zhí)行統(tǒng)計分析并返回結(jié)果
res_df = res.summary()
#結(jié)果如下
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
