O estudo foi conduzido de acordo com a Declaração de Helsinque, e o protocolo foi aprovado pelo Comitê de Ética do Hospital Torácico de Anhui (K2025-007) em 22 de abril de 2025. O consentimento informado foi obtido de todos os participantes envolvidos no estudo.
Extração e normalização de dados
Os perfis transcriptômicos e os conjuntos de dados clínicos correspondentes para LUAD foram obtidos dos cohorts do TCGA e do GEO. O conjunto de dados TCGA-LUAD foi designado como conjunto de treinamento, com GSE72094, GSE31210 e GSE26939 servindo como cohorts para validação externa (Tabela 1). Além disso, 900 MCRGs foram coletados de um estudo anterior12 (Tabela Suplementar 1). Os dados do transcriptoma foram anotados utilizando o GENCODE v36 ou os arquivos de anotação da plataforma GPL correspondentes. Os IDs dos probes foram convertidos em símbolos gênicos, genes duplicados foram mesclados usando a função avereps e apenas genes codificadores de proteínas foram mantidos para gerar matrizes de expressão ao nível gênico. Para o conjunto de treinamento TCGA-LUAD, genes com Fragments per kilobase of exon model per million mapped fragments (FPKM) < 1 em mais de 50% das amostras foram filtrados, e os valores de expressão restantes foram transformados por log2 (log2[FPKM+1]). Para os cohorts de validação do GEO, os dados brutos de expressão foram baixados, os IDs dos probes foram mapeados para símbolos gênicos usando os arquivos de anotação da plataforma respectiva e múltiplos probes correspondentes ao mesmo gene foram agrupados pela média de seus valores de expressão. Esses conjuntos de dados foram transformados por log2 quando necessário. Nenhuma correção de efeito de lote entre plataformas foi aplicada entre TCGA e GEO, pois adotamos uma estratégia de padronização por cohort para garantir comparabilidade relativa. Especificamente, tanto para os cohorts de treinamento quanto para os de validação, os valores de expressão gênica foram centralizados e padronizados (transformação por escore z) utilizando a média e o desvio padrão de cada conjunto de dados individualmente. Os mesmos coeficientes de regressão de Cox derivados do conjunto de treinamento foram então usados para calcular os escores de risco para todos os cohorts. Para manter a aplicabilidade clínica e evitar superajuste a qualquer conjunto de validação, a mediana do escore de risco do cohort de treinamento foi utilizada como um ponto de corte fixo para estratificar os pacientes em grupos de alto e baixo risco em todos os cohorts de validação externa. As informações clínicas, incluindo idade, sexo, estágio patológico, estágio Tumor-Nó-Metástase (TNM), tipo histológico, tempo de sobrevida, status de sobrevida e tipo de tecido, foram extraídas quando disponíveis. O desfecho foi a sobrevida global (OS). Amostras com informações incompletas de sobrevida ou tempo de sobrevida < 30 dias foram excluídas. O tempo de sobrevida foi convertido em anos, e o status de sobrevida foi codificado como 0 para vivo e 1 para morto.
Identificação e análises funcionais de genes candidatos
O pacote Limma identificou genes diferencialmente expressos (DEGs) entre amostras tumorais e normais de LUAD no conjunto de treinamento13. Os seguintes critérios definiram os DEGs: |log2FC| > 0,5 e valor p ajustado < 0,05. Posteriormente, o algoritmo de agrupamento difuso mfuzz no pacote R ClusterGVis foi utilizado para particionar os DEGs em clusters distintos de expressão. A análise da Ontologia Genética – Processo Biológico (GO-BP) foi realizada nos cinco principais genes representativos de cada cluster com base em suas pontuações de pertencimento. Um conjunto de genes compartilhados foi obtido pela interseção dos DEGs com os MCRGs. A análise de enriquecimento funcional utilizando Ontologia Genética/Enciclopédia de Genes e Genomas de Kyoto (GO/KEGG) avaliou a relevância biológica dos genes sobrepostos. As redes de interação proteína-proteína (PPI) foram obtidas a partir da base de dados STRING14. Apenas interações com pontuações de confiança > 0,7 foram mantidas para melhorar a confiabilidade da rede.
Triagem de genes prognósticos
O pacote Survival foi utilizado para realizar uma análise de regressão de Cox univariada a fim de identificar genes provavelmente associados à sobrevida global em LUAD15. Genes com valor de p < 0,05 foram considerados indicadores prognósticos potenciais. A coorte de treinamento TCGA-LUAD incluiu 500 pacientes com dados completos de sobrevida, dos quais 216 (43,2%) apresentaram eventos de morte durante o acompanhamento. A razão entre genes candidatos (n = 108) e eventos (n = 216) foi de aproximadamente 1:2, o que é aceitável para análise de regressão de Cox. Posteriormente, a análise de regressão por Seleção e Encolhimento Absoluto Mínimo (LASSO) e o modelo de Boosting Extremo de Gradiente (XGBoost) selecionaram ainda mais as características. Modelos de riscos proporcionais de Cox foram construídos com family = "cox" por meio da função cv.glmnet do pacote glmnet. O parâmetro de regularização ótimo foi identificado utilizando validação cruzada de 10 dobras, sendo o valor λ.min, que representa o erro mínimo de validação cruzada, selecionado como o valor ótimo de λ. Foram extraídos como características candidatas os genes com coeficientes de regressão não nulos. Para o modelo XGBoost, o tempo de sobrevida e o status de sobrevida foram combinados como variável de desfecho, atribuindo-se valores positivos aos eventos de morte e valores negativos aos casos censurados. Os parâmetros foram definidos como objective = "survival: cox" e eval_metric = "cox-nloglik", com 100 iterações e taxa de aprendizado de 0,1. Após o treinamento do modelo, os escores de importância dos genes foram calculados utilizando os valores de ganho das características. Os 20 principais genes foram mantidos após ordenar os escores de importância em ordem decrescente, a fim de reduzir a dimensionalidade das características e a complexidade do modelo. Os genes comuns entre os resultados do LASSO e do XGBoost foram identificados como genes prognósticos candidatos.
Construção e avaliação de um modelo prognóstico
Um modelo prognóstico foi desenvolvido utilizando análise de regressão de Cox multivariada dos genes candidatos identificados. Os escores de risco foram calculados individualmente da seguinte forma:
.
em que Coefi refere-se ao coeficiente para o gene i, e Expi indica o respectivo valor de expressão gênica. Os indivíduos foram subsequentemente divididos em dois grupos: alto risco e baixo risco, utilizando a pontuação mediana de risco como ponto de corte. Em seguida, foram criadas curvas de característica operacional do receptor (ROC) dependentes do tempo. Para avaliar o potencial de superajuste, foi realizada validação interna por bootstrap com 1.000 iterações de reamostragem para calcular o índice C corrigido para viés e as AUCs dependentes do tempo com intervalos de confiança de 95%. Curvas de calibração foram geradas para avaliar a concordância entre as probabilidades de sobrevida previstas e observadas aos 2, 3 e 5 anos. Além disso, foi realizada análise de curva de decisão (DCA) utilizando o pacote ggDCA no R para avaliar o benefício clínico líquido do modelo nos pontos temporais de 2, 3 e 5 anos, quantificando o valor potencial da pontuação de risco na tomada de decisões clínicas ao longo de diferentes probabilidades de corte. As disparidades de sobrevida entre os grupos estratificados por risco e entre outras categorias clínicas foram comparadas utilizando curvas de sobrevida de Kaplan-Meier (KM) com teste log-rank. Além disso, para elucidar as contribuições individuais dos genes para o desempenho do modelo, foi empregada a análise Shapley Additive exPlanations (SHAP) para interpretação explicativa pós-hoc.
Desenvolvimento e validação externa de nomograma
As relações entre os escores de risco calculados e diversas características clínicas (incluindo sexo, idade e estágio TNM) foram examinadas utilizando testes de soma de postos de Wilcoxon ou Kruskal-Wallis para avaliar a aplicabilidade clínica do modelo. Para avaliar se o escore de risco funcionava como um fator prognóstico independente, variáveis clínicas, juntamente com o escore de risco, foram incorporadas a modelos de regressão de Cox multivariados. Em seguida, foi construído um nomograma prognóstico combinando os fatores de risco clínicos independentes (por exemplo, Estágio) e o escore genético de risco, utilizando o pacote R regplot, para individualizar as previsões de probabilidade de sobrevida. Curvas de calibração foram utilizadas para avaliar a concordância entre a probabilidade de sobrevida prevista pelo nomograma e os desfechos reais de sobrevida. Por fim, a capacidade preditiva final e a generalização do sistema integrado de nomograma foram rigorosamente validadas por meio de curvas ROC dependentes do tempo e análises abrangentes de subgrupos clínicos KM nas coortes.
Análises de infiltração imune e subtipos imunes
O CIBERSORT, utilizando a matriz de assinatura de genes de leucócitos (LM22), foi usado para estimar as proporções relativas de 22 tipos de células imunes e avaliar a infiltração de células imunes em pacientes com LUAD. As relações entre os níveis de expressão gênica prognóstica e a infiltração imunológica foram avaliadas por análise de correlação de Spearman. Os escores de imunidade, estroma, pureza tumoral e ESTIMATE foram obtidos por meio do algoritmo ESTIMATE, e o teste de Wilcoxon avaliou as diferenças entre os grupos de risco. Os pacientes com LUAD foram classificados em seis subtipos imunes utilizando o pacote ImmuneSubtypeClassifier16. O teste de Wilcoxon também foi utilizado para comparar as distribuições dos subtipos imunes entre os grupos de risco.
Análises de pontos de controle imunológico, escore imunofenotípico e ciclo de imunidade contra o câncer
Neste estudo, o teste de soma de postos de Wilcoxon foi utilizado para avaliar 21 genes de pontos de controle imunológico17 entre grupos estratificados por risco, com o objetivo de caracterizar o panorama imunológico do LUAD. A correlação de Spearman associou genes prognósticos candidatos aos genes de pontos de controle imunológico. Para avaliar as diferenças na resposta a inibidores de pontos de controle imunológico (ICIs) em pacientes com LUAD em diferentes níveis de risco, dados de escore imunofenotípico (IPS) para tratamentos anti-PD-1 e anti-CTLA-4 foram obtidos do The Cancer Immunome Atlas (TCIA)18, e a base de dados Tracking Tumor Immunophenotype (TIP)19 foi utilizada para avaliar a atividade do ciclo câncer-imunidade, comparando os escores correspondentes entre os grupos de risco.
Análise de mutação somática e sensibilidade a medicamentos
A ferramenta de mutações do TCGA recuperou perfis de mutações somáticas para os casos do TCGA-LUAD a fim de investigar a variação nos padrões de mutação entre os grupos de risco. O pacote maftools processou e visualizou os dados de mutação. Os níveis de carga mutacional tumoral (TMB) foram determinados para cada espécime e comparados entre as duas categorias de risco. A análise de sensibilidade farmacogenômica foi realizada utilizando o pacote pRRophetic de acordo com a base de dados Genomics of Drug Sensitivity in Cancer (GDSC)20. Valores da concentração inibitória máxima (IC50) para fármacos anticâncer foram previstos para cada paciente com LUAD, e as diferenças entre os grupos de risco foram quantificadas utilizando o teste de soma de postos de Wilcoxon.
Avaliação dos níveis de expressão de genes prognósticos
Cada conjunto de dados serviu para avaliar os níveis de transcritos de genes candidatos selecionados relacionados ao desfecho. Para associar a expressão gênica ao prognóstico do paciente, cortes ideais foram derivados por meio da função surv_cutpoint no pacote R survminer. Com base nesses limiares, os casos de LUAD foram categorizados em subconjuntos de alta e baixa expressão para análise subsequente de sobrevida.
Além disso, cinco pares de tumores de CAUP e tecidos normais adjacentes correspondentes foram obtidos do Hospital Torácico de Anhui e a validação por qPCR foi posteriormente realizada. Cada participante forneceu consentimento informado por escrito. Seis genes prognósticos candidatos (PDGFB, LDHA, ZEB2, FKBP4, DMD e S100B) foram selecionados para validação por qPCR. O RNA foi extraído de amostras de tecido homogeneizadas utilizando reagente de extração de RNA, seguido por extração com clorofórmio e precipitação com isopropanol. O espectrofotômetro permitiu a medição da concentração e pureza do RNA. A validação por qPCR dos seis genes prognósticos candidatos foi realizada utilizando uma mistura mestra de PCR baseada em SYBR Green em um sistema de PCR em tempo real: desnaturação inicial a 95 °C por 30 s, seguida por 40 ciclos de 95 °C por 20 s, 55 °C por 20 s e 72 °C por 20 s. A expressão relativa foi calculada e normalizada à Gliceraldeído-3-Fosfato Desidrogenase (GAPDH) utilizando a técnica de 2-ΔΔCt. Os detalhes de todos os reagentes e instrumentos estão fornecidos na Tabela de Materiais.
Análise estatística
As análises estatísticas foram realizadas utilizando software de computação e visualização estatística. A rede de interação proteína-proteína foi visualizada usando software de análise de redes. Após a avaliação da normalidade, o teste t de Student foi utilizado para variáveis contínuas com distribuição normal, e o teste U de Mann-Whitney para aquelas com distribuição não normal.