Identificação e validação de uma assinatura de expressão gênica BMI1 no sarcoma de Ewing
Data de Publicação
9 de março de 2026
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 09/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:
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).
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
library(oligo)
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
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
O seguinte objeto é mascarado por 'package:dplyr':
select
Carregando pacotes exigidos: org.Hs.eg.db
library(hgu133plus2.db)
library(hugene10sttranscriptcluster.db)
library(hgu133a.db)
setwd("C:/Projetos/Nicole_GEO_DEGS/Nicole")
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
# Busca recursiva para encontrar CEL files do GSE16016all_cels <-list.files(".", pattern ="CEL.gz$",recursive =TRUE, full.names =TRUE,ignore.case =TRUE)cels_16016 <- all_cels[grepl("GSE16016", all_cels)]# Remover duplicatas (mesmo arquivo em pastas diferentes)cels_16016 <- cels_16016[!duplicated(basename(cels_16016))]message(length(cels_16016), " arquivos CEL encontrados para GSE16016")
10 arquivos CEL encontrados para GSE16016
dados_brutos <-read.celfiles(cels_16016)
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.
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.
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.
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.
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.
# Identificar todos os pares de GSEs com sobreposição de IDs de amostrasids_por_gse <-lapply(lista_matrizes, colnames)gse_nomes <-names(ids_por_gse)for (i inseq_along(gse_nomes)) {for (j inseq_along(gse_nomes)) {if (j <= i) next n <-length(intersect(ids_por_gse[[i]], ids_por_gse[[j]]))if (n >0) cat(gse_nomes[i], "x", gse_nomes[j], ":", n, "amostras em comum\n") }}
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.
# Manter apenas coortes com pacientes únicos na matriz de expressãogses_unicos <-c("GSE17674", "GSE17618", "GSE63157", "GSE68776","GSE12102", "GSE45544", "GSE39262", "GSE63156")lista_matrizes_clean <- lista_matrizes[gses_unicos]# Genes presentes em todos os datasetsgenes_comuns <-Reduce(intersect, lapply(lista_matrizes_clean, rownames))message(length(genes_comuns), " genes em comum entre todas as coortes")
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.
message("Correção concluída: ", nrow(mega_matriz_limpa), " genes x ",ncol(mega_matriz_limpa), " amostras")
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.
# Preparar assinatura: um símbolo por gene (maior |logFC| tem prioridade)assinatura_symbols <- assinatura_bmi1[order(abs(assinatura_bmi1$logFC), decreasing =TRUE), ]assinatura_symbols <- assinatura_symbols[!duplicated(assinatura_symbols$Symbol), ]rownames(assinatura_symbols) <- assinatura_symbols$Symbol# Intersecção com os genes disponíveis na mega-matrizgenes_comuns_assinatura <-intersect(rownames(assinatura_symbols),rownames(mega_matriz_limpa))# Score UPgenes_up <- assinatura_symbols$Symbol[assinatura_symbols$logFC >0]genes_up_final <-intersect(genes_up, genes_comuns_assinatura)# Score DOWN (negado para alinhar com direção do BMI1)genes_down <- assinatura_symbols$Symbol[assinatura_symbols$logFC <0]genes_down_final <-intersect(genes_down, genes_comuns_assinatura)message(length(genes_up_final), " genes UP | ", length(genes_down_final), " genes DOWN disponíveis")
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
extrair_clinico_manual <-function(gse_id) { gse_obj <-getGEO(gse_id, destdir =".", getGPL =FALSE)[[1]] meta <-pData(gse_obj) df <-switch(gse_id,# GPL570: tempo em meses, evento = Dead ou AWD"GSE17674"= ,"GSE17618"= ,"GSE17679"=data.frame(ID =rownames(meta),GSE = gse_id,Time =as.numeric(meta[["efs (months):ch1"]]) *30.44,Event =ifelse(grepl("Dead|AWD", meta[["status:ch1"]],ignore.case =TRUE), 1, 0) ),# GPL5175 grupo 1: tempo em dias, evento textual"GSE63157"= ,"GSE63155"=data.frame(ID =rownames(meta),GSE = gse_id,Time =as.numeric(meta[["efs time (days):ch1"]]),Event =ifelse(grepl("No event", meta[["efs event:ch1"]],ignore.case =TRUE), 0, 1) ),# GPL5175 grupo 2: tempo em dias, evento numérico 0/1"GSE63156"=data.frame(ID =rownames(meta),GSE = gse_id,Time =as.numeric(meta[["efs-time:ch1"]]),Event =as.integer(meta[["efs (0=no event, 1=event):ch1"]]) ),NULL# GSE68776, GSE12102, GSE45544, GSE39262: sem dados de pacientes )if (!is.null(df)) df <- df[!is.na(df$Time) & df$Time >0, ]return(df)}dados_clinicos_total <-do.call( rbind,lapply(c("GSE17674", "GSE17618", "GSE17679","GSE63157", "GSE63155", "GSE63156"), extrair_clinico_manual))
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
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.
# Construir dataset com score DOWNfinal_down <-merge( dados_clinicos_total[!duplicated(dados_clinicos_total$ID), ],data.frame(ID =names(zscore_down), Score = zscore_down),by ="ID") |> dplyr::filter(!is.na(Time), Time >0, !is.nan(Score)) |> dplyr::group_by(GSE) |> dplyr::mutate(Risk =ifelse(Score >median(Score), "High", "Low")) |> dplyr::ungroup() |>as.data.frame()df_63156 <- final_down[final_down$GSE =="GSE63156", ]# Teste de Schoenfeld (score contínuo e categórico)cox_cont <-coxph(Surv(Time, Event) ~ Score, data = df_63156)cox_cat <-coxph(Surv(Time, Event) ~ Risk, data = df_63156)cat("=== Modelo contínuo ===\n"); print(cox.zph(cox_cont))
plot(cox.zph(cox_cat),main ="Resíduos de Schoenfeld — Score DOWN categórico (GSE63156)",ylab ="Beta(t)", xlab ="Tempo (dias)")
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)
ggplot(cox_df, aes(x = HR, y = GSE)) +geom_vline(xintercept =1, linetype ="dashed", color ="grey50") +geom_errorbarh(aes(xmin = CI_low, xmax = CI_high), height =0.2) +geom_point(aes(size = N, color = p <0.05), shape =15) +geom_text(aes(label = label), vjust =-1.2, size =3.2) +scale_color_manual(values =c("TRUE"="#E41A1C", "FALSE"="grey40"),labels =c("TRUE"="p < 0,05", "FALSE"="p ≥ 0,05"),name ="" ) +scale_size_continuous(name ="N pacientes", range =c(3, 6)) +scale_x_log10(breaks =c(0.25, 0.5, 1, 2, 4)) +labs(title ="Forest plot — Cox univariada por coorte (Score DOWN)",subtitle ="Variável contínua: BMI1-Signature Z-Score DOWN (EFS)",x ="Hazard Ratio (escala log)",y =NULL ) +theme_bw(base_size =13) +theme(legend.position ="bottom")
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.
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.
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.
Nota metodológica: Uma versão anterior deste pipeline apresentava p-valores menores, porém parcialmente artefatuais. Datasets como GSE17679 e GSE63155 referenciavam os mesmos pacientes de GSE17618/GSE17674 e GSE63157, respectivamente, e ao serem incluídos como lotes separados no ComBat inflavam artificialmente a separação entre grupos.
Próximos passos sugeridos:
Refinar a assinatura DOWN com genes que tenham evidência direta de regulação por BMI1 na literatura (ex.: alvos de RING1B), reduzindo o ruído dos 210 genes incluídos.
Realizar meta-análise de efeito randômico combinando os HRs das quatro coortes para obter uma estimativa consolidada.
Avaliar os datasets de linhagem celular (GSE68776, GSE12102, GSE45544, GSE39262) para validação funcional da assinatura.
