Analisi Trascrittomica basata su dati di RNA-seq in blocco

Summary

Il protocollo attuale stabilisce una pipeline completa per l'analisi del processo di RNA-seq in massa, dai dati grezzi all'analisi dell'arricchimento funzionale.

Abstract

Il fegato grasso non alcolico (NAFL) è solitamente considerato una condizione benigna; tuttavia, una volta che progredisce verso la steatoepatite non alcolica (NASH), i pazienti affrontano un rischio significativamente aumentato di sviluppare malattie epatiche in stadio terminale. Molti studi stanno cercando di chiarire il meccanismo molecolare alla base della transizione da NAFL a NASH. Le tecnologie di sequenziamento ad alto rendimento (come l'RNA-seq in massa) hanno fornito ai ricercatori una comprensione più profonda esaminando il trascrittoma, rivelando l'espressione delle molecole, l'attivazione delle vie di segnalazione e altri fattori associati alla progressione della malattia. Esiste una grande quantità di dati open source che i ricercatori possono analizzare per identificare potenziali bersagli per il trattamento delle malattie. Tuttavia, la ricerca correlata è limitata dalla mancanza di un processo efficiente e affidabile per l'analisi a monte del trascritoma. Qui, viene fornita una pipeline di analisi upstream altamente riproducibile e intuitiva e successivamente correlata di analisi genica differenziale per ottenere un'elaborazione standardizzata e un analisi approfondita di dati privati o pubblici. La pipeline è suddivisa in quattro fasi: (1) controllo qualità dei dati; (2) mappatura genica; (3) analisi genica differenziale; e (4) analisi funzionale. Questo processo mira a scoprire i meccanismi molecolari della trasformazione delle malattie e ad assistere i ricercatori nello screening di potenziali target farmacologici e approcci terapeutici attraverso l'analisi dei dati Bulk RNA-seq.

Introduction

La fegato grasso non alcolica (NAFLD) è la malattia cronica epatica più diffusa a livello globale, colpendo più di un quarto della popolazione. La sua incidenza è aumentata drasticamente negli ultimi decenni ^1,2,3. Il crescente carico di malattie, in particolare la sua forma più avanzata, la steatoepatite non alcolica (NASH), rappresenta una grande sfida sanitaria globale e un pesante onere^{economico 4}. La prima fase della NAFLD è il fegato grasso non alcolico (NAFL), accompagnato da infiammazione e fibrosi che possono progredire in NASH. Quest'ultimo aumenta significativamente il rischio di progressione verso malattie epatiche in stadio terminale, inclusa cirrosi e carcinoma epatocellulare (HCC^)5,6,7. L'incidenza e la mortalità dell'HCC sono associate a un aumento di^{NASH 8,9}, e si prevede che NAFLD/NASH diventeranno la principale indicazione per il trapianto di fegato entro il 2030-10. Tuttavia, la progressione clinica della NAFLD è altamente^{eterogenea 11}, il che ostacola gravemente lo sviluppo dei farmaci rilevanti^, rendendo particolarmente importante esplorare con precisione i meccanismi molecolari coinvolti.

L'acquisizione di informazioni composizionali cellulari basata su RNA-seq in massa può chiarire in modo significativo la patogenesi di varie malattie. Negli ultimi decenni, numerosi studi su RNA-seq in massa sono stati condotti su organismi modello e umani per chiarire differenze di espressione genica nella progressione^NASH 13,14,15, al fine di identificare nuovi bersagli terapeutici per l'intervento. Basandosi sull'analisi bulk RNA-seq, Xiong et al. hanno scoperto che le cellule non parenchimate (NPC) nel fegato sono coinvolte in processi come la formazione della matrice extracellulare e l'adesione cellulare, che contribuiscono alla progressione di NASH¹⁶. Li et al. hanno dimostrato che la proteina associativa del tumore di Wilms 1 (WTAP) epatica negli epatociti regola l'accumulo e l'infiammazione dei lipidi ectopici, favorendo così la formazione^{di NASH 17}. Sebbene l'analisi di RNA-seq in massa sia uno strumento potente per chiarire i meccanismi della NASH, i suoi risultati sono altamente sensibili alla qualità dei dati a monte. L'eterogeneità delle operazioni sperimentali a monte e dei processi di analisi può compromettere seriamente l'affidabilità dei dati, mascherando così le vere informazioni biologiche e interferendo con l'accuratezza delle analisi successive. Pertanto, è importante stabilire un insieme di procedure standardizzate di analisi a monte.

Rispetto al sequenziamento a RNA a singola cellula (scRNA-seq), il segmento di RNA-seq offre diversi vantaggi distinti sia nella progettazione sperimentale che nelle applicazioni pratiche. Sebbene scRNA-seq consenta l'identificazione dell'eterogeneità cellulare a livello di singola cellula e un'analisi precisa delle caratteristiche trascrizionali specifiche per tipo cellulare, è associato ad elevati costi, requisiti complessi di elaborazione dei dati e a una sensibilità limitata per il rilevamento di trascritti a bassa^{abbondanza 18}. Al contrario, il RNA-seq in massa offre una maggiore profondità di sequenziamento, costi inferiori e una maggiore capacità di campionamento, rendendolo particolarmente adatto per analisi di espressione genica differenziale a livello di popolazione e per l'esplorazione dei meccanismi^{molecolari 19}. Pertanto, guidato da flussi di lavoro analitici standardizzati, il RNA-seq globale rimane un approccio efficiente, conveniente e robusto per indagare la base molecolare delle malattie complesse.

Questo protocollo è progettato specificamente per dataset di RNA-seq in blocco derivati da tessuti umani con alta integrità dell'RNA (RIN ≥ 7.0) e sufficiente RNA in input (≥ 500 ng per campione). Per garantire un'esecuzione affidabile dei passaggi di allineamento e quantificazione, si raccomanda una workstation locale dotata di almeno una CPU da 10 core, 32 GB di RAM e un minimo di 200 GB di spazio libero su disco. Sulla base di questi requisiti, il protocollo offre un flusso di lavoro analitico efficiente e facile da usare, includendo istruzioni operative dettagliate e configurazioni standardizzate dei parametri, per soddisfare le esigenze dei ricercatori che analizzano dati trascritomici su larga scala.

Protocol

A scopo dimostrativo, il dataset pubblicamente disponibile PRJNA1023502 generato da Lan Bai et al. è stato utilizzato per illustrare ogni fase sia delle analisi upstream che^{downstream 20}. Poiché questo dataset proviene dal database NCBI SRA ad accesso aperto, non sono necessarie ulteriori autorizzazioni o approvazioni etiche. Consulta la Tabella dei Materiali per verificare tutte le versioni richieste del software e dei R-package. Il dataset pubblicamente disponibile PRJNA1023502 comprende 6 campioni di RNA epatico non NASH, 6 NAFL e 6 NASH RNA-seq. In questo protocollo, il dataset è stato utilizzato per dimostrare tutte le fasi del flusso di lavoro bulk RNA-seq, inclusi il recupero dati dal database SRA, il controllo qualità (fastp), l'allineamento (HISAT2), la quantificazione (featureCounts) e le analisi differenziali e di arricchimento funzionale a valle.

1. Installazione del kit strumenti SRA

1. Visita il sito ufficiale del SRA Toolkit e scarica la versione 3.2.1.

2. Download pubblico dei dati

1. Ottieni il numero SRA.
   1. Accedi nell'articolo^{supplemento 20}, nella sezione Data Availability, oppure cercando parole chiave nel database NCBI SRA.
   2. Digita prefetch nel terminale per scaricarlo.

3. Generazione della matrice del conteggio genico

1. Definisci percorsi e ambientazioni come descritto di seguito:
   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. Scarica il genoma di riferimento umano (assemblaggio primario GRCh38) e il corrispondente file di annotazione genica (gencode.v44, set completo di annotazione genica) dal database ufficiale 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. Rinomina i file SRA.
   1. Assicurarsi che tutti i file SRA utilizzino l'estensione .sra per facilitare un riconoscimento e un'elaborazione accurati da parte degli strumenti a valle.
   2. Esegui i seguenti comandi:
      for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
3. Esegui il seguente comando per convertire i file SRA in formato FASTQ:
   for f in *.sra; do fasterq-dump "$f" --split-files -O $FASTQ_DIR - e 20; done
4. Costruisci l'indice HISAT2 dal genoma di riferimento.
   hisat2-build $REFERENCE $INDEX
5. Elabora file FASTQ con monitoraggio visivo dei progressi.
   for fq in $FASTQ_DIR/*.fastq; do
​sample=$(basename "$fq" .fastq)
   1. Per i dati di sequenziamento a coppia, esegui il seguente comando
      for fq1 in $FASTQ_DIR/*_1.fastq; do
sample=$(basename "$fq1" _1.fastq)
​fq2=$FASTQ_DIR/${sample}_2.fastq
   2. Esegui il controllo qualità e il filtraggio usando fastp e invia i risultati nella directory $OUT_FASTP. Per i dati di sequenziamento a estremità singola, esegui il seguente comando:
      fastp \
-i "${fq}" \
-o $OUT_FASTP/${sample}.clean.fastq \
-h $OUT_FASTP/${sample}.html \
-j $OUT_FASTP/${sample}.json \
-w 20
      Per i dati di sequenziamento a fine accoppiata, esegui il seguente comando:
      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. Per ogni campione, ispezionare il report HTML per il controllo qualità visivo e consultare il report JSON per metriche strutturate, inclusi la qualità di lettura, la distribuzione GC, la ridondanza di lettura e la contaminazione degli adattatori. Assicurarsi che tutti i campioni rispettino le seguenti soglie di qualità: Q30 ≥ 85%, tasso di contaminazione degli adattatori < 5% e nessun pattern GC anomalo. Memorizza tutti i file di output nella cartella $OUT_FASTP.
   4. Esegui l'allineamento delle sequenze usando HISAT2, genera file SAM e inviali nella directory $OUT_HISAT2. Per i dati di sequenziamento a estremità singola, esegui il seguente comando:
      hisat2 -p 20 \ -x $INDEX \-U $OUT_FASTP/${sample}.clean.fastq \
-S $OUT_HISAT2/${sample}.sam
      Per i dati di sequenziamento a fine accoppiata, esegui il seguente comando:
      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. Eseguire l'allineamento HISAT2 delle letture filtrate per qualità al genoma di riferimento umano, utilizzando un indice costruito a partire dal file FASTA dell'assemblaggio primario GRCh38. Fornire comandi di esempio sia per i dati di sequenziamento single-end che paired-end.
   6. Converti SAM in BAM, ordina e indicizza.
      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. Genera un file SAM per ogni campione e convertilo in un file BAM ordinato e indicizzato usando SAMtools. Per campioni trascrittomici umani di alta qualità (numero di integrità RNA, RIN ≥ 8,0), assicurarsi che il tasso complessivo di allineamento superi l'85%; per i campioni tipici di RNA-seq in massa (RIN ≥ 7,0), considerare i tassi di allineamento ≥ 70% accettabili.
   8. Esegui la quantificazione genica usando featureCount.
      featureCounts -T 20 -p -s 0 \
-a $GTF \
-o $OUT_COUNTS /${sample}.counts.txt \
$OUT_HISAT2/${sample}.sorted.bam
Done
   9. Controlla il file di output delimitato per tabulazioni (*.counts.txt) e il suo rapporto riepilogo (*.counts.txt.summary) generato dai featureCounts per ogni campione. Assicurarsi che il tasso di assegnazione delle letture rispetti la soglia tipica del ≥70% per la sequenza di RNA in massa umana; Un tasso nettamente più basso può indicare un orientamento del filo discorrispondente, problemi di annotazione o scarsa qualità dell'allineamento. Usa il parametro -s 0 per questo dataset RNA-seq non specifico per il filamento. Per le librerie specifiche per filo, sostituisci -s 0 con -s 1 o -s 2 nel comando.
6. Genera la matrice del conteggio genico come descritto di seguito.
   1. Inizializzare la matrice dei conteggi con ID geni e conteggi dal primo campione.
      cut -f1 $(ls $OUT_COUNTS/*.counts.txt | head -1) > all_counts.txt
   2. Sommano iterativamente i conteggi di ogni campione in una singola matrice.
      for f in $OUT_COUNTS/*.counts.txt; do
cut -f7 "$f" | paste all_counts.txt - > tmp && mv tmp
all_counts.txt
​done
   3. Aggiungi una linea di intestazione con gli ID dei campioni alla matrice dei conteggi.
      samples=$(ls *.counts.txt | sed 's/.counts.txt//' | paste -sd "\t")
echo -e "Geneid\t$samples" | cat - all_counts.txt > counts_matrix.txt
   4. Estrarre le lunghezze dei geni dal file GTF (somma delle lunghezze degli esoni per gene).
      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. Elaborazione grezze della matrice di conteggio e annotazione genica

1. Avvializzazione e caricamento dell'ambiente R del pacchetto R pertinente.
   1. Usa la funzione install.packages() per installare tidyverse, ggplot2, ggrepel, RColorBrewer, ggridges, FactoMineR. Il pacchetto tidyverse è impiegato per la manipolazione e la tracciatura dei dati; Il pacchetto ggplot2 viene utilizzato per la visualizzazione; Il pacchetto ggrepel è impiegato per etichette di testo non sovrapposte; il pacchetto RColorBrewer viene utilizzato per le palette di colori; Il pacchetto GGRIDGES viene impiegato per i lotti di cresta; il pacchetto FactoMineR è impiegato per l'analisi PCA e multivariata.
   2. Usa la funzione BiocManager::install() per installare biomaRt, DESeq2, clusterProfiler. Il pacchetto biomaRt viene utilizzato per l'annotazione genica da Ensembl; il pacchetto DESeq2 viene impiegato per l'analisi differenziale dell'espressione; il pacchetto clusterProfiler è utilizzato per l'analisi di arricchimento funzionale.
   3. Leggi la matrice di conteggio originale generata nel passo 2, usa biomaRt per mappare l'ID di Ensembl al nome del gene HGNC per l'analisi successiva, rimuovi la colonna di Geneid duplicata e genera la matrice di conteggio pulita. Usa la matrice di conteggio originale (counts_matrix.csv) come input, con gli ID geni di Ensembl come righe e i campioni come colonne. Salva la matrice dei conteggi pulita come output (clean_counts_SRA.csv), con i simboli del gene HGNC come righe e i campioni come colonne.
      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. Quantificazione dell'espressione genica

NOTA: Consulta il File Supplementare 1 per lo script dettagliato.

1. Esegui il seguente comando per calcolare frammenti per kilobase di trascrittazione per milione di letture mappate (FPKM) e trascrizioni per milione (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. Clustering del campione e visualizzazione delle differenze

1. Valuta la qualità dei dati esaminando il clustering dei campioni tramite PCA. Consulta il File Supplementare 1 per la scrittura dettagliata.
2. Per la visualizzazione PCA, usa il codice qui sotto.
   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. Analisi differenziale dell'espressione e visualizzazione dei risultati

NOTA: Consulta il File Supplementare 1 per lo script dettagliato.

1. Esegui i seguenti comandi per costruire il dataset DESeq2 e impostare la soglia. Solo i geni con conteggi totali > 1 su tutti i campioni sono stati mantenuti prima della modellazione. Identificare geni espressi differenzialmente (DEG) utilizzando le soglie di significatività del valore p corretto < 0,05 e |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. Visualizza geni espressi differenzialmente con un grafico vulcanico.
   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. Eseguire analisi e visualizzazione dell'arricchimento funzionale

NOTA: Consulta il File Supplementare 1 per lo script dettagliato.

1. Visualizza i risultati dell'arricchimento 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. Visualizza i risultati dell'arricchimento GO. Personalizza la visualizzazione modificando la risoluzione dell'arco (n=500), le mappe dei colori e l'organizzazione delle facce secondo le specifiche sperimentali.
   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. Visualizza i risultati dell'analisi di arricchimento del set genico (GSEA). Regolare il numero di percorsi, la scala della cresta di densità e la trasparenza degli strati per adattarsi alle diverse distribuzioni della dimensione degli effetti.
   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

Il flusso di lavoro di analisi upstream per il RNA-seq in massa è illustrato nella Figura 1A. Questo flusso di lavoro esegue sequenzialmente i seguenti passaggi chiave su una piattaforma Linux: innanzitutto, un rigoroso controllo qualità dei dati grezzi di sequenziamento viene eseguito usando fastp per rimuovere letture e sequenze di adattatore di bassa qualità; successivamente, HISAT2 allinea le letture di alta qualità al genoma di riferimento, con Samtools che converte e ordina i file di a...

Discussion

L'analisi dei dati RNA-seq in blocco è caratterizzata come un compito interdisciplinare che integra genomica, bioinformatica, statistica e informatica. Un flusso di lavoro analitico completo comprende molteplici fasi a monte e a valle, inclusi la pre-elaborazione dei dati grezzi, il controllo di qualità, l'allineamento delle sequenze, la quantificazione a livello genico, la normalizzazione dei dati, l'analisi dell'espressione differenziale e l'interpretazione biologica. Tra questi passag...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
biomaRt	Bioconduttore	2.64.0	Annotazione genica da Ensembl
clusterProfiler	Bioconduttore	4.16.0	Analisi dell'arricchimento funzionale
DESeq2	Bioconduttore	1.48.1	Analisi differenziale dell'espressione
FactoMineR	AgroParisTech	2.11.0	PCA e analisi multivariata
fastp	OpenGene	1.0.1	Controllo qualità e filtraggio dei dati FASTQ
Conteggio delle funzionalità	Divisione Bioinformatica, Istituto Walter ed Eliza Hall di Ricerca Medica	2.0.0	Contare il numero di letture mappate su ciascun gene per la quantificazione dell'espressione genica
ggplot2	Ipotesi	3.5.2	Visualizzazione dei dati
ggrepel	Kamil Slowikowski	0.9.6	Etichette di testo non sovrapposte
ggridges	Claus O. Wilke	0.5.6	Crea grafici di cresta
HISAT2	Università Johns Hopkins	2.2.1	Allinea le letture filtrate di alta qualità al genoma di riferimento
R	R Core Team	4.5.0	Un ambiente per il calcolo, l'analisi e la visualizzazione dei dati
RColorBrewer	Erich Neuwirth	1.1.3	Palette di colori per la tracciatura
samtools	Flusso di lavoro sulla Genomica su larga scala	1.22.0	Converti ed elabora file SAM per un recupero e un accesso efficienti
Kit degli Attrezzi SRA	Centro Nazionale per l'Informazione Biotecnologica	3.2.1	Ottenere e pre-processare i dati grezzi di sequenziamento dal database NCBI SRA