Identificação e validação de uma assinatura de expressão gênica BMI1 no sarcoma de Ewing

Nota

Este relatório foi gerado com auxílio de IA sob direção humana geral.

Conteúdo revisado e validado por Tainá Schons em 19/03/26.

Introdução

Este documento descreve o pipeline bioinformático utilizado para identificar uma assinatura transcriptômica associada à atividade do gene BMI1 no sarcoma de Ewing e avaliar seu valor prognóstico em coortes independentes obtidas do GEO.

O estudo foi estruturado em duas etapas:

1. Descoberta: análise de expressão diferencial no dataset GSE16016, que compara tumores com alta e baixa expressão de BMI1 em dados de array de exon (Affymetrix HuEx-1.0-ST).
2. Validação: cálculo de um score de assinatura BMI1 em 10 coortes públicas adicionais e avaliação da associação com sobrevida livre de eventos (EFS) via Kaplan-Meier e regressão de Cox.

1. Ambiente e configuração

Carregando pacotes exigidos: BiocGenerics

Carregando pacotes exigidos: generics

Anexando pacote: 'generics'

Os seguintes objetos são mascarados por 'package:base':

    as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
    setequal, union

Anexando pacote: 'BiocGenerics'

Os seguintes objetos são mascarados por 'package:stats':

    IQR, mad, sd, var, xtabs

Os seguintes objetos são mascarados por 'package:base':

    anyDuplicated, aperm, append, as.data.frame, basename, cbind,
    colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
    get, grep, grepl, is.unsorted, lapply, Map, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
    rbind, Reduce, rownames, sapply, saveRDS, table, tapply, unique,
    unsplit, which.max, which.min

Carregando pacotes exigidos: oligoClasses

Welcome to oligoClasses version 1.70.0

Carregando pacotes exigidos: Biobase

Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

Carregando pacotes exigidos: Biostrings

Carregando pacotes exigidos: S4Vectors

Carregando pacotes exigidos: stats4

Anexando pacote: 'S4Vectors'

O seguinte objeto é mascarado por 'package:utils':

    findMatches

Os seguintes objetos são mascarados por 'package:base':

    expand.grid, I, unname

Carregando pacotes exigidos: IRanges

Anexando pacote: 'IRanges'

O seguinte objeto é mascarado por 'package:grDevices':

    windows

Carregando pacotes exigidos: XVector

Carregando pacotes exigidos: GenomeInfoDb

Anexando pacote: 'Biostrings'

O seguinte objeto é mascarado por 'package:base':

    strsplit

================================================================================

Welcome to oligo version 1.72.0

================================================================================

Setting options('download.file.method.GEOquery'='auto')

Setting options('GEOquery.inmemory.gpl'=FALSE)

Warning: pacote 'limma' foi compilado no R versão 4.5.1

Anexando pacote: 'limma'

O seguinte objeto é mascarado por 'package:oligo':

    backgroundCorrect

O seguinte objeto é mascarado por 'package:BiocGenerics':

    plotMA

Carregando pacotes exigidos: mgcv

Carregando pacotes exigidos: nlme

Anexando pacote: 'nlme'

O seguinte objeto é mascarado por 'package:Biostrings':

    collapse

O seguinte objeto é mascarado por 'package:IRanges':

    collapse

This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.

Carregando pacotes exigidos: genefilter

Carregando pacotes exigidos: BiocParallel

Warning: pacote 'survminer' foi compilado no R versão 4.5.2

Carregando pacotes exigidos: ggplot2

Warning: pacote 'ggplot2' foi compilado no R versão 4.5.2

Carregando pacotes exigidos: ggpubr

Warning: pacote 'ggpubr' foi compilado no R versão 4.5.2

Anexando pacote: 'survminer'

O seguinte objeto é mascarado por 'package:survival':

    myeloma

Warning: pacote 'dplyr' foi compilado no R versão 4.5.2

Anexando pacote: 'dplyr'

O seguinte objeto é mascarado por 'package:nlme':

    collapse

O seguinte objeto é mascarado por 'package:oligo':

    summarize

Os seguintes objetos são mascarados por 'package:Biostrings':

    collapse, intersect, setdiff, setequal, union

O seguinte objeto é mascarado por 'package:GenomeInfoDb':

    intersect

O seguinte objeto é mascarado por 'package:XVector':

    slice

Os seguintes objetos são mascarados por 'package:IRanges':

    collapse, desc, intersect, setdiff, slice, union

Os seguintes objetos são mascarados por 'package:S4Vectors':

    first, intersect, rename, setdiff, setequal, union

O seguinte objeto é mascarado por 'package:Biobase':

    combine

Os seguintes objetos são mascarados por 'package:BiocGenerics':

    combine, intersect, setdiff, setequal, union

O seguinte objeto é mascarado por 'package:generics':

    explain

Os seguintes objetos são mascarados por 'package:stats':

    filter, lag

Os seguintes objetos são mascarados por 'package:base':

    intersect, setdiff, setequal, union

Carregando pacotes exigidos: AnnotationDbi

Anexando pacote: 'AnnotationDbi'

O seguinte objeto é mascarado por 'package:dplyr':

    select

Carregando pacotes exigidos: org.Hs.eg.db

2. Dataset de descoberta — GSE16016

O GSE16016 contém 10 amostras de sarcoma de Ewing geradas na plataforma Affymetrix HuEx-1.0-ST (GPL5175), divididas em grupos de alta e baixa expressão de BMI1 por western blot. Os dados brutos (arquivos .CEL.gz) foram processados com normalização RMA.

2.1 Leitura dos arquivos CEL e normalização RMA

10 arquivos CEL encontrados para GSE16016

Carregando pacotes exigidos: pd.huex.1.0.st.v2

Carregando pacotes exigidos: RSQLite

Carregando pacotes exigidos: DBI

Platform design info loaded.

Reading in : ./aff/GSE16016_RAW/GSM400796.CEL.gz
Reading in : ./aff/GSE16016_RAW/GSM400798.CEL.gz
Reading in : ./aff/GSE16016_RAW/GSM400799.CEL.gz
Reading in : ./aff/GSE16016_RAW/GSM400800.CEL.gz
Reading in : ./aff/GSE16016_RAW/GSM546483.CEL.gz
Reading in : ./aff/GSE16016_RAW/GSM546484.CEL.gz
Reading in : ./aff/GSE16016_RAW/GSM546485.CEL.gz
Reading in : ./aff/GSE16016_RAW/GSM546486.CEL.gz
Reading in : ./aff/GSE16016_RAW/GSM546487.CEL.gz
Reading in : ./aff/GSE16016_RAW/GSM546488.CEL.gz

A normalização usa target = "core", que resume as sondas para o nível de gene (transcript cluster), resultando em um valor por gene por amostra.

Normalizando com RMA (nível 'core')...

Background correcting
Normalizing
Calculating Expression

Normalização concluída: 22011 genes, 10 amostras

Nota técnica: O pacote affy não deve ser carregado junto com oligo nesta etapa, pois ambos exportam funções com o mesmo nome (como probeNames e summarize). Caso affy esteja carregado, utilize detach("package:affy", unload = TRUE) antes de rodar o RMA.

2.2 Definição dos grupos BMI1 High vs Low

Os metadados clínicos do GSE16016 foram obtidos diretamente do GEO. A coluna bmi1 status:ch1 contém a classificação de cada amostra.

Found 1 file(s)

GSE16016_series_matrix.txt.gz

Using locally cached version: ./GSE16016_series_matrix.txt.gz

grupo
High  Low 
   5    5 

2.3 Análise de expressão diferencial com limma

2.4 Anotação e filtragem da assinatura BMI1

Os IDs de transcript cluster foram mapeados para gene symbols usando o pacote de anotação huex10sttranscriptcluster.db. A assinatura foi definida com critérios estritos: adj.P.Val < 0,01 e |logFC| > 1.

'select()' returned 1:many mapping between keys and columns

Assinatura BMI1 identificada com 1219 genes

A assinatura resultante inclui genes diferencialmente expressos entre tumores com alta e baixa atividade de BMI1. Genes com logFC positivo são considerados up-regulados pelo BMI1 e serão utilizados para construir o score de atividade nas coortes de validação.

3. Mega-coorte de validação

Dez datasets públicos do GEO foram selecionados para validação. A função processar_validacao() automaticamente detecta a plataforma de cada estudo e aplica o pacote de anotação correspondente.

Plataforma	GSEs	Pacote
GPL570 (HG-U133 Plus 2)	GSE17674, GSE17618, GSE17679, GSE12102	hgu133plus2.db
GPL5175 (HuEx-1.0-ST)	GSE63157, GSE68776, GSE63155, GSE63156	huex10sttranscriptcluster.db
GPL6244 (HuGene-1.0-ST)	GSE45544	hugene10sttranscriptcluster.db
GPL96 (HG-U133A)	GSE39262	hgu133a.db

3.1 Função de processamento por coorte

3.2 Execução do loop e construção da mega-matriz

--- Processando: GSE17674 ---

Found 1 file(s)

GSE17674_series_matrix.txt.gz

Using locally cached version: ./GSE17674_series_matrix.txt.gz

  Plataforma: GPL570

'select()' returned 1:many mapping between keys and columns

--- Processando: GSE17618 ---

Found 1 file(s)

GSE17618_series_matrix.txt.gz

Using locally cached version: ./GSE17618_series_matrix.txt.gz

  Plataforma: GPL570

'select()' returned 1:many mapping between keys and columns

--- Processando: GSE63157 ---

Found 1 file(s)

GSE63157_series_matrix.txt.gz

Using locally cached version: ./GSE63157_series_matrix.txt.gz

  Plataforma: GPL5175

'select()' returned 1:many mapping between keys and columns

--- Processando: GSE17679 ---

Found 1 file(s)

GSE17679_series_matrix.txt.gz

Using locally cached version: ./GSE17679_series_matrix.txt.gz

  Plataforma: GPL570

'select()' returned 1:many mapping between keys and columns

--- Processando: GSE68776 ---

Found 1 file(s)

GSE68776_series_matrix.txt.gz

Using locally cached version: ./GSE68776_series_matrix.txt.gz

  Plataforma: GPL5175

'select()' returned 1:many mapping between keys and columns

--- Processando: GSE12102 ---

Found 1 file(s)

GSE12102_series_matrix.txt.gz

Using locally cached version: ./GSE12102_series_matrix.txt.gz

  Plataforma: GPL570

'select()' returned 1:many mapping between keys and columns

--- Processando: GSE45544 ---

Found 1 file(s)

GSE45544_series_matrix.txt.gz

Using locally cached version: ./GSE45544_series_matrix.txt.gz

  Plataforma: GPL6244

'select()' returned 1:many mapping between keys and columns

--- Processando: GSE39262 ---

Found 1 file(s)

GSE39262_series_matrix.txt.gz

Using locally cached version: ./GSE39262_series_matrix.txt.gz

  Plataforma: GPL96

'select()' returned 1:many mapping between keys and columns

--- Processando: GSE63155 ---

Found 1 file(s)

GSE63155_series_matrix.txt.gz

Using locally cached version: ./GSE63155_series_matrix.txt.gz

  Plataforma: GPL5175

'select()' returned 1:many mapping between keys and columns

--- Processando: GSE63156 ---

Found 1 file(s)

GSE63156_series_matrix.txt.gz

Using locally cached version: ./GSE63156_series_matrix.txt.gz

  Plataforma: GPL5175

'select()' returned 1:many mapping between keys and columns

Found 1 file(s)

GSE17674_series_matrix.txt.gz

Using locally cached version: ./GSE17674_series_matrix.txt.gz

Found 1 file(s)

GSE17618_series_matrix.txt.gz

Using locally cached version: ./GSE17618_series_matrix.txt.gz

Found 1 file(s)

GSE63157_series_matrix.txt.gz

Using locally cached version: ./GSE63157_series_matrix.txt.gz

Found 1 file(s)

GSE17679_series_matrix.txt.gz

Using locally cached version: ./GSE17679_series_matrix.txt.gz

Found 1 file(s)

GSE68776_series_matrix.txt.gz

Using locally cached version: ./GSE68776_series_matrix.txt.gz

Found 1 file(s)

GSE12102_series_matrix.txt.gz

Using locally cached version: ./GSE12102_series_matrix.txt.gz

Found 1 file(s)

GSE45544_series_matrix.txt.gz

Using locally cached version: ./GSE45544_series_matrix.txt.gz

Found 1 file(s)

GSE39262_series_matrix.txt.gz

Using locally cached version: ./GSE39262_series_matrix.txt.gz

Found 1 file(s)

GSE63155_series_matrix.txt.gz

Using locally cached version: ./GSE63155_series_matrix.txt.gz

Found 1 file(s)

GSE63156_series_matrix.txt.gz

Using locally cached version: ./GSE63156_series_matrix.txt.gz

     GSE17674      GSE17618      GSE63157      GSE17679      GSE68776 
"21357 genes" "21357 genes" "17259 genes" "21357 genes" "17259 genes" 
     GSE12102      GSE45544      GSE39262      GSE63155      GSE63156 
"21357 genes" "19974 genes" "13039 genes" "17259 genes" "17259 genes" 

Antes de construir a mega-matriz, é necessário verificar se algum GSE compartilha amostras com outro. Datasets que pertencem à mesma série ou ao mesmo consórcio frequentemente referenciam os mesmos arquivos CEL em múltiplas submissões, o que causaria duplicação de pacientes na matriz de expressão e distorceria o ComBat e os scores.

GSE17674 x GSE17679 : 62 amostras em comum
GSE17618 x GSE17679 : 55 amostras em comum
GSE63157 x GSE63155 : 46 amostras em comum
GSE63157 x GSE63156 : 39 amostras em comum

A inspeção revelou que GSE17679 é subconjunto de GSE17618 + GSE17674 (mesmos pacientes), e GSE63155 é subconjunto de GSE63157. Esses GSEs foram removidos da mega-matriz para garantir que cada paciente contribua com uma única observação na correção de lote e no cálculo do score.

11839 genes em comum entre todas as coortes

Mega-matriz: 11839 genes x 447 amostras únicas

4. Correção de efeito de lote (ComBat)

Os dados provenientes de plataformas e laboratórios distintos apresentam efeitos de lote sistemáticos que precisam ser removidos antes de qualquer análise combinada. O método ComBat (Johnson et al., 2007) foi aplicado com cada GSE como um lote independente.

Executando ComBat...

Found8batches

Adjusting for0covariate(s) or covariate level(s)

Standardizing Data across genes

Fitting L/S model and finding priors

Finding parametric adjustments

Adjusting the Data

Correção concluída: 11839 genes x 447 amostras

5. Cálculo dos scores da assinatura BMI1

Três versões do score foram testadas para avaliar qual subconjunto de genes da assinatura tem maior valor prognóstico:

• Score UP: média dos genes up-regulados pelo BMI1 (logFC > 0). Alta expressão = alta atividade BMI1.
• Score DOWN: média negada dos genes down-regulados pelo BMI1 (logFC < 0). Como o BMI1 suprime esses genes, a negação alinha a direção do score com a atividade BMI1.
• Score Combinado: média(UP) − média(DOWN), capturando o efeito direcional completo da assinatura.

680 genes UP | 210 genes DOWN disponíveis

Scores calculados para 447 amostras

6. Consolidação dos dados clínicos de sobrevida

Após inspeção manual das colunas de metadados de cada GSE, identificou-se que:

• 4 datasets são estudos de linhagem celular (GSE68776, GSE12102, GSE45544, GSE39262) e não possuem dados de sobrevida de pacientes.
• GSE17679 e GSE63155 foram removidos da mega-matriz por compartilharem IDs de amostras com GSE17618/GSE17674 e GSE63157, respectivamente (seção 3.2).
• GSE63156 partilha os arquivos de expressão com parte do GSE63157, mas possui metadados clínicos independentes — os pacientes clínicos são distintos. Seus scores foram calculados a partir da normalização de lote do GSE63157.

As análises de sobrevida foram realizadas nas 4 coortes com dados clínicos independentes: GSE17618, GSE17674, GSE63156 e GSE63157.

As colunas de tempo e evento variam entre os estudos:

GSE	Coluna de tempo	Unidade	Coluna de evento
GSE17674, GSE17618, GSE17679	efs (months):ch1	meses → convertido para dias (×30,44)	status:ch1 (Dead/AWD = 1)
GSE63157, GSE63155	efs time (days):ch1	dias	efs event:ch1 (No event = 0)
GSE63156	efs-time:ch1	dias	efs (0=no event, 1=event):ch1

Found 1 file(s)

GSE17674_series_matrix.txt.gz

Using locally cached version: ./GSE17674_series_matrix.txt.gz

Found 1 file(s)

GSE17618_series_matrix.txt.gz

Using locally cached version: ./GSE17618_series_matrix.txt.gz

Found 1 file(s)

GSE17679_series_matrix.txt.gz

Using locally cached version: ./GSE17679_series_matrix.txt.gz

Found 1 file(s)

GSE63157_series_matrix.txt.gz

Using locally cached version: ./GSE63157_series_matrix.txt.gz

Found 1 file(s)

GSE63155_series_matrix.txt.gz

Using locally cached version: ./GSE63155_series_matrix.txt.gz

Found 1 file(s)

GSE63156_series_matrix.txt.gz

Using locally cached version: ./GSE63156_series_matrix.txt.gz

# A tibble: 6 × 6
  GSE          n eventos pct_evento tempo_med tempo_max
  <chr>    <int>   <dbl>      <dbl>     <dbl>     <dbl>
1 GSE17618    44      27       61.4       638      5851
2 GSE17674    44      27       61.4       638      5851
3 GSE17679    88      54       61.4       638      5851
4 GSE63155    46      18       39.1      1644      3987
5 GSE63156    39      16       41        2062      4261
6 GSE63157    46      18       39.1      1644      3987

Dataset final: 212 pacientes em 4 coortes independentes

GSE17618 GSE17674 GSE63156 GSE63157 
      44       44       78       46 

7. Análise de sobrevida

Uma função auxiliar é definida para facilitar a comparação dos três scores:

7.1 Comparação dos três scores — Kaplan-Meier por coorte

Score UP (genes ativados pelo BMI1):

Score DOWN (genes suprimidos pelo BMI1):

Score Combinado (UP − DOWN):

Comparação dos p-valores por score:

Coorte	Score UP	Score DOWN	Score Combinado
GSE17618	0,30	0,14	0,058
GSE17674	0,30	0,083	0,14
GSE63156	0,70	0,018	0,35
GSE63157	0,94	0,65	0,94

O score DOWN apresenta o melhor desempenho prognóstico. Em GSE63156, a separação é estatisticamente significativa (p = 0,018), com a direção biologicamente esperada: pacientes com score alto (= baixa expressão dos genes suprimidos pelo BMI1, ou seja, alta atividade de BMI1) apresentam pior EFS. O score combinado melhora ligeiramente GSE17618, mas atenua o sinal em GSE63156.

7.2 Verificação dos hazards proporcionais — Score DOWN no GSE63156

Antes de reportar o HR do score DOWN em GSE63156, a premissa de hazards proporcionais foi verificada com o teste de Schoenfeld.

=== Modelo contínuo ===

       chisq df    p
Score  0.659  1 0.42
GLOBAL 0.659  1 0.42

=== Modelo categórico ===

       chisq df    p
Risk   0.154  1 0.69
GLOBAL 0.154  1 0.69

A premissa de hazards proporcionais não foi violada em nenhum modelo (p = 0,15 para o modelo categórico; p = 0,36 para o contínuo). Os resíduos de Schoenfeld não mostram tendência sistemática ao longo do tempo.

7.3 Forest plot — Cox univariada por coorte (score DOWN)

       GSE  N        HR    CI_low   CI_high           p
1 GSE17618 44 0.3353480 0.1540770 0.7298836 0.005896447
2 GSE17674 44 0.3410660 0.1434029 0.8111833 0.014960103
3 GSE63156 78 1.1215864 0.8616485 1.4599412 0.393660986
4 GSE63157 46 0.6844936 0.3054951 1.5336792 0.357074616
                         label
1 HR=0.34 (0.15–0.73), p=0.006
2 HR=0.34 (0.14–0.81), p=0.015
3 HR=1.12 (0.86–1.46), p=0.394
4 HR=0.68 (0.31–1.53), p=0.357

Warning: `geom_errorbarh()` was deprecated in ggplot2 4.0.0.
ℹ Please use the `orientation` argument of `geom_errorbar()` instead.

`height` was translated to `width`.

Interpretação: Nenhuma coorte apresenta HR significativo no modelo contínuo. A discrepância entre o p-valor do log-rank (p = 0,018 em GSE63156) e o do modelo de Cox contínuo (p > 0,05) reflete a instabilidade das estimativas com apenas 16 eventos: a dicotomização pela mediana cria separação clara de grupos, mas a relação entre o score contínuo e o log-hazard não é bem capturada com tamanho amostral tão reduzido. Os ICs amplos em todas as coortes indicam poder estatístico insuficiente para conclusões definitivas.

8. Análise de sobrevida gene a gene

Para identificar quais genes da assinatura contribuem individualmente para o prognóstico, foi realizada uma análise de Cox univariada por gene, com a expressão como variável contínua, na mega-coorte combinada. A abordagem segue a metodologia descrita por Richter et al. (referência do artigo): Cox univariada por gene + Kaplan-Meier com corte na mediana, com correção para múltiplos testes pelo método de Benjamini-Hochberg (FDR).

8.1 Cox univariada por gene e correção FDR

Coorte combinada: 173 pacientes | 88 eventos | 890 genes a testar

Genes testados: 890

FDR < 0,05: 13

FDR < 0,20: 86

8.2 Verificação dos hazards proporcionais

Antes de reportar qualquer HR, a premissa de hazards proporcionais (PH) foi verificada com o teste de Schoenfeld (cox.zph) para todos os genes com FDR < 0,05.

       Gene  logFC    HR p_cox    FDR p_zph PH_valido
1   ATP5F1E -1.439 1.931 7e-04 0.0452 0.070         ✓
2     ASCC1  1.072 0.548 6e-04 0.0452 0.012 ⚠ VIOLADO
3     FBXO3  1.099 0.499 2e-04 0.0370 0.186         ✓
4    PPP3CB  1.231 0.558 6e-04 0.0452 0.035 ⚠ VIOLADO
5  CALCOCO2  1.264 0.536 5e-04 0.0452 0.001 ⚠ VIOLADO
6     DZIP3  1.300 0.612 2e-04 0.0370 0.261         ✓
7    KLHDC2  1.370 0.546 2e-04 0.0370 0.013 ⚠ VIOLADO
8  TOGARAM1  1.411 0.598 7e-04 0.0452 0.015 ⚠ VIOLADO
9       SRR  1.474 0.548 4e-04 0.0452 0.244         ✓
10     CBR4  1.524 0.598 3e-04 0.0417 0.036 ⚠ VIOLADO
11   NHLRC2  1.687 0.489 0e+00 0.0142 0.025 ⚠ VIOLADO
12     IPO5  1.711 0.660 2e-04 0.0370 0.007 ⚠ VIOLADO
13      RDX  2.256 0.716 7e-04 0.0452 0.800         ✓

Dos 13 genes com FDR < 0,05, apenas 5 têm PH válido e podem ter o HR reportado de forma confiável: ATP5F1E, FBXO3, DZIP3, SRR, RDX. Os demais 8 genes apresentam efeito tempo-variante e devem ser reportados apenas pelo log-rank.

8.3 Heatmap — z de Cox e logFC por gene

Warning: pacote 'pheatmap' foi compilado no R versão 4.5.1

O heatmap está ordenado primeiramente por validade do PH e, dentro de cada bloco, por logFC crescente. A linha separadora divide os genes com HR reportável (verde, bloco superior) dos genes com efeito tempo-variante (laranja, bloco inferior).

8.4 Efeito tempo-variante do NHLRC2

O gene com maior significância (menor FDR) no grupo com PH violado é NHLRC2. Para esse gene, foi ajustado um modelo de Cox com coeficiente tempo-variante da forma β(t) = β₀ + β₁·log(t), usando a função tt().

Call:
coxph(formula = Surv(Time, Event) ~ tt(Expr), data = df_nhlrc2, 
    tt = function(x, t, ...) cbind(x, x * log(t)))

  n= 173, number of events= 88 

             coef exp(coef) se(coef)      z Pr(>|z|)   
tt(Expr)x -3.8302    0.0217   1.4241 -2.690  0.00715 **
tt(Expr)   0.5027    1.6532   0.2282  2.203  0.02760 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

          exp(coef) exp(-coef) lower .95 upper .95
tt(Expr)x    0.0217    46.0724  0.001332    0.3538
tt(Expr)     1.6532     0.6049  1.057002    2.5857

Concordance= 0.659  (se = 0.03 )
Likelihood ratio test= 23.16  on 2 df,   p=9e-06
Wald test            = 22.91  on 2 df,   p=1e-05
Score (logrank) test = 23.93  on 2 df,   p=6e-06

O efeito protetor do NHLRC2 é restrito aos primeiros anos de seguimento: HR ≈ 0,30 aos 180 dias, cruzando HR = 1 em torno de 2037 dias (~5.6 anos). Após esse ponto, o efeito se inverte. Por esse motivo, o HR único do modelo de Cox padrão não deve ser reportado para esse gene.

8.5 Curvas KM — genes com PH válido

Para os 5 genes com premissa de PH satisfeita, as curvas de Kaplan-Meier foram geradas com divisão pela mediana de expressão. Esses são os candidatos mais robustos da assinatura para reportar.

Ignoring unknown labels:
• colour : ""
Ignoring unknown labels:
• colour : ""
Ignoring unknown labels:
• colour : ""
Ignoring unknown labels:
• colour : ""
Ignoring unknown labels:
• colour : ""

[[1]]
NULL

Os 4 genes UP (DZIP3, RDX, FBXO3, SRR) mostram alta expressão associada a melhor EFS (HR = 0,50–0,72), e o gene DOWN (ATP5F1E) apresenta a direção biologicamente esperada: alta expressão (= baixa supressão pelo BMI1) associa-se a pior EFS (HR = 1,93). Todos os cinco genes têm HRs reportáveis com confiança metodológica.

9. Validação transcriptômica em datasets públicos

Esta seção avalia se os 5 genes candidatos (ATP5F1E, FBXO3, DZIP3, SRR, RDX) mostram padrões de co-expressão com BMI1 consistentes com a assinatura identificada no GSE16016, utilizando quatro datasets públicos independentes: GSE68776, GSE12102, GSE45544 e GSE39262.

9.1 BMI1 não está uniformemente elevado em linhagens EWS

Uma premissa implícita da assinatura seria que linhagens de Ewing sarcoma, por terem alta atividade de BMI1, expressariam o gene em níveis superiores a outros tipos celulares. Isso não se confirmou nos dados públicos.

No GSE39262 (linhagens celulares de sarcoma), a expressão de BMI1 em linhagens EWS (mediana = 8,90) é indistinguível da encontrada em outros sarcomas (mediana = 9,02; teste de Wilcoxon, p = 0,78). Da mesma forma, no GSE12102 (tumores primários), BMI1 não difere entre grupos de recaída, NED e metástase (p = 0,39).

Isso é biologicamente esperado: a assinatura foi derivada da variação intra-EWS de BMI1 (alta vs. baixa expressão dentro de tumores de Ewing), não da premissa de que EWS tem BMI1 uniformemente alto em comparação a outros tipos tumorais. A validação relevante, portanto, é verificar se os genes da assinatura co-expressam com BMI1 dentro de coortes independentes.

9.2 Co-expressão com BMI1 em coortes independentes

A correlação de Spearman entre BMI1 e cada gene candidato foi calculada em GSE12102 (N = 37 tumores) e GSE39262 (N = 51 linhagens celulares de sarcoma).

No GSE12102, os quatro genes UP da assinatura (DZIP3, RDX, SRR, FBXO3) apresentam correlação positiva significativa com BMI1 (rho = 0,37–0,58), enquanto ATP5F1E mostra correlação negativa na direção esperada, porém sem significância estatística (rho = −0,25, p = 0,14). No GSE39262, DZIP3 e FBXO3 mantêm correlação positiva significativa (rho = 0,41 e 0,30, respectivamente), mas SRR e RDX não atingem significância nesse dataset — possivelmente pela heterogeneidade dos tipos de sarcoma que confunde a relação.

9.3 Scatter plots: BMI1 × DZIP3 e BMI1 × RDX (GSE12102)

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

As correlações são visualmente consistentes e não parecem ser conduzidas por outliers. A relação positiva entre BMI1 e DZIP3/RDX em tumores primários de Ewing corrobora que esses genes co-variam com a atividade de BMI1 in vivo, fortalecendo sua pertinência como componentes da assinatura.

Nota

Nota sobre terminologia: Os datasets avaliados nesta seção permitem validação associativa (co-expressão transcriptômica), não validação funcional causal. Para estabelecer que BMI1 regula diretamente esses genes, seriam necessários experimentos de ganho e perda de função (knockdown/overexpression) de BMI1 em linhagens EWS — o que está além do escopo dos dados públicos analisados.

10. Meta-análise de efeito randômico — genes candidatos

Para quantificar o efeito prognóstico combinado dos 5 genes com premissa de hazards proporcionais satisfeita, foi realizada uma meta-análise de efeito randômico (método REML) combinando os HRs da Cox univariada nas quatro coortes independentes de validação (GSE17618, GSE17674, GSE63156, GSE63157; N total = 173, 88 eventos).

Warning: pacote 'metafor' foi compilado no R versão 4.5.2

Carregando pacotes exigidos: Matrix

Anexando pacote: 'Matrix'

O seguinte objeto é mascarado por 'package:S4Vectors':

    expand

Carregando pacotes exigidos: metadat

Warning: pacote 'metadat' foi compilado no R versão 4.5.2

Carregando pacotes exigidos: numDeriv

Loading the 'metafor' package (version 4.8-0). For an
introduction to the package please type: help(metafor)

Anexando pacote: 'metafor'

Os seguintes objetos são mascarados por 'package:oligo':

    rma, se

10.1 Resultados poolados

Meta-análise de efeito randômico — Cox univariada por coorte (N=173)

Gene	Direção	HR poolado	IC 95% low	IC 95% high	p-valor	I² (%)	Q heterog. (p)
ATP5F1E	DOWN	3.33	2.00	5.54	0.0000	14.3	0.193
FBXO3	UP	0.25	0.10	0.60	0.0018	55.4	0.079
DZIP3	UP	0.48	0.35	0.65	0.0000	0.0	0.930
SRR	UP	0.44	0.30	0.65	0.0000	0.0	0.953
RDX	UP	0.61	0.49	0.77	0.0000	0.0	0.461

Todos os 5 genes atingem significância estatística no modelo poolado. Os genes DZIP3, SRR e RDX apresentam heterogeneidade nula entre coortes (I² = 0%), indicando efeito altamente consistente. FBXO3 apresenta heterogeneidade moderada (I² = 55%), com o teste Q borderline (p = 0,08), sugerindo variação real do efeito entre as coortes — o modelo de efeito randômico é, portanto, a escolha mais conservadora e adequada. ATP5F1E (gene DOWN, HR = 3,33) mantém baixa heterogeneidade (I² = 14%).

10.2 Forest plots

A direção do efeito é consistente com a assinatura original em todas as coortes: genes UP (DZIP3, RDX, SRR, FBXO3) associam-se a melhor EFS (HR < 1); o gene DOWN (ATP5F1E) associa-se a pior EFS (HR > 1). O efeito de FBXO3 no GSE63156 (HR = 0,49) e GSE63157 (HR = 0,59) aponta na mesma direção, mas com ICs amplos dado o menor número de eventos nessas coortes.

10.3 Investigação da heterogeneidade de FBXO3

FBXO3 é o único gene com I² elevado (55%), justificando uma análise diagnóstica. As quatro coortes utilizam duas plataformas distintas: GPL570 (Affymetrix U133 Plus 2.0) nas coortes GSE17618 e GSE17674, e GPL5175 (Affymetrix Exon 1.0 ST) nas coortes GSE63156 e GSE63157. As plataformas diferem na estratégia de medição: enquanto o U133 Plus 2.0 utiliza probes de oligonucleotídeos selecionados por gene, o Exon Array mede expressão exon a exon e agrega para o nível do gene.

Para identificar se a heterogeneidade de FBXO3 é proporcional à compressão do range de expressão entre plataformas ou desproporcional a ela, foi calculada a razão SD(GPL5175)/SD(GPL570) (compressão de variância) e a razão |logHR|(GPL5175)/|logHR|(GPL570) (atenuação do efeito) para cada gene.

DZIP3, SRR e RDX situam-se acima ou sobre a diagonal — a atenuação do seu efeito entre plataformas é proporcional ou menor do que a compressão de variância, indicando comportamento estável. FBXO3, por contraste, apresenta compressão modesta de SD (razão = 0,77) mas atenuação acentuada do efeito (razão = 0,29): o Exon Array retém 77% da variância de FBXO3, mas o HR reduz para apenas 29% do valor observado no U133 Plus 2.0.

A hipótese mais provável é que os probes do GPL570 capturam seletivamente um transcrito ou isoforma de FBXO3 com maior relevância prognóstica, amplificando o sinal em relação à medida agregada do Exon Array. O HR real de FBXO3 provavelmente situa-se entre os extremos das duas plataformas (0,15–0,49), e o estimado pelo modelo de efeito randômico (HR poolado = 0,25 [0,10–0,60]) é a síntese metodologicamente mais conservadora disponível. A direção protetora é consistente em todas as quatro coortes.

Aviso

Cautela na interpretação de FBXO3: O IC do HR poolado é amplo (0,10–0,60) e o I² = 55% indica heterogeneidade real entre estudos. Para reportar FBXO3 como biomarcador, recomendam-se estudos em plataformas adicionais (preferencialmente RNA-seq) para arbitrar entre os dois níveis de efeito observados.

Conclusão

O pipeline identificou uma assinatura de 1219 genes associados à atividade de BMI1 no sarcoma de Ewing, dos quais 680 up-regulados e 210 down-regulados estavam disponíveis na mega-coorte de validação.

Três versões do score foram testadas em quatro coortes independentes (173 pacientes):

• Score UP: sem separação prognóstica em nenhuma coorte (p = 0,30–0,94).
• Score DOWN: melhor desempenho — separação significativa em GSE63156 (p = 0,018, log-rank) e tendência em GSE17674 (p = 0,083). Premissa de hazards proporcionais não violada (p = 0,15). Porém, nenhuma coorte atinge significância no modelo de Cox contínuo, e GSE63157 (maior coorte, N = 85) permanece sem separação.
• Score Combinado: desempenho intermediário, sem vantagem clara. http://127.0.0.1:16957/graphics/plot_zoom_png?width=2174&height=793 O resultado mais relevante é que os genes suprimidos pelo BMI1 carregam mais sinal prognóstico do que os ativados. A consistência da direção (High Score → pior EFS) em três das quatro coortes é biologicamente coerente, porém as amostras são pequenas demais para conclusões definitivas.

A análise gene a gene identificou 13 genes com FDR < 0,05 na Cox univariada (N = 173, 88 eventos). Desses, apenas 5 têm a premissa de hazards proporcionais satisfeita. A meta-análise de efeito randômico nas quatro coortes independentes fornece os HRs definitivos:

Gene	Direção	HR poolado	IC 95%	p	I²
ATP5F1E	DOWN	3,33	2,00–5,54	< 0,001	14%
DZIP3	UP	0,48	0,35–0,65	< 0,001	0%
SRR	UP	0,44	0,30–0,65	< 0,001	0%
RDX	UP	0,61	0,49–0,77	< 0,001	0%
FBXO3	UP	0,25	0,10–0,60	0,002	55%

DZIP3, SRR e RDX apresentam I² = 0%, indicando efeito altamente consistente entre coortes e plataformas — os candidatos mais robustos da assinatura. ATP5F1E tem o maior HR absoluto (3,33) com baixa heterogeneidade (14%), sendo biologicamente relevante como componente da ATP sintase mitocondrial. FBXO3, embora significativo, apresenta heterogeneidade moderada (I² = 55%) explicada por diferença sistemática entre plataformas de array (ver Seção 10.3); seu HR poolado deve ser interpretado com cautela.

Os 8 genes restantes apresentam efeito tempo-variante — em particular NHLRC2, cujo HR cruza o valor 1 em ~2037 dias, tornando o HR único inadequado. Para esses genes, apenas o log-rank pode ser reportado.