Análise Transcriptômica Baseada em Dados de RNA-seq em massa

Summary

O protocolo atual estabelece um pipeline completo para analisar o processo de RNA-seq em massa, desde dados brutos até análise de enriquecimento funcional.

Abstract

Fígado gordududo não alcoólico (NAFL) geralmente é considerado uma condição benigna; no entanto, uma vez que progride para esteatohepatite não alcoólica (NASH), os pacientes enfrentam um risco significativamente aumentado de desenvolver doença hepática em estágio terminal. Muitos estudos estão tentando elucidar o mecanismo molecular subjacente à transição da NAFL para a NASH. Tecnologias de sequenciamento de alto débito (como o RNA-seq em massa) proporcionaram aos pesquisadores um entendimento mais profundo ao examinar o transcriptoma, revelar a expressão de moléculas, ativação de vias de sinalização e outros fatores associados à progressão da doença. Há uma grande quantidade de dados de código aberto disponíveis para pesquisadores analisarem a fim de identificar potenciais alvos para o tratamento de doenças. No entanto, pesquisas relacionadas são limitadas pela falta de um processo eficiente e confiável para a análise a montante do transcriptoma. Aqui, é oferecido um pipeline altamente reprodutível e fácil de usar de análise upstream e subsequente análise diferencial de genes relacionados, para alcançar processamento padronizado e análise profunda de dados privados ou públicos. O pipeline é dividido em quatro etapas: (1) controle de qualidade dos dados; (2) mapeamento gênico; (3) análise diferencial de genes; e (4) análise funcional. Esse processo visa descobrir os mecanismos moleculares da transformação de doenças e auxiliar pesquisadores na triagem de potenciais alvos medicamentosos e abordagens terapêuticas por meio da análise de dados de RNA-seq em massa.

Introduction

A doença hepática gordurosa não alcoólica (NAFLD) é a doença hepática crônica mais prevalente no mundo, afetando mais de um quarto da população. Sua incidência aumentou dramaticamente nas últimas^{décadas 1, 2, 3}. A crescente carga de doenças, especialmente sua forma mais avançada, a esteatohepatite não alcoólica (NASH), representa um grande desafio global à saúde e um pesado ônus^{econômico 4}. O primeiro estágio da NAFLD é o fígado gorduroso não alcoólico (NAFL), que é acompanhado de inflamação e fibrose que podem evoluir para NASH. Este último aumenta significativamente o risco de progressão para doença hepática terminal, incluindo cirrose e carcinoma hepatocelular (HCC^)5,6,7. A incidência e mortalidade por HCC estão associadas a um aumento no^{NASH 8,9}, e espera-se que a NAFLD/NASH se torne a principal indicação para transplante de fígado até 2030^. No entanto, a progressão clínica da NAFLD é altamente^{heterogênea 11}, o que dificulta severamente o desenvolvimento de medicamentos^relevantes, tornando particularmente importante explorar com precisão os mecanismos moleculares envolvidos.

A aquisição em massa de informações composicionais celulares baseada em RNA-seq pode elucidar significativamente a patogênese de várias doenças. Nas últimas décadas, diversos estudos em grande escala de RNA-seq foram realizados em organismos modelo e humanos para elucidar diferenças de expressão gênica na progressão 13,14,15 do NASH, identificando novos alvos terapêuticos para intervenção. Com base na análise de RNA-seq em massa, Xiong et al. descobriram que células não parênquimatosas (NPCs) no fígado estão envolvidas em processos como formação de matriz extracelular e adesão celular, que contribuem para a progressão da NASH¹⁶. Li et al. demonstraram que a proteína associativa do tumor de Wilms hepático (WTAP) nos hepatócitos regula o acúmulo e inflamação de lipídios ectópicos, promovendo assim a formação^{de NASH 17}. Embora a análise de RNA-seq em massa seja uma ferramenta poderosa para elucidar os mecanismos do NASH, seus resultados são altamente sensíveis à qualidade dos dados a montante. A heterogeneidade das operações experimentais a montante e dos processos de análise pode prejudicar seriamente a confiabilidade dos dados, mascarando assim informações biológicas verdadeiras e interferindo na precisão das análises subsequentes. Portanto, é importante estabelecer um conjunto de procedimentos padronizados de análise upstream.

Comparado ao sequenciamento de RNA de célula única (scRNA-seq), o RNA-seq em massa oferece várias vantagens distintas tanto no design experimental quanto nas aplicações práticas. Embora o scRNA-seq permita a identificação da heterogeneidade celular no nível de célula única e a análise precisa das características transcricionais específicas do tipo celular, ele está associado a altos custos, requisitos complexos de processamento de dados e sensibilidade limitada para detectar transcritos de baixa^{abundância 18}. Em contraste, o RNA-seq em grande escala oferece maior profundidade de sequenciamento, menor custo e maior capacidade de amostra, tornando-o particularmente adequado para análises de expressão gênica diferencial em nível populacional e para a exploração de mecanismos^{moleculares 19}. Portanto, quando guiado por fluxos de trabalho analíticos padronizados, o RNA-seq em massa permanece uma abordagem eficiente, econômica e robusta para investigar a base molecular de doenças complexas.

Este protocolo foi projetado especificamente para conjuntos de dados de RNA-seq em grande volume derivados de tecidos humanos com alta integridade de RNA (RIN ≥ 7,0) e RNA de entrada suficiente (≥ 500 ng por amostra). Para garantir a execução confiável das etapas de alinhamento e quantificação, recomenda-se uma estação de trabalho local equipada com pelo menos uma CPU de 10 núcleos, 32 GB de RAM e no mínimo 200 GB de espaço livre em disco. Com base nesses requisitos, o protocolo oferece um fluxo de trabalho analítico eficiente e amigável para o usuário, incluindo instruções operacionais detalhadas e configurações padronizadas de parâmetros, para atender às necessidades de pesquisadores que analisam dados transcriptômicos em grande escala.

Protocol

Para fins de demonstração, o conjunto de dados público PRJNA1023502 gerado por Lan Bai et al. foi usado para ilustrar cada etapa tanto das análises upstream quanto^{downstream 20}. Como este conjunto de dados se origina do banco de dados de acesso aberto NCBI SRA, não são necessárias permissões adicionais ou aprovações éticas. Veja a Tabela de Materiais para verificar todas as versões necessárias de software e R-package. O conjunto de dados público PRJNA1023502 compreende 6 amostras não-NASH, 6 NAFL e 6 amostras de RNA-seq hepática NASH. Neste protocolo, o conjunto de dados foi usado para demonstrar todas as etapas do fluxo de trabalho de RNA-seq em massa, incluindo recuperação de dados do banco de dados SRA, controle de qualidade (fastp), alinhamento (HISAT2), quantificação (featureCounts) e análises de expressão diferencial e enriquecimento funcional a jusante.

1. Instalação do kit de ferramentas SRA

1. Visite o site oficial do SRA Toolkit e baixe a versão 3.2.1.

2. Download de dados públicos

1. Obtenha o número SRA.
   1. Acesse-o no suplemento^{do artigo 20}, na seção de Disponibilidade de Dados, ou buscando palavras-chave no banco de dados NCBI SRA.
   2. Digite o prefetch no terminal para baixá-lo.

3. Geração da matriz de contagem gênica

1. Defina caminhos e configurações conforme descrito abaixo:
   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. Baixe o genoma de referência humano (conjunto primário GRCh38) e o arquivo correspondente de anotação genética (gencode.v44, conjunto abrangente de anotação genética) do banco de dados 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. Renomeie arquivos SRA.
   1. Garanta que todos os arquivos SRA utilizem a extensão .sra para facilitar o reconhecimento e processamento precisos por ferramentas posteriores.
   2. Execute os seguintes comandos:
      for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
3. Execute o seguinte comando para converter arquivos SRA no formato FASTQ:
   for f in *.sra; do fasterq-dump "$f" --split-files -O $FASTQ_DIR - e 20; done
4. Construa o índice HISAT2 a partir do genoma de referência.
   hisat2-build $REFERENCE $INDEX
5. Processe arquivos FASTQ com monitoramento visual do progresso.
   for fq in $FASTQ_DIR/*.fastq; do
​sample=$(basename "$fq" .fastq)
   1. Para dados de sequenciamento de extremidade pareada, execute o seguinte comando
      for fq1 in $FASTQ_DIR/*_1.fastq; do
sample=$(basename "$fq1" _1.fastq)
​fq2=$FASTQ_DIR/${sample}_2.fastq
   2. Realize controle de qualidade e filtragem usando o fastp e envie os resultados para o diretório $OUT_FASTP. Para dados de sequenciamento de extremidade única, execute o seguinte comando:
      fastp \
-i "${fq}" \
-o $OUT_FASTP/${sample}.clean.fastq \
-h $OUT_FASTP/${sample}.html \
-j $OUT_FASTP/${sample}.json \
-w 20
      Para dados de sequenciamento de extremidade pareada, execute o seguinte 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 amostra, inspecione o relatório HTML para controle visual de qualidade e consulte o relatório JSON para métricas estruturadas, incluindo qualidade de leitura, distribuição do GC, redundância de leitura e contaminação do adaptador. Certifique-se de que todas as amostras atendam aos seguintes limites de qualidade: Q30 ≥ 85%, taxa de contaminação do adaptador < 5% e sem padrões anormais de GC. Armazene todos os arquivos de saída no diretório $OUT_FASTP.
   4. Realize o alinhamento de sequências usando HISAT2, gere arquivos SAM e envie-os para o diretório $OUT_HISAT2. Para dados de sequenciamento de extremidade única, execute o seguinte comando:
      hisat2 -p 20 \ -x $INDEX \-U $OUT_FASTP/${sample}.clean.fastq \
-S $OUT_HISAT2/${sample}.sam
      Para dados de sequenciamento de extremidade pareada, execute o seguinte 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. Realize o alinhamento HISAT2 das leituras filtradas por qualidade ao genoma de referência humano, usando um índice construído a partir do arquivo FASTA da montagem primária GRCh38. Forneça comandos de exemplo tanto para dados de sequenciamento de extremidade única quanto de extremidade pareada.
   6. Converta SAM para BAM, ordene e indexe.
      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. Gerar um arquivo SAM para cada amostra e convertê-lo em um arquivo BAM ordenado e indexado usando SAMtools. Para amostras transcriptômicas humanas de alta qualidade (número de integridade de RNA, RIN ≥ 8,0), certifique-se de que a taxa geral de alinhamento ultrapasse 85%; para amostras típicas de RNA-seq em massa (RIN ≥ 7,0), considere taxas de alinhamento ≥ 70% aceitáveis.
   8. Realize a quantificação de genes usando featureCount.
      featureCounts -T 20 -p -s 0 \
-a $GTF \
-o $OUT_COUNTS /${sample}.counts.txt \
$OUT_HISAT2/${sample}.sorted.bam
Done
   9. Verifique o arquivo de saída delimitado por tabulação (*.counts.txt) e seu relatório resumo (*.counts.txt.summary) gerado por featureCounts para cada amostra. Garantir que a taxa de atribuição de leitura atinja o limiar típico de ≥70% para RNA-seq humano em massa; Uma taxa consideravelmente menor pode indicar orientação desalinhada dos fios, problemas de anotação ou baixa qualidade de alinhamento. Use o parâmetro -s 0 para este conjunto de dados de RNA-seq não específico para fias. Para bibliotecas específicas de fios, substitua -s 0 por -s 1 ou -s 2 no comando.
6. Gerar a matriz de contagem genética conforme descrito abaixo.
   1. Inicialize a matriz de contagens com IDs de genes e contagens da primeira amostra.
      cut -f1 $(ls $OUT_COUNTS/*.counts.txt | head -1) > all_counts.txt
   2. Somam iterativamente contagens de cada amostra em uma única matriz.
      for f in $OUT_COUNTS/*.counts.txt; do
cut -f7 "$f" | paste all_counts.txt - > tmp && mv tmp
all_counts.txt
​done
   3. Adicione uma linha de cabeçalho com IDs de amostra à matriz de contagens.
      samples=$(ls *.counts.txt | sed 's/.counts.txt//' | paste -sd "\t")
echo -e "Geneid\t$samples" | cat - all_counts.txt > counts_matrix.txt
   4. Extraia os comprimentos dos genes do arquivo GTF (soma dos comprimentos dos exons por 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. Processamento bruto da matriz de contagem e anotação genética

1. Inicialização do ambiente R e carregamento do pacote R relevante.
   1. Use a função install.packages() para instalar tidyverse, ggplot2, ggrepel, RColorBrewer, ggridges, FactoMineR. O pacote tidyverse é empregado para manipulação e plotagem de dados; O pacote ggplot2 é empregado para visualização; O pacote ggrepel é utilizado para rótulos de texto que não se sobrepõem; o pacote RColorBrewer é utilizado para paletas de cores; O pacote GGRIDGES é utilizado para parcelas de cume; o pacote FactoMineR é empregado para análise PCA e multivariada.
   2. Use a função BiocManager::install() para instalar biomaRt, DESeq2, clusterProfiler. O pacote biomaRt é empregado para anotação gênica do Ensembl; o pacote DESeq2 é empregado para análise diferencial de expressão; o pacote clusterProfiler é usado para análise de enriquecimento funcional.
   3. Leia a matriz de contagem original gerada no passo 2, use o biomaRt para mapear o ID do Ensembl ao nome do gene HGNC para análise subsequente, remova a coluna duplicada do geneid e gere a matriz de contagem limpa. Use a matriz de contagens original (counts_matrix.csv) como entrada, com IDs de genes Ensembl como linhas e amostras como colunas. Salve a matriz de contagens limpa como saída (clean_counts_SRA.csv), com símbolos do gene HGNC como linhas e amostras como colunas.
      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. Quantificação da expressão gênica

NOTA: Consulte o Arquivo Suplementar 1 para o script detalhado.

1. Execute o seguinte comando para calcular fragmentos por kilobase de transcrito por milhão de leituras mapeadas (FPKM) e transcrições por milhão (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. Agrupamento de amostras e visualização de diferenças

1. Avalie a qualidade dos dados examinando o agrupamento de amostras por meio da PCA. Consulte o Arquivo Suplementar 1 para o script detalhado.
2. Para visualização da PCA, use o código abaixo.
   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álise diferencial de expressão e visualização de resultados

NOTA: Consulte o Arquivo Suplementar 1 para o script detalhado.

1. Execute os seguintes comandos para construir o conjunto de dados DESeq2 e definir o limiar. Apenas genes com contagens totais > 1 em todas as amostras foram retidos antes da modelagem. Identifique genes diferencialmente expressos (DEGs) usando os limiares de significância do valor p ajustado < 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. Visualize genes diferencialmente expressos com um gráfico vulcâ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álise e visualização de enriquecimento funcional

NOTA: Consulte o Arquivo Suplementar 1 para o script detalhado.

1. Visualize os resultados do enriquecimento 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. Visualize os resultados do enriquecimento do GO. Personalize a visualização modificando a resolução do arco (n=500), mapeamentos de cores e arranjo de facetas de acordo com especificações experimentais.
   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. Visualize os resultados da análise de enriquecimento de conjuntos gênicos (GSEA). Ajuste o número de caminhos, a escala de dorsal de densidade e a transparência das camadas para acomodar diferentes distribuições de tamanho de efeito.
   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

O fluxo de trabalho de análise upstream para RNA-seq em massa é ilustrado na Figura 1A. Esse fluxo de trabalho executa sequencialmente os seguintes passos-chave em uma plataforma Linux: primeiro, um rigoroso controle de qualidade dos dados brutos de sequenciamento é realizado usando fastp para remover leituras e sequências de adaptadores de baixa qualidade; subsequentemente, o HISAT2 alinha leituras de alta qualidade ao genoma de referência, com Samtools convertendo e ordenando os arquivos d...

Discussion

A análise de dados de RNA-seq em massa é caracterizada como uma tarefa interdisciplinar que integra genômica, bioinformática, estatística e ciência da computação. Um fluxo de trabalho analítico completo abrange múltiplas etapas upstream e downstream, incluindo pré-processamento de dados brutos, controle de qualidade, alinhamento de sequências, quantificação em nível de gene, normalização de dados, análise de expressão diferencial e interpretação biológica. Entre esses passos, converter c...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
biomaRt	Biocondutor	2.64.0	Anotação genética do Ensembl
clusterProfiler	Biocondutor	4.16.0	Análise de enriquecimento funcional
DESeq2	Biocondutor	1.48.1	Análise diferencial de expressão
FactoMineR	AgroParisTech	2.11.0	ACP e análise multivariada
fastp	OpenGene	1.0.1	Controle de qualidade e filtragem dos dados do FASTQ
Contagens de Destaques	Divisão de Bioinformática, Instituto Walter e Eliza Hall de Pesquisa Médica	2.0.0	Conte o número de leituras mapeadas para cada gene para quantificação da expressão gênica
ggplot2	Postular	3.5.2	Visualização de dados
ggrepel	Kamil Slowikowski	0.9.6	Rótulos de texto não sobrepostos
ggridges	Claus O. Wilke	0.5.6	Criar plots de crista
HISAT2	Universidade Johns Hopkins	2.2.1	Alinhe as leituras filtradas de alta qualidade ao genoma de referência
R	Equipe Principal R	4.5.0	Um ambiente para computação, análise e visualização de dados
RColorBrewer	Erich Neuwirth	1.1.3	Paletas de cores para plotar
samtools	Fluxo de trabalho em Genômica em Grande Escala	1.22.0	Converter e processar arquivos SAM para recuperação e acesso eficientes
Kit de Ferramentas da SRA	Centro Nacional de Informação em Biotecnologia	3.2.1	Obtenha e pré-processe dados brutos de sequenciamento a partir do banco de dados NCBI SRA