Análisis transcriptómico basado en datos de ARN-seq a granel

Summary

El protocolo actual establece una cadena completa para analizar el proceso de RNA-seq a granel desde los datos brutos hasta el análisis de enriquecimiento funcional.

Abstract

El hígado graso no alcohólico (NAFL) suele considerarse una condición benigna; sin embargo, una vez que progresa a esteatohepatitis no alcohólica (NASH), los pacientes enfrentan un riesgo significativamente mayor de desarrollar enfermedad hepática en fase terminal. Muchos estudios intentan esclarecer el mecanismo molecular que subyace a la transición de NAFL a NASH. Las tecnologías de secuenciación de alto rendimiento (como el ARN-seq a granel) han proporcionado a los investigadores una comprensión más profunda al examinar el transcriptoma, revelar la expresión de moléculas, la activación de vías de señalización y otros factores asociados a la progresión de la enfermedad. Existe una gran cantidad de datos de código abierto disponibles para que los investigadores los analicen con el fin de identificar posibles objetivos para el tratamiento de enfermedades. Sin embargo, la investigación relacionada está limitada por la falta de un proceso eficiente y fiable para el análisis aguas arriba del transcriptoma. Aquí, se proporciona una línea de análisis upstream altamente reproducible y fácil de usar y posterior análisis diferencial de genes relacionados para lograr un procesamiento estandarizado y un análisis profundo de datos privados o públicos. La cadena se divide en cuatro pasos: (1) control de calidad de los datos; (2) mapeo génico; (3) análisis diferencial de genes; y (4) análisis funcional. Este proceso tiene como objetivo descubrir los mecanismos moleculares de transformación de enfermedades y ayudar a los investigadores a detectar posibles objetivos y enfoques terapéuticos mediante el análisis de datos de ARN-seq a granel.

Introduction

La enfermedad hepática grasa no alcohólica (NAFLD) es la enfermedad hepática crónica más prevalente a nivel mundial, afectando a más de una cuarta parte de la población. Su incidencia ha aumentado drásticamente en las últimas^{décadas 1,2,3}. La creciente carga de enfermedades, especialmente su forma más avanzada, la esteatohepatitis no alcohólica (NASH), supone un gran desafío para la salud global y una gran carga^{económica 4}. La primera etapa de la NAFLD es el hígado graso no alcohólico (NAFL), que va acompañado de inflamación y fibrosis que pueden progresar a NASH. Este último incrementa significativamente el riesgo de progresión hacia enfermedad hepática en fase terminal, incluyendo cirrosis y carcinoma hepatocelular (HCC^)5,6,7. La incidencia y mortalidad por HCC están asociadas con un aumento de^{NASH 8,9}, y se espera que NAFLD/NASH se convierta en la principal indicación para el trasplante hepático para 2030-10. Sin embargo, la progresión clínica de la NAFLD es altamente^{heterogénea 11}, lo que dificulta gravemente el desarrollo de fármacos^relevantes, por lo que es especialmente importante explorar con precisión los mecanismos moleculares implicados.

La adquisición masiva de información composicional celular basada en ARN-seq puede elucidar significativamente la patogénesis de diversas enfermedades. En las últimas décadas, se han realizado numerosos estudios masivos de RNA-seq en organismos modelo y humanos para dilucidar diferencias de expresión génica en la progresión de^NASH 13,14,15, con el fin de identificar nuevos objetivos terapéuticos para la intervención. Basándose en el análisis masivo de RNA-seq, Xiong et al. encontraron que las células no parenquimales (NPCs) en el hígado están implicadas en procesos como la formación de matrices extracelulares y la adhesión celular, que contribuyen a la progresión de NASH¹⁶. Li et al. demostraron que la proteína asociadora del tumor de Wilms hepático (WTAP) en hepatocitos regula la acumulación de lípidos ectópicos y la inflamación, promoviendo así la formación^{de NASH 17}. Aunque el análisis masivo de RNA-seq es una herramienta poderosa para esclarecer los mecanismos de NASH, sus resultados son muy sensibles a la calidad de los datos ascendentes. La heterogeneidad de las operaciones experimentales y procesos de análisis aguas arriba puede perjudicar gravemente la fiabilidad de los datos, enmascarando así la verdadera información biológica e interferiendo con la precisión de los análisis posteriores. Por lo tanto, es importante establecer un conjunto de procedimientos estandarizados de análisis upstream.

En comparación con la secuenciación de ARN unicelular (scRNA-seq), el ARN-seq en masa ofrece varias ventajas tanto en el diseño experimental como en aplicaciones prácticas. Aunque el scRNA-seq permite identificar la heterogeneidad celular a nivel de célula única y permite un análisis preciso de características transcripcionales específicas de cada tipo celular, se asocia a altos costes, requisitos complejos de procesamiento de datos y una sensibilidad limitada para detectar transcritos de baja^{abundancia 18}. En cambio, el RNA-seq a granel proporciona mayor profundidad de secuenciación, menor coste y mayor rendimiento de muestras, lo que lo hace especialmente adecuado para análisis de expresión génica diferencial a nivel poblacional y la exploración de mecanismos^{moleculares 19}. Por lo tanto, guiado por flujos de trabajo analíticos estandarizados, el ARN-seq a granel sigue siendo un enfoque eficiente, rentable y robusto para investigar la base molecular de enfermedades complejas.

Este protocolo está diseñado específicamente para conjuntos de datos de ARN-seq a granel derivados de tejidos humanos con alta integridad de ARN (RIN ≥ 7,0) y suficiente ARN de entrada (≥ 500 ng por muestra). Para garantizar la ejecución fiable de los pasos de alineación y cuantificación, se recomienda una estación de trabajo local equipada con al menos una CPU de 10 núcleos, 32 GB de RAM y un mínimo de 200 GB de espacio libre en disco. Basándose en estos requisitos, el protocolo proporciona un flujo de trabajo analítico eficiente y fácil de usar, incluyendo instrucciones operativas detalladas y configuraciones de parámetros estandarizadas, para satisfacer las necesidades de los investigadores que analizan datos transcriptómicos a gran escala.

Protocol

Para fines demostrativos, se utilizó el conjunto de datos público PRJNA1023502 generado por Lan Bai et al. para ilustrar cada paso tanto de análisis upstream como^{downstream 20}. Como este conjunto de datos se origina en la base de datos de acceso abierto NCBI SRA, no se requieren permisos adicionales ni aprobaciones éticas. Consulte la Tabla de Materiales para verificar todas las versiones requeridas de software y R-package. El conjunto de datos disponible públicamente PRJNA1023502 comprende 6 muestras de RNA-seq de hígado que no son NASH, 6 NAFL y 6 muestras de RNA de hígado NASH. En este protocolo, el conjunto de datos se utilizó para demostrar todos los pasos del flujo de trabajo de RNA-seq en masa, incluyendo la recuperación de datos de la base de datos SRA, control de calidad (fastp), alineación (HISAT2), cuantificación (featureCounts) y análisis posteriores de expresión diferencial y enriquecimiento funcional.

1. Instalación del kit de herramientas SRA

1. Visita la página web oficial de SRA Toolkit y descarga la versión 3.2.1.

2. Descarga de datos públicos

1. Obtén el número SRA.
   1. Accede a él en el artículo^{suplemento 20}, en la sección de Disponibilidad de Datos, o buscando palabras clave en la base de datos NCBI SRA.
   2. Escribe prefetch en el terminal para descargarlo.

3. Generación de la matriz de recuento génico

1. Define los caminos y configuraciones como se describe a continuación:
   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. Descarga el genoma de referencia humano (ensamblaje primario GRCh38) y el archivo de anotación genética correspondiente (gencode.v44, conjunto completo de anotaciones génicas) de la base de datos oficial 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. Renombra los archivos SRA.
   1. Asegúrese de que todos los archivos SRA usen la extensión .sra para facilitar el reconocimiento y procesamiento precisos por parte de las herramientas posteriores.
   2. Ejecuta los siguientes comandos:
      for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
3. Ejecuta el siguiente comando para convertir archivos SRA al formato FASTQ:
   for f in *.sra; do fasterq-dump "$f" --split-files -O $FASTQ_DIR - e 20; done
4. Construye el índice HISAT2 a partir del genoma de referencia.
   hisat2-build $REFERENCE $INDEX
5. Procesa archivos FASTQ con monitorización visual del progreso.
   for fq in $FASTQ_DIR/*.fastq; do
​sample=$(basename "$fq" .fastq)
   1. Para los datos de secuenciación de extremos emparejados, ejecute el siguiente comando
      for fq1 in $FASTQ_DIR/*_1.fastq; do
sample=$(basename "$fq1" _1.fastq)
​fq2=$FASTQ_DIR/${sample}_2.fastq
   2. Realiza control de calidad y filtrado usando fastp y envía los resultados al directorio $OUT_FASTP. Para datos de secuenciación de extremo único, ejecute el siguiente comando:
      fastp \
-i "${fq}" \
-o $OUT_FASTP/${sample}.clean.fastq \
-h $OUT_FASTP/${sample}.html \
-j $OUT_FASTP/${sample}.json \
-w 20
      Para datos de secuenciación de extremos emparejados, ejecute el siguiente 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. Para cada muestra, inspecciona el informe HTML para control visual de calidad y consulta el informe JSON para métricas estructuradas, incluyendo calidad de lectura, distribución de GC, redundancia de lectura y contaminación de adaptadores. Asegúrese de que todas las muestras cumplan los siguientes umbrales de calidad: Q30 ≥ 85%, tasa de contaminación del adaptador < 5% y sin patrones anormales de GC. Almacena todos los archivos de salida en el directorio $OUT_FASTP.
   4. Realiza alineaciones de secuencias usando HISAT2, genera archivos SAM y los envía al directorio $OUT_HISAT2. Para datos de secuenciación de extremo único, ejecute el siguiente comando:
      hisat2 -p 20 \ -x $INDEX \-U $OUT_FASTP/${sample}.clean.fastq \
-S $OUT_HISAT2/${sample}.sam
      Para datos de secuenciación de extremos emparejados, ejecute el siguiente 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. Realizar la alineación HISAT2 de las lecturas filtradas por calidad con el genoma de referencia humano, utilizando un índice construido a partir del archivo FASTA del ensamblaje primario GRCh38. Proporciona comandos de ejemplo tanto para datos de secuenciación de extremo único como de extremo emparejado.
   6. Convierte SAM a BAM, ordena e indexa.
      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 archivo SAM para cada muestra y conviértelo en un archivo BAM ordenado e indexado usando SAMtools. Para muestras transcriptómicas humanas de alta calidad (número de integridad de ARN, RIN ≥ 8,0), asegúrese de que la tasa global de alineación supere el 85%; para muestras típicas de ARN-seq a granel (RIN ≥ 7,0), considerar tasas de alineación ≥ 70% aceptables.
   8. Realiza cuantificación génica usando featureCount.
      featureCounts -T 20 -p -s 0 \
-a $GTF \
-o $OUT_COUNTS /${sample}.counts.txt \
$OUT_HISAT2/${sample}.sorted.bam
Done
   9. Revisa el archivo de salida delimitado por tabulación (*.counts.txt) y su informe resumen (*.counts.txt.summary) generado por featureCounts para cada muestra. Asegurarse de que la tasa de asignación de lectura cumpla con el umbral típico del ≥70% para la secuenciación de ARN a granel humana; Una tasa notablemente más baja puede indicar una orientación de la hebra desajustada, problemas de anotación o mala calidad de alineación. Utiliza el parámetro -s 0 para este conjunto de datos de RNA-seq no específico para cadenas. Para bibliotecas específicas de la hebra, sustituye -s 0 por -s 1 o -s 2 en el comando.
6. Genera la matriz de recuento génico como se describe a continuación.
   1. Inicializar la matriz de recuentos con identificadores génicos y conteos de la primera muestra.
      cut -f1 $(ls $OUT_COUNTS/*.counts.txt | head -1) > all_counts.txt
   2. Suman iterativamente los recuentos de cada muestra en una sola matriz.
      for f in $OUT_COUNTS/*.counts.txt; do
cut -f7 "$f" | paste all_counts.txt - > tmp && mv tmp
all_counts.txt
​done
   3. Añade una línea de encabezado con IDs de muestra a la matriz de conteo.
      samples=$(ls *.counts.txt | sed 's/.counts.txt//' | paste -sd "\t")
echo -e "Geneid\t$samples" | cat - all_counts.txt > counts_matrix.txt
   4. Extrae las longitudes de los genes del archivo GTF (suma de longitudes de exones por gen).
      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. Procesamiento de matrices de conteo en bruto y anotación genética

1. Inicialización y carga del entorno R correspondiente.
   1. Usa la función install.packages() para instalar tidyverse, ggplot2, ggrepel, RColorBrewer, ggridges, FactoMineR. El paquete tidyverse se emplea para la manipulación y trazado de datos; El paquete ggplot2 se emplea para la visualización; El paquete ggrepel se emplea para etiquetas de texto que no se solapan; el paquete RColorBrewer se emplea para paletas de color; El paquete GGRIDGES se emplea para parcelas de cresta; el paquete FactoMineR se emplea para análisis PCA y multivariante.
   2. Usa la función BiocManager::install() para instalar biomaRt, DESeq2, clusterProfiler. El paquete biomaRt se emplea para la anotación genética de Ensembl; el paquete DESeq2 se emplea para el análisis de expresión diferencial; el paquete clusterProfiler se utiliza para el análisis de enriquecimiento funcional.
   3. Lee la matriz de recuento original generada en el paso 2, usa biomaRt para mapear el ID de Ensembl al nombre del gen HGNC para análisis posteriores, elimina la columna de geneides duplicada y genera la matriz de recuento limpia. Utiliza la matriz de conteo original (counts_matrix.csv) como entrada, con los identificadores de genes de Ensembl como filas y las muestras como columnas. Guarda la matriz de conteo limpia como salida (clean_counts_SRA.csv), con símbolos del gen HGNC como filas y muestras como columnas.
      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. Cuantificación de expresión génica

NOTA: Consulte el Archivo Suplementario 1 para el guion detallado.

1. Ejecute el siguiente comando para calcular fragmentos por kilobase de transcrito por millón de lecturas mapeadas (FPKM) y transcripciones por millón (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. Agrupamiento de muestras y visualización de diferencias

1. Evalúa la calidad de los datos examinando el agrupamiento de muestras a través de PCA. Consulte el Archivo Suplementario 1 para el guion detallado.
2. Para la visualización de PCA, utiliza el código que aparece a continuación.
   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. Análisis de expresión diferencial y visualización de resultados

NOTA: Consulte el Archivo Suplementario 1 para el guion detallado.

1. Ejecuta los siguientes comandos para construir el conjunto de datos DESeq2 y establecer el umbral. Solo se conservaron genes con conteo total > 1 en todas las muestras antes del modelado. Identificar genes diferencialmente expresados (DEGs) utilizando los umbrales de significación del valor p ajustado < 0,05 y |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. Visualiza genes expresados diferencialmente con un gráfico volcánico.
   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. Realizar análisis y visualización de enriquecimiento funcional

NOTA: Consulte el Archivo Suplementario 1 para el guion detallado.

1. Visualiza los resultados del enriquecimiento de 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. Visualiza los resultados del enriquecimiento GO. Personaliza la visualización modificando la resolución del arco (n=500), los mapeos de color y la disposición de facetas según especificaciones experimentales.
   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. Visualiza los resultados del análisis de enriquecimiento de conjuntos génicos (GSEA). Ajusta el número de caminos, la escala de crestas de densidad y la transparencia de las capas para adaptarse a diferentes distribuciones de tamaño de efecto.
   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

El flujo de trabajo de análisis upstream para la secuencia masiva de ARN se ilustra en la Figura 1A. Este flujo de trabajo ejecuta secuencialmente los siguientes pasos clave en una plataforma Linux: primero, se realiza un riguroso control de calidad de los datos de secuenciación en bruto usando fastp para eliminar lecturas y secuencias adaptadoras de baja calidad; posteriormente, HISAT2 alinea las lecturas de alta calidad con el genoma de referencia, con Samtools convirtiendo y ordenando los...

Discussion

El análisis masivo de datos de RNA-seq se caracteriza como una tarea interdisciplinar que integra genómica, bioinformática, estadística e informática. Un flujo de trabajo analítico completo abarca múltiples pasos aguas arriba y posterior, incluyendo preprocesamiento de datos en bruto, control de calidad, alineación de secuencias, cuantificación a nivel génico, normalización de datos, análisis de expresión diferencial e interpretación biológica. Entre estos pasos, convertir con precisión ...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
biomaRt	Bioconductor	2.64.0	Anotación génica de Ensembl
clusterProfiler	Bioconductor	4.16.0	Análisis de enriquecimiento funcional
DESeq2	Bioconductor	1.48.1	Análisis de expresión diferencial
FactoMineR	AgroParisTech	2.11.0	ACP y análisis multivariante
fastp	OpenGene	1.0.1	Control de calidad y filtrado de datos FASTQ
RecuentosCaracterísticas	División de Bioinformática, Instituto Walter y Eliza Hall de Investigación Médica	2.0.0	Cuenta el número de lecturas asignadas a cada gen para la cuantificación de la expresión génica
ggplot2	Postul	3.5.2	Visualización de datos
ggrepel	Kamil Slowikowski	0.9.6	Etiquetas de texto que no se solapan
ggridges	Claus O. Wilke	0.5.6	Crear gráficos de crestas
HISAT2	Universidad Johns Hopkins	2.2.1	Alinea las lecturas filtradas de alta calidad con el genoma de referencia
R	Equipo Principal R	4.5.0	Un entorno para la computación, análisis y visualización de datos
RColorBrewer	Erich Neuwirth	1.1.3	Paletas de colores para trazar gráficos
samtools	Corriente de trabajo de Genómica a Gran Escala	1.22.0	Convertir y procesar archivos SAM para una recuperación y acceso eficientes
Kit de herramientas SRA	Centro Nacional de Información Biotecnológica	3.2.1	Obtener y preprocesar los datos de secuenciación en bruto de la base de datos NCBI SRA