pydeseq2

Par mkurman · zorai

Analyse d'expression différentielle des gènes (Python DESeq2). Identifie les gènes différentiellement exprimés à partir de comptages RNA-seq en vrac, tests de Wald, correction FDR, graphiques volcano/MA, pour l'analyse RNA-seq.

npx skills add https://github.com/mkurman/zorai --skill pydeseq2

PyDESeq2

Présentation générale

PyDESeq2 est une implémentation Python de DESeq2 pour l'analyse d'expression différentielle avec des données RNA-seq en vrac. Concevez et exécutez des flux de travail complets allant du chargement des données à l'interprétation des résultats, incluant les designs à un seul facteur et multi-facteurs, les tests de Wald avec correction des tests multiples, le rétrécissement optionnel apeGLM et l'intégration avec pandas et AnnData.

Quand utiliser cette compétence

Cette compétence doit être utilisée quand :

  • Analyser des données de comptage RNA-seq en vrac pour l'expression différentielle
  • Comparer l'expression génique entre conditions expérimentales (p. ex., traité vs contrôle)
  • Effectuer des designs multi-facteurs tenant compte des effets de lot ou des covariables
  • Convertir des flux de travail DESeq2 basés sur R en Python
  • Intégrer l'analyse d'expression différentielle dans des pipelines Python
  • Les utilisateurs mentionnent « DESeq2 », « expression différentielle », « analyse RNA-seq » ou « PyDESeq2 »

Flux de travail de démarrage rapide

Pour les utilisateurs qui souhaitent effectuer une analyse d'expression différentielle standard :

import pandas as pd
from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats

# 1. Charger les données
counts_df = pd.read_csv("counts.csv", index_col=0).T  # Transposer en échantillons × gènes
metadata = pd.read_csv("metadata.csv", index_col=0)

# 2. Filtrer les gènes à faible comptage
genes_to_keep = counts_df.columns[counts_df.sum(axis=0) >= 10]
counts_df = counts_df[genes_to_keep]

# 3. Initialiser et ajuster DESeq2
dds = DeseqDataSet(
    counts=counts_df,
    metadata=metadata,
    design="~condition",
    refit_cooks=True
)
dds.deseq2()

# 4. Effectuer des tests statistiques
ds = DeseqStats(dds, contrast=["condition", "treated", "control"])
ds.summary()

# 5. Accéder aux résultats
results = ds.results_df
significant = results[results.padj < 0.05]
print(f"Found {len(significant)} significant genes")

Étapes principales du flux de travail

Étape 1 : Préparation des données

Exigences d'entrée :

  • Matrice de comptage : DataFrame Échantillons × gènes avec des comptages de lectures entiers non-négatifs
  • Métadonnées : DataFrame Échantillons × variables avec des facteurs expérimentaux

Modèles courants de chargement de données :

# À partir d'un CSV (format typique : gènes × échantillons, nécessite une transposition)
counts_df = pd.read_csv("counts.csv", index_col=0).T
metadata = pd.read_csv("metadata.csv", index_col=0)

# À partir d'un TSV
counts_df = pd.read_csv("counts.tsv", sep="\t", index_col=0).T

# À partir d'AnnData
import anndata as ad
adata = ad.read_h5ad("data.h5ad")
counts_df = pd.DataFrame(adata.X, index=adata.obs_names, columns=adata.var_names)
metadata = adata.obs

Filtrage des données :

# Supprimer les gènes à faible comptage
genes_to_keep = counts_df.columns[counts_df.sum(axis=0) >= 10]
counts_df = counts_df[genes_to_keep]

# Supprimer les échantillons avec métadonnées manquantes
samples_to_keep = ~metadata.condition.isna()
counts_df = counts_df.loc[samples_to_keep]
metadata = metadata.loc[samples_to_keep]

Étape 2 : Spécification du design

La formule de design spécifie comment l'expression génique est modélisée.

Designs à un seul facteur :

design = "~condition"  # Simple comparaison entre deux groupes

Designs multi-facteurs :

design = "~batch + condition"  # Contrôler les effets de lot
design = "~age + condition"     # Inclure une covariable continue
design = "~group + condition + group:condition"  # Effets d'interaction

Directives pour la formule de design :

  • Utiliser la notation de formule Wilkinson (style R)
  • Mettre les variables d'ajustement (p. ex., lot) avant la variable d'intérêt principal
  • S'assurer que les variables existent en tant que colonnes dans le DataFrame de métadonnées
  • Utiliser les types de données appropriés (catégorique pour les variables discrètes)

Étape 3 : Ajustement DESeq2

Initialiser le DeseqDataSet et exécuter le pipeline complet :

from pydeseq2.dds import DeseqDataSet

dds = DeseqDataSet(
    counts=counts_df,
    metadata=metadata,
    design="~condition",
    refit_cooks=True,  # Réajuster après suppression des valeurs aberrantes
    n_cpus=1           # Traitement parallèle (ajuster selon les besoins)
)

# Exécuter le pipeline DESeq2 complet
dds.deseq2()

Ce que fait deseq2() :

  1. Calcule les facteurs de taille (normalisation)
  2. Ajuste les dispersions géniques
  3. Ajuste la courbe de tendance de dispersion
  4. Calcule les a priori de dispersion
  5. Ajuste les dispersions MAP (rétrécissement)
  6. Ajuste les changements de repli logarithmiques
  7. Calcule les distances de Cook (détection de valeurs aberrantes)
  8. Réajuste si des valeurs aberrantes sont détectées (optionnel)

Étape 4 : Tests statistiques

Effectuer des tests de Wald pour identifier les gènes exprimés différentiellement :

from pydeseq2.ds import DeseqStats

ds = DeseqStats(
    dds,
    contrast=["condition", "treated", "control"],  # Tester traité vs contrôle
    alpha=0.05,                # Seuil de significativité
    cooks_filter=True,         # Filtrer les valeurs aberrantes
    independent_filter=True    # Filtrer les tests à faible puissance
)

ds.summary()

Spécification du contraste :

  • Format : [variable, test_level, reference_level]
  • Exemple : ["condition", "treated", "control"] teste traité vs contrôle
  • Si None, utilise le dernier coefficient du design

Colonnes du DataFrame de résultats :

  • baseMean : Comptage normalisé moyen entre les échantillons
  • log2FoldChange : Changement de repli Log2 entre les conditions
  • lfcSE : Erreur standard du LFC
  • stat : Statistique du test de Wald
  • pvalue : P-valeur brute
  • padj : P-valeur ajustée (FDR-corrigée via Benjamini-Hochberg)

Étape 5 : Rétrécissement LFC optionnel

Appliquer le rétrécissement pour réduire le bruit dans les estimations de changement de repli :

ds.lfc_shrink()  # Applique le rétrécissement apeGLM

Quand utiliser le rétrécissement LFC :

  • Pour la visualisation (graphiques de volcan, cartes thermiques)
  • Pour classer les gènes par taille d'effet
  • Lors de la priorisation de gènes pour des expériences de suivi

Important : Le rétrécissement affecte uniquement les valeurs log2FoldChange, pas les résultats des tests statistiques (les p-valeurs restent inchangées). Utilisez les valeurs rétrécies pour la visualisation mais signalez les p-valeurs non rétrécies pour la significativité.

Étape 6 : Export des résultats

Enregistrer les résultats et les objets intermédiaires :

import pickle

# Exporter les résultats en CSV
ds.results_df.to_csv("deseq2_results.csv")

# Enregistrer uniquement les gènes significatifs
significant = ds.results_df[ds.results_df.padj < 0.05]
significant.to_csv("significant_genes.csv")

# Enregistrer DeseqDataSet pour une utilisation ultérieure
with open("dds_result.pkl", "wb") as f:
    pickle.dump(dds.to_picklable_anndata(), f)

Modèles d'analyse courants

Comparaison entre deux groupes

Comparaison cas-contrôle standard :

dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~condition")
dds.deseq2()

ds = DeseqStats(dds, contrast=["condition", "treated", "control"])
ds.summary()

results = ds.results_df
significant = results[results.padj < 0.05]

Comparaisons multiples

Tester plusieurs groupes de traitement par rapport au contrôle :

dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~condition")
dds.deseq2()

treatments = ["treatment_A", "treatment_B", "treatment_C"]
all_results = {}

for treatment in treatments:
    ds = DeseqStats(dds, contrast=["condition", treatment, "control"])
    ds.summary()
    all_results[treatment] = ds.results_df

    sig_count = len(ds.results_df[ds.results_df.padj < 0.05])
    print(f"{treatment}: {sig_count} significant genes")

Tenir compte des effets de lot

Contrôler la variation technique :

# Inclure le lot dans le design
dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~batch + condition")
dds.deseq2()

# Tester la condition tout en contrôlant le lot
ds = DeseqStats(dds, contrast=["condition", "treated", "control"])
ds.summary()

Covariables continues

Inclure des variables continues comme l'âge ou la dose :

# S'assurer que la variable continue est numérique
metadata["age"] = pd.to_numeric(metadata["age"])

dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~age + condition")
dds.deseq2()

ds = DeseqStats(dds, contrast=["condition", "treated", "control"])
ds.summary()

Utilisation du script d'analyse

Cette compétence inclut un script en ligne de commande complet pour les analyses standard :

# Utilisation basique
python scripts/run_deseq2_analysis.py \
  --counts counts.csv \
  --metadata metadata.csv \
  --design "~condition" \
  --contrast condition treated control \
  --output results/

# Avec des options supplémentaires
python scripts/run_deseq2_analysis.py \
  --counts counts.csv \
  --metadata metadata.csv \
  --design "~batch + condition" \
  --contrast condition treated control \
  --output results/ \
  --min-counts 10 \
  --alpha 0.05 \
  --n-cpus 4 \
  --plots

Caractéristiques du script :

  • Chargement et validation automatiques des données
  • Filtrage des gènes et des échantillons
  • Exécution complète du pipeline DESeq2
  • Tests statistiques avec paramètres personnalisables
  • Export des résultats (CSV, pickle)
  • Visualisation optionnelle (graphiques de volcan et MA)

Orienter les utilisateurs vers scripts/run_deseq2_analysis.py quand ils ont besoin d'un outil d'analyse autonome ou souhaitent traiter par lot plusieurs ensembles de données.

Interprétation des résultats

Identification des gènes significatifs

# Filtrer par p-valeur ajustée
significant = ds.results_df[ds.results_df.padj < 0.05]

# Filtrer par significativité et taille d'effet
sig_and_large = ds.results_df[
    (ds.results_df.padj < 0.05) &
    (abs(ds.results_df.log2FoldChange) > 1)
]

# Séparer les gènes surexprimés et sous-exprimés
upregulated = significant[significant.log2FoldChange > 0]
downregulated = significant[significant.log2FoldChange < 0]

print(f"Upregulated: {len(upregulated)}")
print(f"Downregulated: {len(downregulated)}")

Classement et tri

# Trier par p-valeur ajustée
top_by_padj = ds.results_df.sort_values("padj").head(20)

# Trier par changement de repli absolu (utiliser les valeurs rétrécies)
ds.lfc_shrink()
ds.results_df["abs_lfc"] = abs(ds.results_df.log2FoldChange)
top_by_lfc = ds.results_df.sort_values("abs_lfc", ascending=False).head(20)

# Trier par une métrique combinée
ds.results_df["score"] = -np.log10(ds.results_df.padj) * abs(ds.results_df.log2FoldChange)
top_combined = ds.results_df.sort_values("score", ascending=False).head(20)

Métriques de qualité

# Vérifier la normalisation (les facteurs de taille doivent être proches de 1)
print("Size factors:", dds.obsm["size_factors"])

# Examiner les estimations de dispersion
import matplotlib.pyplot as plt
plt.hist(dds.varm["dispersions"], bins=50)
plt.xlabel("Dispersion")
plt.ylabel("Frequency")
plt.title("Dispersion Distribution")
plt.show()

# Vérifier la distribution des p-valeurs (doit être principalement plate avec un pic près de 0)
plt.hist(ds.results_df.pvalue.dropna(), bins=50)
plt.xlabel("P-value")
plt.ylabel("Frequency")
plt.title("P-value Distribution")
plt.show()

Directives de visualisation

Graphique de volcan

Visualiser la significativité par rapport à la taille d'effet :

import matplotlib.pyplot as plt
import numpy as np

results = ds.results_df.copy()
results["-log10(padj)"] = -np.log10(results.padj)

plt.figure(figsize=(10, 6))
significant = results.padj < 0.05

plt.scatter(
    results.loc[~significant, "log2FoldChange"],
    results.loc[~significant, "-log10(padj)"],
    alpha=0.3, s=10, c='gray', label='Not significant'
)
plt.scatter(
    results.loc[significant, "log2FoldChange"],
    results.loc[significant, "-log10(padj)"],
    alpha=0.6, s=10, c='red', label='padj < 0.05'
)

plt.axhline(-np.log10(0.05), color='blue', linestyle='--', alpha=0.5)
plt.xlabel("Log2 Fold Change")
plt.ylabel("-Log10(Adjusted P-value)")
plt.title("Volcano Plot")
plt.legend()
plt.savefig("volcano_plot.png", dpi=300)

Graphique MA

Afficher le changement de repli par rapport à l'expression moyenne :

plt.figure(figsize=(10, 6))

plt.scatter(
    np.log10(results.loc[~significant, "baseMean"] + 1),
    results.loc[~significant, "log2FoldChange"],
    alpha=0.3, s=10, c='gray'
)
plt.scatter(
    np.log10(results.loc[significant, "baseMean"] + 1),
    results.loc[significant, "log2FoldChange"],
    alpha=0.6, s=10, c='red'
)

plt.axhline(0, color='blue', linestyle='--', alpha=0.5)
plt.xlabel("Log10(Base Mean + 1)")
plt.ylabel("Log2 Fold Change")
plt.title("MA Plot")
plt.savefig("ma_plot.png", dpi=300)

Résolution des problèmes courants

Problèmes de format de données

Problème : « Index mismatch between counts and metadata »

Solution : S'assurer que les noms des échantillons correspondent exactement

print("Counts samples:", counts_df.index.tolist())
print("Metadata samples:", metadata.index.tolist())

# Prendre l'intersection si nécessaire
common = counts_df.index.intersection(metadata.index)
counts_df = counts_df.loc[common]
metadata = metadata.loc[common]

Problème : « All genes have zero counts »

Solution : Vérifier si les données doivent être transposées

print(f"Counts shape: {counts_df.shape}")
# Si gènes > échantillons, la transposition est nécessaire
if counts_df.shape[1] < counts_df.shape[0]:
    counts_df = counts_df.T

Problèmes de matrice de design

Problème : « Design matrix is not full rank »

Cause : Variables confondues (p. ex., tous les échantillons traités dans un seul lot)

Solution : Supprimer la variable confondante ou ajouter un terme d'interaction

# Vérifier la confusion
print(pd.crosstab(metadata.condition, metadata.batch))

# Soit simplifier le design, soit ajouter une interaction
design = "~condition"  # Supprimer le lot
# OU
design = "~condition + batch + condition:batch"  # Modéliser l'interaction

Aucun gène significatif

Diagnostics :

# Vérifier la distribution de dispersion
plt.hist(dds.varm["dispersions"], bins=50)
plt.show()

# Vérifier les facteurs de taille
print(dds.obsm["size_factors"])

# Regarder les gènes principaux par p-valeur brute
print(ds.results_df.nsmallest(20, "pvalue"))

Causes possibles :

  • Petites tailles d'effet
  • Variabilité biologique élevée
  • Taille d'échantillon insuffisante
  • Problèmes techniques (effets de lot, valeurs aberrantes)

Documentation de référence

Pour des détails complets au-delà de ce guide orienté flux de travail :

  • Référence API (references/api_reference.md) : Documentation complète des classes, méthodes et structures de données de PyDESeq2. À utiliser pour obtenir des informations détaillées sur les paramètres ou pour comprendre les attributs d'objets.

  • Guide de flux de travail (references/workflow_guide.md) : Guide détaillé couvrant les flux de travail d'analyse complets, les modèles de chargement de données, les designs multi-facteurs, la résolution de problèmes et les meilleures pratiques. À utiliser pour gérer les designs expérimentaux complexes ou rencontrer des problèmes.

Charger ces références dans le contexte quand les utilisateurs ont besoin de :

  • Documentation API détaillée : Read references/api_reference.md
  • Exemples de flux de travail complets : Read references/workflow_guide.md
  • Directives de résolution de problèmes : Read references/workflow_guide.md (voir section Troubleshooting)

Points clés à retenir

  1. L'orientation des données compte : Les matrices de comptage se chargent généralement sous la forme gènes × échantillons mais doivent être échantillons × gènes. Toujours transposer avec .T si nécessaire.

  2. Filtrage des échantillons : Supprimer les échantillons avec métadonnées manquantes avant l'analyse pour éviter les erreurs.

  3. Filtrage des gènes : Filtrer les gènes à faible comptage (p. ex., < 10 lectures totales) pour améliorer la puissance et réduire le temps de calcul.

  4. Ordre de la formule de design : Mettre les variables d'ajustement avant la variable d'intérêt (p. ex., "~batch + condition" et non "~condition + batch").

  5. Timing du rétrécissement LFC : Appliquer le rétrécissement après les tests statistiques et uniquement pour la visualisation/classement. Les p-valeurs restent basées sur les estimations non rétrécies.

  6. Interprétation des résultats : Utiliser padj < 0.05 pour la significativité, pas les p-valeurs brutes. La procédure Benjamini-Hochberg contrôle le taux de fausses découvertes.

  7. Spécification du contraste : Le format est [variable, test_level, reference_level] où test_level est comparé à reference_level.

  8. Enregistrer les objets intermédiaires : Utiliser pickle pour enregistrer les objets DeseqDataSet pour une utilisation ultérieure ou d'analyses supplémentaires sans réexécuter l'étape d'ajustement coûteuse.

Installation et exigences

uv pip install pydeseq2

Exigences système :

  • Python 3.10-3.11
  • pandas 1.4.3+
  • numpy 1.23.0+
  • scipy 1.11.0+
  • scikit-learn 1.1.1+
  • anndata 0.8.0+

Optionnel pour la visualisation :

  • matplotlib
  • seaborn

Ressources supplémentaires

Skills similaires

shap
mkurman / zorai

Expliquer les prédictions d'un modèle ML via les valeurs SHAP et visualisations associées.

309
gget
mkurman / zorai

Interroger 20+ bases de données génomiques via une interface Python ou ligne de commande unifiée.

309
scanpy
mkurman / zorai

Analyser des données single-cell RNA-seq avec clustering, visualisation et identification de gènes marqueurs.

309
depmap
mkurman / zorai

Interroger DepMap pour identifier les vulnérabilités génétiques et sensibilités médicamenteuses des cancers.

309
reactome-database
mkurman / zorai

Analyser des voies biologiques et mapper des identifiants via la base Reactome.

309