1. De GWAS à l'annotation fonctionnelle : pourquoi l'analyse de l'héritabilité partitionnée ?

Si vous avez déjà réalisé une étude d'association à l'échelle du génome (GWAS), vous connaissez probablement cette situation : vous terminez votre analyse avec esfuerzo, obtenez une série de loci significatifs, et ensuite ? De nombreux chercheurs s'arrêtent après avoir identifié quelques gènes et publié un article. Cependant, après des années d'analyse en bioinformatique, j'ai découvert que la véritable valeur du GWAS réside souvent dans les loci "non significatifs". C'est précisément là qu'intervient l'analyse de l'héritabilité partitionnée.

En termes simples, l'analyse de l'héritabilité partitionnée vise à répondre à cette question : l'héritabilité de notre phénotype d'intérêt (comme la taille, le risque de maladie) est-elle due à quelles régions fonctionnelles du génome ? Est-ce que les régions codantes contribuent davantage, ou les régions régulatrices ? Les Enhancers actifs dans les neurones sont-ils plus importants, ou les promoteurs dans les cellules immunitaires ? Le GWAS traditionnel peut uniquement indiquer "quels loci sont associés", tandis que l'analyse de l'héritabilité partitionnée révèle "quels types de régions fonctionnelles du génome sont globalement plus importants".

J'ai rencontré cette méthode pour la première fois lors de l'analsye de données de schizophrénie. À l'époque, nous avons découvert que bien que le GWAS ait identifié des centaines de loci, ces loci expliquaient moins de 10% de l'héritabilité totale. Où allait le reste des 90% ? L'analyse partitionnée avec LDSC a révélé qu'une grande partie de l'héritabilité était enrichie dans les régions d'enhancers conservées, particulièrement dans les régions de chromatine ouverte spécifiques aux neurones. Cette découverte a directement influencé la direction des expériences suivantes.

LDSC (Linkage Disequilibrium Score Regression) présente l'avantage de ne pas nécessiter de données génotypiques au niveau individuel, fonctionnant uniquement avec les statistiques sommaires du GWAS. Pour de nombreuses études où seules les données sommaires sont disponibles, c'est véritablement une bénédiction. De plus, cette méthode prend en compte le déséquilibre de liaison (LD), évitant les faux positifs liés au LD.

2. Préparation avant la pratique : données et environnement

2.1 Comprendre les données nécessaires

Avant de commencer, il est crucial de comprendre quelles données sont requises. D'après mon expérience, les débutants commettent souvent l'erreur de ne pas préparer complètement les données. Vous aurez besoin de trois éléments :

Le premier est le fichier de scores LD de la population de référence. LDSC fournit des scores LD précalculés basés sur le projet 1000 Genomes, différenciés par population. Si vos données GWAS concernent une population d'Asie de l'Est, utilisez la version EAS (East Asian) ; si c'est une population européenne, utliisez la version EUR. Utiliser les données de référence d'une population différente peut complètement fausser les résultats.

Le second élément est les statistiques sommaires GWAS, généralement au format .sumstats.gz. Ce fichier doit contenir des informations de base comme l'identifiant SNP, l'allèle effet, l'allèle non-effet, la valeur P, la taille d'échantillon, etc. Si le format n'est pas adapté, vous devrez le traiter avec le script munge_sumstats.py fourni avec LDSC.

Le troisième élément est le fichier d'annotations fonctionnelles. LDSC propose un modèle de base (baseline model) contenant 53 catégories d'annotations fonctionnelles, comme les régions codantes, UTR, promoteurs, enhancers, diverses régions de modifications des histones, etc. Vous pouvez également ajouter vos propres annotations, comme les pics ATAC-seq spécifiques à un type cellulaire.

2.2 Procédure de téléchargement des fichiers

Laissez-moi vous guider à travers le processus complet de téléchargement. Supposons que vos données GWAS concernent une population d'Asie de l'Est, nous utiliserons donc les données EAS :

# Créer le répertoire du projet et y accéder
mkdir -p analyse_heritabilite_partitionnee
cd analyse_heritabilite_partitionnee

# Télécharger le fichier de scores LD de base pour population EAS
echo "Téléchargement du fichier de scores LD de base..."
wget -c https://storage.googleapis.com/broad-alkesgroup-public/LDSCORE/1000G_Phase3_EAS_baseline_v1.2_ldscores.tgz

# Télécharger les fichiers plink pour population EAS (pour le calcul des fréquences)
echo "Téléchargement des fichiers plink..."
wget -c https://storage.googleapis.com/broad-alkesgroup-public/LDSCORE/1000G_Phase3_EAS_plinkfiles.tgz

# Télécharger le fichier de poids (liste SNP HM3, sans région MHC)
echo "Téléchargement du fichier de poids..."
wget -c https://storage.googleapis.com/broad-alkesgroup-public/LDSCORE/1000G_Phase3_EAS_weights_hm3_no_MHC.tgz

# Décompresser tous les fichiers
echo "Décompression en cours..."
tar -zxvf 1000G_Phase3_EAS_baseline_v1.2_ldscores.tgz
tar -zxvf 1000G_Phase3_EAS_plinkfiles.tgz  
tar -zxvf 1000G_Phase3_EAS_weights_hm3_no_MHC.tgz

# Vérifier la structure des fichiers
echo "Vérification des fichiers téléchargés..."
ls -lh 1000G_Phase3_EAS_baseline_v1.2_ldscores/

Après le téléchargement, vous devriez voir une structure de répertoire similaire à :

1000G_Phase3_EAS_baseline_v1.2_ldscores/
├── baseline.1.l2.ldscore.gz
├── baseline.1.annot.gz
├── baseline.2.l2.ldscore.gz
├── baseline.2.annot.gz
...
└── baseline.22.l2.ldscore.gz

Chaque chromosome possède son fichier .l2.ldscore.gz (fichier de scores LD) et .annot.gz (fichier d'annotations). Le préfixe baseline. correspond au préfixe du fichier, et le chiffre suivant représente le numéro de chromosome.

2.3 Configuration de l'environnement et installation de LDSC

L'installation de LDSC est relativement simple, mais certaines dépendances méritent attention. Je recommande d'utiliser conda pour gérer l'environnement, ce qui permet d'éviter les conflits de versions :

# Cloner le dépôt LDSC
git clone https://github.com/bulik/ldsc.git
cd ldsc

# Créer et activer l'environnement conda
conda env create -f environment.yml
conda activate ldsc

# Vérifier si l'installation est réussie
python ldsc.py -h

Si vous voyez s'afficher l'aide, l'installation est réussie. Cependant, je dois vous avertir que des problèmes de version Python peuvent survenir dans la pratique. LDSC recommande officiellement Python 2.7, mais de nombreux systèmes utilisent Python 3 par défaut. J'ai testé avec Python 3.7+ et cela fonctionne églaement, bien qu'il puisse être nécessaire d'ajuster légèrement le code. Si vous rencontrez des erreurs concernant les instructions print, c'est généralement dû aux différences de syntaxe entre Python 2 et 3, et quelques modifications simples suffisent pour résoudre le problème.

Un autre problème fréquent est la mémoire insuffisante. LDSC требуbeaucoup de mémoire lors du traitement de grandes données, surtout lorsque votre GWAS contient des millions de SNPs. Je recommande d'exécuter sur un serveur, en allouant au moins 8-16 Go de mémoire. Si vous travaillez sur un ordinateur personnel, envisagez d'utiliser le paramètre --chunksize pour traiter par morceaux.

3. Préparation des statistiques sommaires GWAS : conversion de format et contrôle qualité

3.1 Comprendre les exigences de format sumstats

LDSC a des exigences strictes concernant le format du fichier d'entrée. Vos statistiques sommaires GWAS doivent être converties dans un format sumstats spécifique. Laissez-moi vous expliquer en détail la signification de chaque champ :

• SNP : L'identifiant SNP, par exemple "rs123456"
• A1 : L'allèle effet (effect allele), correspondant généralement à la direction de l'effet dans le GWAS
• A2 : L'allèle non-effet (other allele)
• N : La taille d'échantillon, élément crucial utilisé par LDSC pour la pondération
• Z : Le score Z, calculable à partir de la valeur P et de la valeur beta
• P : La valeur P (facultatif mais recommandé)

Pourquoi avoir besoin du score Z plutôt que directement beta et SE ? Car le modèle central de LDSC est basé sur la statistique du χ², et le carré du score Z correspond au χ². De plus, les unités d'effet peuvent varier selon les études, et le score Z normalisé est plus comparable.

3.2 Utilisation de munge_sumstats.py pour la conversion de format

La plupart des résultats GWAS ne sont pas directement utilisables au format requis, c'est là que le script munge_sumstats.py entre en jeu. Ce script peut traiter divers formats d'entrée, reconnaître automatiquement les noms de colonnes et effectuer les contrôles qualité nécessaires.

Supposez que vous possédez un fichier de résultats GWAS typique mon_gwas.txt contenant ces colonnes : SNP, CHR, BP, A1, A2, BETA, SE, P, N. La commande de conversion est la suivante :

python munge_sumstats.py \
    --sumstats mon_gwas.txt \
    --out ma_trait \
    --merge-alleles w_hm3.snplist \
    --a1 A1 \
    --a2 A2 \
    --snp SNP \
    --p P \
    --frq freq \
    --N-col N \
    --signed-sumstats BETA,SE \
    --chunksize 500000

Voici l'explication des paramètres clés :

• --merge-alleles w_hm3.snplist : Conserver uniquement les SNPs HapMap3, pratique standard LDSC car les SNPs HM3 sont de haute qualité et bien couverts
• --chunksize 500000 : Traitement par blocs pour éviter les débordements mémoire
• --signed-sumstats BETA,SE : Indiquer au script d'utiliser beta et SE pour calculer le score Z

Après exécution, le script génère un rapport détaillé de contrôle qualité, par exemple :

Read 2,500,000 SNPs from --sumstats file.
Removed 800,000 SNPs not in --merge-alleles.
Removed 50,000 SNPs with missing values.
Removed 20,000 SNPs with INFO <= 0.9.
Removed 30,000 SNPs with MAF <= 0.01.
1,600,000 SNPs remain.

3.3 Gestion des problèmes de format courants

Dans la pratique, j'ai rencontré divers problèmes de format étranges. Une fois, un collaborateur a fourni un fichier où les colonnes d'allèles contenaient des représentations "I/D" (insertion/délétion) au lieu des bases standard "A/C/G/T". Une autre fois, la colonne de taille d'échantillon contenait des valeurs NA sur certaines lignes. Ces problèmes peuvent tous provoquer l'échec de la conversion.

Mon expérience m'a appris : commencer par un test à petite échelle. Utilisez head -1000 pour extraire les 1000 premières lignes afin de tester la conversion, et traitez le fichier complet uniquement après vérification. De plus, consultez toujours le fichier .log généré, qui enregistre en détail le traitement de chaque étape.

Un autre problème fréquent est l'ambiguïté de la chaîne (strand). Certains SNPs sont A/T ou C/G, ces paires de bases symétriques ne permettent pas de déterminer le sens positif ou négatif. LDSC filtre par défaut ces SNPs à chaîne ambiguë. Si vous constatez que le nombre de SNPs filtrés est particulièrement élevé, vous devrez peut-être vérifier les informations de chaîne des données originales.

4. Exécution de l'analyse de l'héritabilité partitionnée : explanation détaillée des commandes et optimisation des paramètres

4.1 Analyse de la commande de base

Une fois les données préparées, vous pouvez exécuter l'analyse de l'héritabilité partitionnée. La structure de la commande de base est la suivante :

conda activate ldsc
python ldsc.py \
    --h2 scz.sumstats.gz \
    --ref-ld-chr /chemin/vers/1000G_Phase3_EAS_baseline_v1.2_ldscores/baseline. \
    --w-ld-chr /chemin/vers/1000G_Phase3_EAS_weights_hm3_no_MHC/weights.EAS.hm3_noMHC. \
    --overlap-annot \
    --print-coefficients \
    --frqfile-chr /chemin/vers/1000G_Phase3_EAS_plinkfiles/1000G.EAS.QC. \
    --out scz_baseline

Laissez-moi vous expliquer chaque paramètre :

• --h2 : Spécifie votre fichier de statistiques sommaires GWAS, indiquant à LDSC de calculer l'héritabilité
• --ref-ld-chr : Fichier de scores LD de référence, le fichier du modèle de base téléchargé précédemment. Notez que le chemin doit terminer par baseline., LDSC ajoutera automatiquement le numéro de chromosome
• --w-ld-chr : Fichier de poids de régression, utilisé pour corriger la structure LD