Analyse transcriptomique basée sur des données de RNA-seq en vrac

Summary

Le protocole actuel établit un pipeline complet pour analyser le processus de séquence ARN en vrac depuis les données brutes jusqu’à l’analyse d’enrichissement fonctionnel.

Abstract

Le foie gras non alcoolique (NAFL) est généralement considéré comme une affection bénigne ; cependant, une fois qu’elle évolue vers une stéatohépatite non alcoolique (NASH), les patients présentent un risque significativement accru de développer une maladie hépatique terminale. De nombreuses études tentent d’élucider le mécanisme moléculaire sous-jacent à la transition de NAFL à NASH. Les technologies de séquençage à haut débit (telles que l’ARN-seq en vrac) ont permis aux chercheurs d’approfondir leur compréhension en examinant le transcriptome, en révélant l’expression des molécules, l’activation des voies de signalisation et d’autres facteurs liés à la progression de la maladie. Il existe une mine de données open source à disposition pour les chercheurs afin d’identifier des cibles potentielles pour le traitement des maladies. Cependant, la recherche associée est limitée par l’absence d’un processus efficace et fiable pour l’analyse en amont du transcriptome. Ici, une chaîne d’analyse en amont hautement reproductible et conviviale ainsi qu’une analyse génétique différentielle associée sont fournies pour permettre un traitement standardisé et un analyse approfondie des données privées ou publiques. Le pipeline est divisé en quatre étapes : (1) contrôle qualité des données ; (2) cartographie génétique ; (3) analyse différentielle des gènes ; et (4) analyse fonctionnelle. Ce processus vise à découvrir les mécanismes moléculaires de la transformation des maladies et à aider les chercheurs à dépister les cibles potentielles des médicaments et les approches thérapeutiques grâce à l’analyse des données Bulk RNA-seq.

Introduction

La stéatose hépatique non alcoolique (NAFLD) est la maladie chronique du foie la plus répandue au monde, touchant plus d’un quart de la population. Son incidence a augmenté de façon spectaculaire ces dernières^{décennies 1, 2, 3}. La charge croissante des maladies, en particulier sa forme plus avancée, la stéatohépatite non alcoolique (NASH), représente un défi sanitaire mondial majeur et un lourd fardeau économique⁴. Le premier stade de la NAFLD est le foie gras non alcoolique (NAFL), accompagné d’inflammation et de fibrose pouvant évoluer vers la NASH. Cette dernière augmentation significative le risque de progression vers une maladie hépatique terminale, incluant la cirrhose et le carcinome hépatocellulaire (HCC^)5,6,7. L’incidence et la mortalité du CHC sont associées à une augmentation de^{NASH 8,9}, et il est prévu que NAFLD/NASH deviendra l’indicateur principal de la transplantation hépatique d’ici 203010^. Cependant, la progression clinique de la NAFLD est très^{hétérogène 11}, ce qui entrave gravement le développement des médicaments^{pertinents 12}, rendant particulièrement important d’explorer précisément les mécanismes moléculaires impliqués.

L’acquisition massive d’informations compositionnelles basées sur l’ARN-seq peut élucider de manière significative la pathogenèse de diverses maladies. Ces dernières décennies, de nombreuses études en vrac sur l’ARN-seq ont été menées chez des organismes modèles et des humains afin d’élucider les différences d’expression génique dans la progression^{13, 14, 15} de ^NASH, afin d’identifier de nouvelles cibles thérapeutiques pour l’intervention. Sur la base d’une analyse globale de l’ARN-seq, Xiong et al. ont constaté que les cellules non parenchymateuses (NPC) du foie sont impliquées dans des processus tels que la formation de matrices extracellulaires et l’adhésion cellulaire, qui contribuent à la progression de NASH¹⁶. Li et al. ont démontré que la protéine associante à la tumeur de Wilms hépatique (WTAP) dans les hépatocytes régule l’accumulation et l’inflammation des lipides ectopiques, favorisant ainsi la formation^{de NASH 17}. Bien que l’analyse globale d’ARN-seq soit un outil puissant pour élucider les mécanismes de la NASH, ses résultats sont très sensibles à la qualité des données en amont. L’hétérogénéité des opérations expérimentales en amont et des processus d’analyse peut sérieusement nuire à la fiabilité des données, masquant ainsi les véritables informations biologiques et nuisant à la précision des analyses ultérieures. Il est donc important d’établir un ensemble de procédures standardisées d’analyse en amont.

Comparé au séquençage à ARN unicellulaire (scRNA-seq), le séquençage ARN en vrac offre plusieurs avantages distincts tant dans la conception expérimentale que dans les applications pratiques. Bien que scRNA-seq permette d’identifier l’hétérogénéité cellulaire au niveau de la cellule unique et d’analyser précisément les caractéristiques transcriptionnelles spécifiques à chaque type cellulaire, il est associé à un coût élevé, des exigences complexes de traitement des données, et une sensibilité limitée à la détection de transcrits à faible^{abondance 18}. En revanche, le séq ARN en vrac offre une profondeur de séquençage plus élevée, un coût moindre et un débit d’échantillons plus élevé, ce qui le rend particulièrement adapté aux analyses d’expression différentielle génique au niveau de la population et à l’exploration des mécanismes^{moléculaires 19}. Ainsi, guidé par des flux de travail analytiques standardisés, le séquence ARN en vrac reste une approche efficace, rentable et robuste pour étudier la base moléculaire des maladies complexes.

Ce protocole est conçu spécifiquement pour les ensembles de données ARN-seq en vrac dérivés de tissus humains avec une intégrité ARN élevée (RIN ≥ 7,0) et suffisamment d’ARN d’entrée (≥ 500 ng par échantillon). Pour garantir une exécution fiable des étapes d’alignement et de quantification, il est recommandé d’avoir une station de travail locale équipée d’au moins un processeur de 10 cœurs, 32 Go de RAM et un minimum de 200 Go d’espace disque libre. S’appuyant sur ces exigences, le protocole offre un flux de travail analytique efficace et convivial, incluant des instructions opérationnelles détaillées et des configurations de paramètres standardisées, afin de répondre aux besoins des chercheurs analysant des données transcriptomiques à grande échelle.

Protocol

À des fins de démonstration, le jeu de données public PRJNA1023502 généré par Lan Bai et al. a été utilisé pour illustrer chaque étape des analyses amont et aval²⁰. Comme ce jeu de données provient de la base de données en libre accès NCBI SRA, aucune autorisation supplémentaire ni approbation éthique n’est requise. Consultez le tableau des matériaux pour vérifier toutes les versions requises des logiciels et des paquets R. L’ensemble de données disponible publiquement PRJNA1023502 comprend 6 échantillons non NASH, 6 NAFL et 6 échantillons d’ARN-seq hépatique NASH. Dans ce protocole, l’ensemble de données a été utilisé pour démontrer toutes les étapes du flux de travail en masse RNA-seq, y compris la récupération des données à partir de la base de données SRA, le contrôle qualité (fastp), l’alignement (HISAT2), la quantification (featureCounts), ainsi que les analyses d’expression différentielle et d’enrichissement fonctionnel en aval.

1. Installation de la boîte à outils SRA

1. Visitez le site officiel du kit à outils de la SRA et téléchargez la version 3.2.1.

2. Téléchargement de données publiques

1. Obtenez le numéro SRA.
   1. Accédez-y dans l’article complémentaire²⁰, section Disponibilité des données, ou en recherchant des mots-clés dans la base de données NCBI SRA.
   2. Tapez le préfetch dans le terminal pour le télécharger.

3. Génération de la matrice de décompte des gènes

1. Définissez les chemins et les paramètres comme décrit ci-dessous :
   REFERENCE=~/reference/human/GRCh38/GRCh38.primary_assembly.genome.fa
GTF=~/reference/human/GRCh38/gencode.v44.annotation.gtf
INDEX=~/reference/human/GRCh38/GRCh38_index
FASTQ_DIR=~/SRA_tutorial/fastq
OUT_FASTP=~/RNAseq/fastp
OUT_HISAT2=~/RNAseq/hisat2
OUT_COUNTS=~/RNAseq/counts
​mkdir -p $FASTQ_DIR $OUT_FASTP $OUT_HISAT2 $OUT_COUNTS
   1. Téléchargez le génome de référence humain (assemblage primaire GRCh38) et le fichier d’annotation génétique correspondant (gencode.v44, ensemble complet d’annotation génétique) depuis la base de données officielle GENCODE (https://www.gencodegenes.org/human/).
      for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
​for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
2. Renommer les fichiers SRA.
   1. Assurez-vous que tous les fichiers SRA utilisent l’extension .sra pour faciliter la reconnaissance et le traitement précis par les outils en aval.
   2. Exécutez les commandes suivantes :
      for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
3. Exécutez la commande suivante pour convertir les fichiers SRA en format FASTQ :
   for f in *.sra; do fasterq-dump "$f" --split-files -O $FASTQ_DIR - e 20; done
4. Construis l’indice HISAT2 à partir du génome de référence.
   hisat2-build $REFERENCE $INDEX
5. Traiter les fichiers FASTQ avec un suivi visuel de l’avancement.
   for fq in $FASTQ_DIR/*.fastq; do
​sample=$(basename "$fq" .fastq)
   1. Pour les données de séquençage à extrémités jumelées, exécutez la commande suivante
      for fq1 in $FASTQ_DIR/*_1.fastq; do
sample=$(basename "$fq1" _1.fastq)
​fq2=$FASTQ_DIR/${sample}_2.fastq
   2. Effectuez le contrôle qualité et le filtrage à l’aide de fastp et envoyez les résultats dans le répertoire $OUT_FASTP. Pour les données de séquençage à une seule extrémité, exécutez la commande suivante :
      fastp \
-i "${fq}" \
-o $OUT_FASTP/${sample}.clean.fastq \
-h $OUT_FASTP/${sample}.html \
-j $OUT_FASTP/${sample}.json \
-w 20
      Pour les données de séquençage par bout, exécutez la commande suivante :
      fastp \
-i "${fq}" \ -I "$fq2" \
-o $OUT_FASTP/${sample}_1.clean.fastq \
-O $OUT_FASTP/${sample}_2.clean.fastq \
-h $OUT_FASTP/${sample}.html \
-j $OUT_FASTP/${sample}.json \
​-w 20
   3. Pour chaque échantillon, inspectez le rapport HTML pour le contrôle qualité visuelle et consultez le rapport JSON pour des métriques structurées, incluant la qualité de lecture, la distribution du GC, la redondance de lecture et la contamination de l’adaptateur. Assurez-vous que tous les échantillons respectent les seuils de qualité suivants : Q30 ≥ 85 %, taux de contamination de l’adaptateur < 5 %, et aucun schéma GC anormal. Stockez tous les fichiers de sortie dans le dossier $OUT_FASTP.
   4. Effectuez l’alignement des séquences à l’aide de HISAT2, générez des fichiers SAM, puis envoyez-les dans le répertoire $OUT_HISAT2. Pour les données de séquençage à une seule extrémité, exécutez la commande suivante :
      hisat2 -p 20 \ -x $INDEX \-U $OUT_FASTP/${sample}.clean.fastq \
-S $OUT_HISAT2/${sample}.sam
      Pour les données de séquençage par bout, exécutez la commande suivante :
      hisat2 -p 20 \-x $INDEX \-1 $OUT_FASTP/${sample}_1.clean.fastq \
-2 $OUT_FASTP/${sample}_2.clean.fastq \
​-S $OUT_HISAT2/${sample}.sam
   5. Effectuez l’alignement HISAT2 des lectures filtrées par qualité avec le génome de référence humain, en utilisant un index construit à partir du fichier FASTA de l’assemblage primaire GRCh38. Fournir des commandes exemples pour les données de séquençage à extrémité unique et à extrémité appariée.
   6. Convertir SAM en BAM, trier et indexer.
      samtools view -@ 20 -bS $OUT_HISAT2/${sample}.sam \
| samtools sort -@ 20 -o $OUT_HISAT2/${sample}.sorted.bam
samtools index $OUT_HISAT2/${sample}.sorted.bam
​done
   7. Générez un fichier SAM pour chaque échantillon et convertissez-le en un fichier BAM trié et indexé à l’aide de SAMtools. Pour des échantillons transcriptomiques humains de haute qualité (nombre d’intégrité ARN, RIN ≥ 8,0), assurez-vous que le taux d’alignement global dépasse 85 % ; pour les échantillons ARN-seq typiques en vrac (RIN ≥ 7,0), considérez les taux d’alignement ≥ 70 % acceptables.
   8. Effectuez la quantification des gènes à l’aide de featureCount.
      featureCounts -T 20 -p -s 0 \
-a $GTF \
-o $OUT_COUNTS /${sample}.counts.txt \
$OUT_HISAT2/${sample}.sorted.bam
Done
   9. Vérifiez le fichier de sortie délimité par tabulation (*.counts.txt) et son rapport de résumé (*.counts.txt.summary) généré par les featureCounts pour chaque échantillon. S’assurer que le taux d’affectation de lecture atteint le seuil typique de ≥70 % pour l’ARN-seq en vrac humain ; Un taux nettement plus bas peut indiquer une orientation des brins déséquilibrée, des problèmes d’annotation ou une mauvaise qualité d’alignement. Utilisez le paramètre -s 0 pour ce jeu de données ARN-seq non spécifique aux brins. Pour les bibliothèques spécifiques à chaque bris, substituez -s 0 par -s 1 ou -s 2 dans la commande.
6. Générez la matrice de comptage des gènes comme décrit ci-dessous.
   1. Initialiser la matrice de comptage avec les identifiants des gènes et les comptages du premier échantillon.
      cut -f1 $(ls $OUT_COUNTS/*.counts.txt | head -1) > all_counts.txt
   2. Additionnez itérativement les comptes de chaque échantillon dans une seule matrice.
      for f in $OUT_COUNTS/*.counts.txt; do
cut -f7 "$f" | paste all_counts.txt - > tmp && mv tmp
all_counts.txt
​done
   3. Ajoutez une ligne d’en-tête avec des identifiants d’échantillons à la matrice des comptages.
      samples=$(ls *.counts.txt | sed 's/.counts.txt//' | paste -sd "\t")
echo -e "Geneid\t$samples" | cat - all_counts.txt > counts_matrix.txt
   4. Extraire les longueurs des gènes du fichier GTF (somme des longueurs des exons par gène).
      awk '$3=="exon"{match($0,/gene_id "([^"]+)"/,a); if(a[1]!=""){len=$5-$4+1; gene_len[a[1]]+=len}} END{print "GENE_ID\tLENGTH"; for(g in gene_len) print g"\t"gene_len[g]}' \$GTF > gene_length.txt

4. Traitement brut de la matrice de comptage et annotation des gènes

1. Initialisation et chargement de l’environnement R pertinent.
   1. Utilisez la fonction install.packages() pour installer tidyverse, ggplot2, ggrepel, RColorBrewer, ggridges, FactoMineR. Le package tidyverse est utilisé pour la manipulation et le tracé de données ; Le package ggplot2 est utilisé pour la visualisation ; Le package ggrepel est utilisé pour les étiquettes de texte non chevauchantes ; le package RColorBrewer est utilisé pour les palettes de couleurs ; Le package GGRIDGES est utilisé pour les parcelles de crête ; le package FactoMineR est utilisé pour l’analyse PCA et multivariée.
   2. Utilisez la fonction BiocManager ::install() pour installer biomaRt, DESeq2, clusterProfiler. Le package biomaRt est utilisé pour l’annotation génétique d’Ensembl ; le package DESeq2 est utilisé pour l’analyse d’expression différentielle ; le package clusterProfiler est utilisé pour l’analyse d’enrichissement fonctionnel.
   3. Lisez la matrice de comptage originale générée à l’étape 2, utilisez biomaRt pour associer Ensembl ID au nom du gène HGNC pour une analyse ultérieure, retirez la colonne Geneid dupliquée, et générez la matrice de comptage nettoyée. Utilisez la matrice de comptage originale (counts_matrix.csv) en entrée, avec les identifiants de gènes Ensembl comme lignes et les échantillons comme colonnes. Sauvegardez la matrice des comptes nettoyée en sortie (clean_counts_SRA.csv), avec les symboles du gène HGNC comme lignes et les échantillons comme colonnes.
      mart <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
id_map <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = exprSet$GeneID,
mart = mart)
exprSet <- exprSet %>%
left_join(id_map, by = c("GeneID" = "ensembl_gene_id")) %>%
filter(!is.na(hgnc_symbol), hgnc_symbol != "") %>%
distinct(hgnc_symbol, .keep_all = TRUE) %>%
column_to_rownames("hgnc_symbol")

5. Quantification de l’expression génique

REMARQUE : Consultez le Fichier Supplémentaire 1 pour le script détaillé.

1. Exécutez la commande suivante pour calculer les fragments par kilobase de transcription par million de lectures cartographiées (FPKM) et les transcriptions par million (TPM).
   counts <- read.csv("output/clean_counts_SRA.csv", header=TRUE, row.names=1)
gene_len <- read.delim("data/gene_length.txt", header=FALSE, col.names=c("gene_symbol","length"))
gene_len <- gene_len %>% distinct(gene_symbol, .keep_all=TRUE)
rownames(gene_len) <- gene_len$gene_symbol
gene_len <- gene_len[match(rownames(counts), gene_len$gene_symbol),]
length_bp <- gene_len$length
fpkm <- (counts / length_bp) * 1e9 / colSums(counts)
write.csv(fpkm, "output/clean_fpkm_SRA.csv")
tpm <- (counts / length_bp) / colSums(counts / length_bp) * 1e6
write.csv(tpm, "output/clean_tpm_SRA.csv")

6. Regroupement d’échantillons et visualisation des différences

1. Évaluez la qualité des données en examinant le regroupement d’échantillons via PCA. Consultez le fichier supplémentaire 1 pour le script détaillé.
2. Pour la visualisation PCA, utilisez le code ci-dessous.
   gene.pca <- PCA(exprSet, ncp = 2, scale.unit = TRUE, graph = FALSE)
ggplot(pca_sample, aes(x = Dim.1, y = Dim.2)) +
geom_point(aes(color = group)) +
labs(x = paste('PC1:', pca_eig1, '%'),
y = paste('PC2:', pca_eig2, '%'))

7. Analyse différentielle de l’expression et visualisation des résultats

REMARQUE : Consultez le Fichier Supplémentaire 1 pour le script détaillé.

1. Exécutez les commandes suivantes pour construire le jeu de données DESeq2 et fixer le seuil. Seuls les gènes avec un total de décomptes > 1 sur tous les échantillons ont été conservés avant la modélisation. Identifier les gènes différenciellement exprimés (DEG) en utilisant les seuils de signification des valeurs p ajustées < 0,05 et |log2FC| > 0,5.
   dds <- DESeq(DESeqDataSetFromMatrix(countData = exprSet, colData = colData, design = ~group)); sizeFactors(dds); res <- results(dds); dds <- dds[rowSums(counts(dds)) > 1,]
dd1 <- results(dds, contrast = contrast, alpha = 0.05)
dd2 <- lfcShrink(dds, contrast = contrast, res = dd1, type = "ashr")
2. Visualisez des gènes exprimés différemment avec un graphique volcanique.
   ggplot(data = data, aes(x = log2FoldChange, y = -log10(padj))) +
geom_point(aes(color = group), alpha = 1, size = 1.2) +
geom_hline(yintercept = -log10(0.05), lty = 4) +
geom_vline(xintercept = c(-0.5, 0.5), lty = 4) +
geom_text_repel(data = subset(data, abs(log2FoldChange) >= 1.5 & padj < 0.05),
aes(label = gene_id))

8. Réaliser une analyse et une visualisation d’enrichissement fonctionnel

REMARQUE : Consultez le Fichier Supplémentaire 1 pour le script détaillé.

1. Visualisez les résultats d’enrichissement KEGG.
   EGG <- enrichKEGG(gene = gene$ENTREZID, organism = 'hsa',
pvalueCutoff = 0.05, qvalueCutoff = 0.05)
ggplot(symboldata, aes(richFactor, Description)) +
geom_point(aes(color = p.adjust, size = Count))
2. Visualisez les résultats de l’enrichissement GO. Personnalisez la visualisation en modifiant la résolution de l’arc (n=500), les mappages de couleurs et l’agencement des facettes selon les spécifications expérimentales.
   ego <- enrichGO(gene = gene$ENTREZID, OrgDb = "org.Hs.eg.db", ont = "ALL",
pvalueCutoff = 0.05, qvalueCutoff = 0.05, pAdjustMethod = "BH")
ggplot(df) +
ggforce::geom_link(aes(x = 0, y = Description, xend = -log10(p.adjust),
yend = Description, color = ONTOLOGY), n = 500, show.legend = FALSE) +
facet_wrap(~ONTOLOGY, scales = "free", ncol = 1)
3. Visualisez les résultats de l’analyse d’enrichissement des ensembles de gènes (GSEA). Ajustez le nombre de voies, l’échelle des crêtes de densité et la transparence des couches pour s’adapter aux différentes distributions de taille d’effet.
   genelist <- sort(res$log2FoldChange, decreasing = TRUE)
names(genelist) <- rownames(res)
hallmarks <- read.gmt('resource/h.all.v2023.2.Hs.symbols.gmt')
y <- GSEA(genelist, TERM2GENE = hallmarks, pvalueCutoff = 0.05)
gsearesult <- yd %>% arrange(desc(NES)) %>% slice_head(n = 10)
ggplot(gsearesult, aes(x = logFC, y = Description, fill = -log10(pvalue))) +
geom_density_ridges(alpha = 0.8, scale = 0.8) +
geom_point(aes(size = abs(NES), x = -0.4, color = NES)) +
scale_fill_distiller(palette = 'Spectral') +
scale_color_distiller(palette = 'Reds') +
scale_size_continuous(range = c(2, 6))

Results

Le flux de travail d’analyse en amont pour le seq d’ARN en vrac est illustré à la Figure 1A. Ce flux de travail exécute séquentiellement les étapes clés suivantes sur une plateforme Linux : premièrement, un contrôle de qualité rigoureux des données brutes de séquençage est effectué à l’aide de fastp pour supprimer les lectures de faible qualité et les séquences d’adaptateurs ; par la suite, HISAT2 aligne les lectures de haute qualité au génome de référence, Samtools convertissant et triant l...

Discussion

L’analyse massive de données ARN-seq est caractérisée comme une tâche interdisciplinaire qui intègre la génomique, la bioinformatique, les statistiques et l’informatique. Un flux de travail analytique complet comprend plusieurs étapes en amont et en aval, incluant le prétraitement des données brutes, le contrôle qualité, l’alignement des séquences, la quantification au niveau des gènes, la normalisation des données, l’analyse d’expression différentielle et l’interprétation biologique. Pa...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
biomaRt	Bioconducteur	2.64.0	Annotation génétique d’Ensembl
clusterProfiler	Bioconducteur	4.16.0	Analyse de l’enrichissement fonctionnel
DESeq2	Bioconducteur	1.48.1	Analyse différentielle de l’expression
FactoMineR	AgroParisTech	2.11.0	ACP et analyse multivariée
fastp	OpenGene	1.0.1	Contrôle qualité et filtrage des données FASTQ
Nombre de fonctionnalités	Division de bioinformatique, Institut Walter et Eliza Hall de recherche médicale	2.0.0	; Compter le nombre de lectures mappées à chaque gène pour la quantification de l’expression génique
ggplot2	Postul	3.5.2	Visualisation des données
ggrepel	Kamil Slowikowski	0.9.6	Étiquettes textuelles non chevauchantes
ggridges	Claus O. Wilke	0.5.6	Créer des graphiques de crête
HISAT2	Université Johns Hopkins	2.2.1	Aligner les lectures filtrées de haute qualité au génome de référence
R	Équipe principale R & nbsp ;	4.5.0	Un environnement pour le calcul, l’analyse et la visualisation des données
RColorBrewer	Erich Neuwirth	1.1.3	Palettes de couleurs pour le tracé
samtools	Flux de travail sur la génomique à grande échelle	1.22.0	Convertir et traiter les fichiers SAM pour une récupération et un accès efficaces
Boîte à outils SRA	Centre national d’information biotechnologique	3.2.1	Obtenir et prétraiter les données brutes de séquençage à partir de la base de données NCBI SRA