$$\rightleftharpoonup{xx}$$
$$\longleftharp{xx}$$,
$$\longrightharp{xx}$$,
Recuperação de dados e pré-processamento
Este estudo foi realizado utilizando dados transcriptômicos e genéticos disponíveis publicamente; nenhum indivíduo humano ou animal esteve diretamente envolvido. Conjuntos de dados transcriptômicos relevantes para a resistência à terapia do câncer de mama HER2+ foram retirados do banco de dados NCBI Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/)11. Dois conjuntos de dados de RNA-seq, GSE231524 e GSE231525, foram selecionados devido ao foco específico na resistência impulsionada por HER3 e inibição de DUSP6 em linhagens celulares de câncer de mama HER2+ (BT474 e MDA-MB-453). Esses conjuntos de dados incluíam fenótipos parentais, tolerantes a medicamentos e resistentes a medicamentos derivados da exposição a Lapatibebe (1 μM) e do knockdown a DUSP6. Matrizes de contagem bruta e arquivos de metadados correspondentes foram acessados usando os pacotes GEOquery (v2.70.0) e Biobase (v2.62.0) no RStudio (v4.3.2)12. Os metadados foram selecionados para definir dois contrastes principais para cada conjunto de dados: GSE231524 compararam o controle (BT474 parental, Dia 0) com amostras tolerantes e resistentes a medicamentos (Dia 9–Mês 9), enquanto GSE231525 compararam o controle (Scrambled siRNA) com knockdown DUSP6 (DUSP6-KD). O controle de qualidade e a normalização dos dados foram realizados usando o framework DESeq2 (v1.42.0), que aplica uma transformação estabilizadora de variância (VST) para reduzir a heteroscedasticidade e garantir a comparabilidade entre amostras. A distribuição dos dados e os padrões de agrupamento foram avaliados visualmente usando ggplot2 (v3.5.0) e pheatmap (v1.0.12) para confirmar a uniformidade dos dados e identificar potenciais valores aberrantes da análise diferencialde expressão 13,14.
Neste estudo, manteve uma distinção clara entre os achados derivados de conjuntos de dados transcribômicos de linhagens celulares e aqueles obtidos de conjuntos clínicos derivados de pacientes. Os dados de linhagens celulares foram usados principalmente para análises exploratórias, incluindo a identificação de genes expressos diferencialmente e a geração de insights mecanicistas preliminares em modelos experimentais controlados. Em contraste, conjuntos de dados derivados dos pacientes foram empregados para validação externa dos padrões de expressão gênica e avaliação da relevância clínica, incluindo avaliação prognóstica. Assim, os resultados de modelos de linhagens celulares e coortes clínicas são interpretados separadamente para evitar generalizações excessivas e garantir o contexto translacional adequado para todos os achados.
Análise diferencial de expressão gênica
Foi realizada uma análise diferencial de expressão para identificar genes que foram significativamente modulados entre as condições de controle e de tratamento. Contagens normalizadas foram processadas usando o modelo de regressão ajustada (ARM) integrado ao DESeq2 para estimar com precisão as variações log₂ de fold e a significância estatística. A fórmula de desenho foi definida como ~condição, representando os grupos controle versus grupos tratados. Genes com valor p ajustado (FDR) < 0,05 e variação absoluta log₂ ≥ 1 foram considerados significativamente expressos de forma diferencial. A redução da mudança de dobra log₂ foi realizada usando o método apeglm para aumentar a robustez na estimativa do tamanho do efeito. Os resultados da análise foram visualizados usando o EnhancedVolcano (v1.22.0)15 e ggplot216, que geraram gráficos vulcânicos e gráficos MA exibindo a relação entre magnitude de expressão e confiança estatística. Estimativas de dispersão também foram avaliadas dentro do DESeq2 para garantir modelagem precisa de variância e normalização consistente entre réplicasbiológicas 17.
Recuperação e identificação de genes diferencialmente expressos relacionados ao estresse oxidativo mitocondrial (MOS-DEGs)
Para investigar a ligação entre metabolismo energético, estresse oxidativo e resistência a medicamentos, foi compilada uma lista abrangente de genes associados ao estresse mitocondrial e oxidativo a partir de múltiplos bancos de dados, incluindo Human MitoCarta3.018 (https://personal.broadinstitute.org/scalvo/MitoCarta3.0/human.mitocarta3.0.html), Gene Ontology (GO:0006979, response to oxidative stress) (http://geneontology.org/), a via de fosforilação oxidativa da Kyoto Encyclopedia of Genes and Genomes (KEGG) (https://www.genome.jp/kegg/), e o Banco de Dados de Genes de Estresse Oxidativo Humano (HOSGDB) (http://hosgdb.com/). Todos os genes recuperados foram padronizados para símbolos genéticos aprovados pela HGNC usando org. Hs.eg.db (v3.18.0) e AnnotationDbi (v1.64.0), enquanto entradas duplicadas, pseudogenes e RNAs não codificantes foram removidos para garantir a precisão das anotações. O painel de genes de estresse oxidativo mitocondrial curado (genes MOS) foi posteriormente usado como conjunto de referência para integração com os genes diferencialmente expressos identificados em ambos os conjuntos de dados transcriptômicos.
A interseção da lista curada de genes MOS com os DEGs obtidos de GSE231524 e GSE231525 foi realizada em R usando funções dplyr (v1.1.3)19 e base R intersect(). Essa abordagem integrativa permitiu a identificação de MOS-DEGs, representando genes funcionalmente ligados ao metabolismo mitocondrial, regulação redox e adaptação ao estresse oxidativo. A sobreposição entre os conjuntos de dados foi visualizada usando o pacote VennDiagram (v1.7.3) em R para ilustrar genes compartilhados e únicos entre os modelosexperimentais 20. A lista refinada de MOS-DEGs foi usada para análises posteriores, fornecendo insights mecanicistas sobre a reprogramação transcricional e metabólica subjacente à resistência à terapia direcionada ao HER2.
Perfilamento de expressões e visualização de MOS-DEGs
O perfilamento de expressão dos MOS-DEGs identificados foi realizado usando os pacotes ComplexHeatmap (v2.18.0)21 e pheatmap (v1.0.12) no RStudio para visualizar padrões globais de expressão em condições parentais, tolerantes a medicamentos e resistentes. Os dados normalizados de contagem foram transformados usando escala de escores z para padronizar a matriz de expressão gênica entre amostras. O agrupamento foi realizado usando distância euclidiana e métodos completos de ligação para detectar padrões de co-expressão e distinguir perfis transcricionais específicos de cada condição. Mapas de calor e gráficos de agrupamento foram gerados usando ggplot2 para garantir uma diferenciação visual clara entre as condições. Essa abordagem de visualização facilitou a identificação de grupos gênicos associados à atividade mitocondrial, modulação do estresse oxidativo e reprogramação metabólica sob estados resistentes a medicamentos.
Enriquecimento funcional e anotação de vias
Para explorar a importância biológica e os mecanismos regulatórios dos MOS-DEGs identificados, foram realizadas análises de enriquecimento da Gene Ontology (GO) e da Kyoto Encyclopedia of Genes and Genomes (KEGG) utilizando o R Studio (versão 4.3.1). As análises foram realizadas no ambiente tidyverse, utilizando múltiplos pacotes Bioconductor para computação e visualização reprodutíveis. A anotação genética e o mapeamento de identificadores eram realizados usando a organização. Hs.eg.db banco de dados (https://bioconductor.org/packages/org.Hs.eg.db/) baseado no genoma de referência do Homo sapiens (GRCh38). A análise de enriquecimento GO foi realizada usando o pacote clusterProfiler (versão 4.8.1; https://bioconductor.org/packages/clusterProfiler/), que categoriza os genes em três ontologias principais — Processo Biológico (BP), Componente Celular (CC) e Função Molecular (MF). A funçãoenrichGO 22 foi usada com parâmetros definidos para o valor p < 0,05 e o valor p ajustado. (FDR) < 0,05, aplicando o método de correção de Benjamini–Hochberg. Visualizações, incluindo gráficos de barras, gráficos de pontos e diagramas de cordas, foram geradas usando enrichplot (https://bioconductor.org/packages/enrichplot/), ggplot213 (https://cran.r-project.org/web/packages/ggplot2/) e GOplot (https://cran.r-project.org/web/packages/GOplot/). Essas ferramentas forneceram uma visão estruturada dos termos GO enriquecidos e suas associações genéticas.
O enriquecimento de caminhos KEGG foi realizado usando a função enrichKEGG() dentro do pacote clusterProfer, referenciando o banco de dados humano KEGG (https://www.genome.jp/kegg/). O pacote KEGGREST (https://bioconductor.org/packages/KEGGREST/) era usado para recuperação e anotação de dados de caminhos. Caminhos com um valor p ajustado. (valor-q) < 0,05 foram considerados significativos. A visualização e o mapeamento de caminhos foram realizados usando pathview (https://bioconductor.org/packages/pathview/), ggplot2 e enrichplot, enquanto igraph e ggraph foram usados para representaçãode rede 23. Todas as análises e visualizações de enriquecimento foram implementadas no R Studio (v4.3.1) usando código reproduzível e fluxos de trabalho padronizados do Bioconductor, garantindo a identificação confiável de categorias funcionais enriquecidas e vias biológicas associadas aos MOS-DEGs.
Validação baseada em ROC de biomarcadores preditivos em câncer de mama
Para validar o poder preditivo clínico dos MOS-DEGs, foi realizada a análise da curva característica operacional do receptor (ROC) usando a ferramenta online ROCplotter (https://www.rocplot.org/)24. ROCplotter é uma plataforma integrada baseada na web que combina dados de expressão gênica com conjuntos de dados clínicamente anotados de resposta ao tratamento de 3.104 pacientes com câncer de mama, incluindo aquelas tratadas com quimioterapia, terapia hormonal ou agentes anti-HER2.
A análise foi realizada utilizando os parâmetros "resposta completa patológica" como variável de desfecho e "qualquer quimioterapia" como categoria de tratamento. Os valores de expressão gênica derivados dos conjuntos de dados microarray do Affymetrix foram automaticamente estratificados em grupos de respondentes e não respondentes com base em anotações clínicas dentro da plataforma.
A curva da característica operacional do receptor (ROC) (AUC), o teste U de Mann–Whitney, a mudança de dobra e o teste qui-quadrado foram aplicados para avaliar a capacidade de cada gene de distinguir respondentes de não respondedores. A área sob a curva (AUC) foi usada como métrica principal para avaliar o desempenho discriminativo. Valores de AUC superiores a 0,55 com valores de p ROC < 0,05 foram considerados significativos, representando desempenho preditivo moderado típico de biomarcadores transcriptômicos, enquanto a correção da taxa de falsa descoberta (FDR) foi aplicada para manter o rigor analítico.
Todos os MOS-DEGs selecionados foram consultados usando seus correspondentes IDs de sonda Affymetrix. O potencial discriminativo de cada gene foi avaliado entre coortes clínicas de câncer de mama, onde os dados de expressão foram estratificados em grupos respondentes e não respondentes. Curvas ROC, boxplots e resultados estatísticos associados eram gerados diretamente pela plataforma ROCplotter e exportados para visualização e comparação a jusante. A análise quantificou o valor preditivo dos reguladores redox e metabólicos envolvidos no estresse oxidativo mitocondrial. Genes que apresentaram significância preditiva consistente entre amostras clínicas foram retidos para inclusão no painel preditivo final.
Análise de expressão diferencial em tecidos tumorais, normais e metastáticos (análise TNMplot)
Os padrões de expressão dos MOS-DEGs superiores foram analisados em tecidos mamários normais, tumorais e metastáticos usando a ferramenta web TNMplot v2 (https://tnmplot.com/analysis/)25. Tanto os conjuntos de dados de RNA-Seq (TCGA + GTEx + MET500) quanto os de chips genéticos foram examinados para garantir validação multiplataforma. O módulo "Análise Múltipla de Genes" foi utilizado com o Carcinoma Invasivo de Mama como o tipo de tecidoselecionado 26. Os valores de expressão foram transformados log₂ e comparados entre os grupos Tumor vs. Normal (TvsN), Metastático vs. Tumor (MvsT) e Metastático vs. Normal (MvsN). O TNMplot calculou automaticamente a mudança de dobra (FC) e os valores p usando o teste Mann–Whitney U para avaliar a significância estatística. As distribuições de expressão foram visualizadas como boxplots e gráficos de densidade gerados diretamente a partir da interface TNMplot, com verde, vermelho e cinza representando os tecidos normais, tumorais e metastáticos, respectivamente. Todos os números foram exportados em alta resolução para integração na seção de resultados. Essa análise de duas plataformas permitiu a identificação robusta e validação dos principais reguladores metabólicos mitocondriais, associados à progressão do câncerde mama 27.
Análise de sobrevivência e prognóstica usando o plotter de Kaplan–Meier
Para avaliar a relevância prognóstica dos MOS-DEGs no câncer de mama, foi realizada uma análise de sobrevivência utilizando a ferramenta online Kaplan–Meier Plotter (https://kmplot.com/analysis/)28. Este banco de dados integra dados de expressão gênica e sobrevivência de mais de 4.900 pacientes com câncer de mama, derivados de múltiplos conjuntos de dados GEO, EGA e TCGA. A análise foi realizada para sobrevivência livre de recidiva (RFS) utilizando IDs individuais de sonda Affymetrix correspondentes aos genes priorizados: 225609_at (GSR), 201761_at (MTHFD2), 201619_at (PRDX3/AOP1) e 201128_s_at (ACLY). Os pacientes foram divididos em grupos de alta e baixa expressão com base no corte mediano de expressão, e as probabilidades de sobrevivência foram estimadas usando o método de Kaplan–Meier. O teste log-rank foi usado para avaliar significância estatística entre curvas de sobrevivência, e razões de risco (HRs) com intervalos de confiança (ICs) de 95% foram calculados automaticamente pela ferramenta. Todas as análises foram realizadas usando o endpoint do RFS, sem restrição baseada no receptor hormonal ou no status HER2 (ER, PR, HER2 = todos). Amostras redundantes foram removidas e as suposições de riscos proporcionais foram verificadas para garantir robustez estatística. Filtros de controle de qualidade excluíam microarrays tendenciosos. Nenhuma seleção manual da sonda ou correção de valor p para múltiplos testes foi aplicada, de acordo com as configurações padrão do KM Plotter. A significância estatística foi definida como p. < 0,05. Gráficos de sobrevivência foram visualizados e baixados em alta resolução para interpretação adicional, comparando os resultados entre expressores altos e baixos de cada candidato ao geneMOS 29,30.
A presente análise foi iniciada utilizando conjuntos de dados que compreendem exclusivamente amostras de câncer de mama HER2+ para a identificação de genes diferencialmente expressos (DEGs) e genes hub. Posteriormente, a análise de sobrevivência foi realizada sem restrições ao status HER2 (ER, PR, HER2 = todos) para avaliar a relevância prognóstica mais ampla e a generalização dos genes identificados. Essa abordagem foi empregada como uma etapa secundária de validação, em vez de redefinir o foco do estudo. Portanto, as implicações prognósticas dos genes hub identificados são interpretadas com cautela, com conclusões primárias permanecendo específicas para o câncer de mama HER2+ .
As sequências de transcrição canônicas de MTHFD2-201 (ENST00000394053.7) e PRDX3-201 (ENST00000298510.4) foram recuperadas do Ensembl Genome Browser (https://www.ensembl.org)31,32. A anotação e classificação de variantes foram realizadas usando o Ensembl Variant Effect Effector (VEP) (https://www.ensembl.org/vep), que forneceu contexto genômico detalhado, alterações de codons e substituições de aminoácidos para cada variante identificada. Apenas variantes de missense (SNPs não sinônimas) foram selecionadas para análise a jusante.
Previsão de patogenicidade e priorização de variantes
As consequências funcionais de cada nsSNP foram avaliadas usando uma combinação de ferramentas computacionais de previsão. O SIFT (https://sift.bii.a-star.edu.sg) foi aplicado para avaliar a conservação de aminoácidos, classificando variantes com escore ≤ 0,05 comodeletérias 33. O PolyPhen-2 (http://genetics.bwh.harvard.edu/pph2) estimou o impacto estrutural e evolutivo das substituições, onde escores ≥ 0,85 indicaram danoprovável 34. O CADD (https://cadd.gs.washington.edu) forneceu uma pontuação composta de deletério integrando múltiplas anotações, com valores ≥ 20 indicando alto potencialpatogênico 35. Métricas complementares do MetaLR36, Mutation Assessor e REVEL foram integradas da interface VEP para maior confiabilidade daprevisão 37. Variantes que atenderam aos limiares MetaLR ≥ 0,70, Mutation Assessor ≥ 3,5 e REVEL ≥ 0,75 foram priorizadas como provavelmente patogênicas.
Previsão de impacto estrutural e mecanicista
Para avaliar como substituições de aminoácidos influenciam a integridade estrutural e a função bioquímica, cada nsSNP de primeiro nível foi analisado adicionalmente usando MutPred2 (http://mutpred.mutdb.org)38 e DynaMut (http://biosig.unimelb.edu.au/dynamut)39. O MutPred2 estimou a probabilidade de interrupção funcional, incluindo alteração da atividade catalítica, ganho ou perda de resíduos de ligação a metais, alterações na acessibilidade ao solvente e modulação alostérica, com escores ≥ 0,80 classificados como altamente patogênicos. O DynaMut calculou a mudança de energia livre de Gibbs (ΔΔG) entre proteínas do tipo selvagem e mutantes, avaliando a direção e a magnitude da alteração de estabilidade, e gerou visualizações de deslocamentos atômicos e rearranjos de ligações de hidrogênio.
Modelagem estrutural secundária e 3D e perfilamento de acessibilidade de solventes
As estruturas cristalinas resolvidas experimentalmente de MTHFD2 e PRDX3 foram obtidas do Banco de Dados de Proteínas (PDB) e processadas usando PyMOL V:3.1 (https://pymol.org)40 para visualizar a distribuição espacial de resíduos deletérios. Modelos mutantes foram criados introduzindo as substituições correspondentes de aminoácidos, seguidas por refinamento estrutural e minimização de energia. A inspeção comparativa 3D destacou mudanças nos elementos secundários, alteração dos contatos interatômicos e a proximidade espacial dos NSSNPs a domínios catalíticos e de ligação a cofatores, revelando possíveis interrupções nas funções redox e metabólicas.
Análises secundárias de estrutura e exposição a solventes foram realizadas usando PSIPRED V: 3.2 (http://bioinf.cs.ucl.ac.uk/psipred)41,42 e NetSurfP 3.0 (https://services.healthtech.dtu.dk/service.php?NetSurfP-2.0)43. Essas ferramentas previam α-hélices, β-fias, bobinas e regiões desordenadas, além de escores de acessibilidade relativa aos solventes (RSA). Resíduos que apresentavam valores de RSA moderados a altos e ordem estrutural foram mapeados para identificar posições expostas a solventes e funcionalmente críticas. Os sítios afetados foram visualizados em diagramas topológicos 2D para determinar se mutações deletérias ocorreram em núcleos catalíticos rígidos ou regiões de laços flexíveis, prevendo assim seus prováveis efeitos na dinâmica do dobramento de proteínas e na eficiência enzimática.
Validação cruzada de banco de dados, integração funcional e validação de estabilidade
Cada nsSNP priorizado foi cruzado com bancos de dados genômicos em nível populacional, incluindo dbSNP, 1000 Genomes, ExAC e gnomAD para confirmar frequência de variantes, distribuição global de alélicos e associações clínicas previamente relatadas. A integração da conservação evolutiva, modelagem estrutural e predição funcional baseada em aprendizado de máquina permitiu a identificação de variantes deletérias de alta confiança em MTHFD2 e PRDX3. Essas mutações de alto impacto foram posteriormente mapeadas para domínios funcionais para elucidar seu papel potencial no desequilíbrio do estresse oxidativo mitocondrial, sinalização metabólica alterada e resistência terapêutica no câncer de mama. Para confirmar ainda mais as consequências termodinâmicas de cada substituição deletéria, o iMutant 3.0 (https://folding.biofold.org/i-mutant/i-mutant3.0.html)44 foi usado para prever os efeitos das mutações na estabilidade das proteínas usando dados de sequência e estruturas. A análise calculou valores de ΔΔG (kcal/mol) representando a mudança de energia livre entre proteínas do tipo selvagem e mutantes. Variantes que apresentavam valores negativos de ΔΔG foram classificadas como mutações desestabilizadoras, indicando redução da estabilidade proteica e aumento da probabilidade de desdobramento. A integração das previsões do iMutant com os resultados do DynaMut e MutPred2 proporcionou validação cruzada para identificar resíduos estruturalmente críticos que provavelmente afetariam a função redox, a integridade catalítica e a estabilidade conformacional geral da proteína.
Interpretação funcional integrada e relevância terapêutica
Todos os NSSNPs prejudiciais identificados foram validados por meio de cruzamento com bancos de dados populacionais dbSNP, gnomAD e ExAC para verificar frequências alelares menores e associações previamente relatadas com fenótipos do câncer. Interpretações integrativas de dados de conservação evolutiva, modelagem estrutural e estabilidade indicaram que as mutações de alto impacto rs1471336772 (MTHFD2) e rs747786383 (PRDX3) exercem os efeitos prejudiciais mais fortes na conformação proteica e na eficiência catalítica. Os resultados computacionais sugerem coletivamente que mutações em MTHFD2 desestabilizam o metabolismo redox dependente de NADPH, enquanto mutações em PRDX3 prejudicam a defesa contra o estresse oxidativo mediada pela peroxidase, contribuindo para a disfunção mitocondrial e a agressividade tumoral. Essa análise estrutural e funcional baseada em nsSNP fornece uma base computacional para futuras triagens terapêuticas e validação mutacional, destacando MTHFD2 e PRDX3 como biomarcadores de precisão para a terapia do câncer de mama direcionada por redox. Para melhorar a clareza e fornecer uma visão abrangente da estratégia analítica, um fluxo de trabalho esquemático resumindo as principais etapas do estudo é apresentado na Figura 2. O fluxo de trabalho integra análise diferencial de expressão gênica, filtragem gênica mitocondrial, construção de redes de interação proteína-proteína, validação clínica por meio de análise ROC e caracterização estrutural baseada em nsSNP. Essa estrutura passo a passo destaca a progressão lógica do processamento de dados transcriptômicos até a identificação de biomarcadores e interpretação funcional.
Figura 2. Fluxo de trabalho integrativo em múltiplas etapas para identificação e validação de biomarcadores relacionados ao estresse oxidativo mitocondrial no câncer de mama HER2+ . Este esquema resume o pipeline analítico empregado no estudo. Primeiro, foi realizada uma análise diferencial de expressão gênica (DEG) em conjuntos de dados de RNA-seq (GSE231524 e GSE231525) para identificar genes significativamente alterados. Esses DEGs foram intersectados com genes curados relacionados ao estresse oxidativo mitocondrial para obter MOS-DEGs. Em seguida, foi realizada uma análise de interação proteína-proteína (PPI) usando STRING e Cytoscape para identificar genes centrais e módulos funcionais. Posteriormente, a análise da curva de característica operacional do receptor (ROC) foi aplicada usando a plataforma ROCplotter para avaliar o desempenho preditivo de genes selecionados em coortes clínicas. Por fim, foram realizadas análises de SNPs não sinônimos (nsSNP) e modelagem estrutural para avaliar o impacto funcional e estrutural potencial de variantes-chave em genes priorizados (MTHFD2 e PRDX3). Esse fluxo de trabalho integrativo conecta análises transcriptômicas, de rede, clínicas e estruturais para identificar potenciais biomarcadores e alvos terapêuticos. Por favor, clique aqui para ver uma versão ampliada desta figura.
