Modelos gráficos espacialmente variables para redes de interacción célula-célula en imágenes tisulares multiplexadas

El microambiente tumoral (MAT) no es un telón de fondo pasivo para el crecimiento maligno, sino un ecosistema activo y espacialmente organizado en el que las poblaciones de células inmunitarias, estromales, vasculares y epiteliales participan en una interacción continua y dependiente de la ubicación, que moldea la vigilancia inmunitaria, la resistencia terapéutica y el pronóstico del paciente ([Hanahan, 2022]; [Hinshaw y Shevde, 2019]). Los estudios recientes que caracterizan el MAT mediante transcriptómica global y de células individuales han establecido la composición celular y el amplio repertorio de señalización en juego; lo que ha resultado más difícil de obtener es la estructura espacial de estas interacciones, es decir, qué tipos de células se co-localizan condicionalmente con otros en compartimentos tisulares específicos, y si estas dependencias espaciales difieren entre estados de enfermedad patológicamente distintos. Esta dimensión espacial no es meramente descriptiva. Existe evidencia convincente de que la organización geométrica de las células inmunitarias y estromales, el posicionamiento relativo de las poblaciones tumorales y efectoras, y la fragmentación o el acoplamiento espacial de los vecindarios inmunitarios son predictores independientes del resultado clínico en múltiples tipos de cáncer ([Schürch et al., 2020]; [Bhadury et al., 2026]; [Feng et al., 2023]).

Las tecnologías de imagen multiplexada han hecho posible, por primera vez, cuantificar la abundancia de tipos de células y la ubicación espacial simultáneamente en docenas de marcadores proteicos en secciones de tejido intacto. Las plataformas, que incluyen CODEX ([Black et al., 2021]), citometría de masas por imagen ([Giesen et al., 2014]) e inmunofluorescencia cíclica ([Lin et al., 2018]), ahora perfilan rutinariamente de 20 a 60 tipos de células por sección a resolución de célula única, generando conjuntos de datos en los que cada célula tiene tanto un fenotipo molecular como una coordenada espacial. Estas plataformas ya han generado información biológica que habría sido invisible para los enfoques basados en la disociación: la identificación de vecindarios celulares conservados en el cáncer colorrectal ([Schürch et al., 2020]), el descubrimiento de subtipos inmunitarios espacialmente distintos en el cáncer de mama triple negativo ([Gruosso et al., 2019]) y la demostración de que la proximidad espacial de las células T CD8+ a las células tumorales es un marcador pronóstico más potente que la abundancia de células T CD8+ por sí sola ([Feng et al., 2023]). A medida que estas plataformas se convierten en herramientas estándar en la oncología traslacional, la demanda de métodos estadísticos rigurosos para extraer redes de interacción condicionales de los datos espaciales resultantes está creciendo rápidamente.

La respuesta computacional a esta demanda ha seguido en gran medida dos vías que, de diferentes maneras, son insuficientes para el problema. La primera vía consiste en métodos de coocurrencia de tipos de células y enriquecimiento de vecindarios, incluidos los enfoques ampliamente utilizados implementados en Squidpy ([Palla et al., 2022]), Giotto ([Dries et al., 2021]), histoCAT ([Schapiro et al., 2017]) y herramientas relacionadas. Estos métodos caracterizan la estructura de proximidad espacial de los tipos de células mediante pruebas de permutación en matrices de recuento de vecindarios o estadísticas de coocurrencia bivariadas. Son computacionalmente ligeros y ampliamente aplicables, pero operan sobre relaciones marginales por pares y no tienen un mecanismo para la inferencia condicional. Debido a que no controlan la presencia de terceros tipos de células, una señal de coocurrencia entre, por ejemplo, macrófagos y células T reguladoras (Tregs), puede surgir de forma espúrea siempre que ambas poblaciones se enriquezcan conjuntamente en el mismo compartimento tisular junto con un tercer tipo de célula, como el estroma, que recluta a ambos de forma independiente. La independencia condicional, que es lo que realmente representa una arista de un modelo gráfico, no se puede recuperar a partir de la coocurrencia marginal.

La segunda vía consiste en métodos que modelan explícitamente las dependencias condicionales entre los tipos de células utilizando modelos gráficos. El análisis de varianza espacial (SVA; [Arnol et al., 2019]) descompone la variabilidad individual de la expresión de genes y proteínas en componentes intrínsecos, ambientales y de interacción entre células mediante un modelo de efecto aleatorio del proceso gaussiano, proporcionando estimaciones a nivel de gen de cuánta varianza de expresión es atribuible a la composición del vecindario celular. Si bien el SVA fue un importante avance metodológico al enmarcar el problema de la interacción espacial como un problema de estimación estadística en lugar de una prueba de proximidad, no recupera una red: asigna una fracción de varianza escalar a cada gen en lugar de estimar un gráfico disperso de dependencias condicionales entre los tipos de células. SpaCeNet ([Schrod et al., 2024]) extiende esta idea hacia la estimación de modelos gráficos modelando los potenciales de interacción por pares dependientes de la distancia entre los tipos de células utilizando la regularización de independencia condicional, pero las intensidades de interacción que infiere son resúmenes globales que no varían en las ubicaciones del tejido. Esta es una limitación importante, porque el MAT no es espacialmente homogéneo: las interacciones entre Tregs y macrófagos que caracterizan el margen invasor de un tumor pueden estar completamente ausentes en el núcleo del tumor, y un método que agrupa la evidencia en todas las ubicaciones del tejido para estimar un único coeficiente de interacción global caracterizará sistemáticamente de forma incorrecta la estructura de interacción en los tejidos heterogéneos. Patrones espacialmente informados para datos de inmunofluorescencia multiplexada (ISPAT; [Bhadury et al., 2026]) modela el perfil de interacción entre células que varía a lo largo de la trayectoria del tumor utilizando un modelo de efecto mixto del proceso gaussiano a través de la estructura de independencia condicional, pero la estimación de la covarianza se basa en una forma de análisis factorial multiestudio post hoc.

Independientemente de la literatura sobre la interacción espacial, se han desarrollado modelos gráficos bayesianos con priors de contracción continuos para datos ómicos no espaciales. El estimador de herradura gráfica de [Li et al. (2019)] coloca parámetros de contracción local de media Cauchy independientes en los elementos fuera de la diagonal de la matriz de precisión, proporcionando una inferencia dispersa adaptativa para las redes de coexpresión de genes. El enfoque de regresión a nivel de nodo de [Meinshausen y Bühlmann (2006)] descompone el problema del modelo gráfico p-dimensional en p problemas de regresión paralelizables, cada uno de los cuales es susceptible de una regularización estándar. Estos métodos se han aplicado eficazmente a la inferencia de redes genéticas y la integración multiómica ([Ni et al., 2022]), pero están diseñados para datos gaussianos multivariados estacionarios y no tienen un componente espacial: el coeficiente de regresión que vincula dos variables es un escalar, no una función de la ubicación del tejido. Aplicarlos a datos estructurados espacialmente de imágenes multiplexadas confunde el promedio espacial de un campo de interacción heterogéneo con una correlación parcial global, que no es ni la cantidad de interés científico ni la cantidad necesaria para recuperar el verdadero gráfico de independencia condicional.

La brecha que motiva este trabajo es, por lo tanto, específica y metodológica: no existe ningún método que estime simultáneamente (i) un gráfico de independencia condicional disperso sobre los tipos de células a partir de datos de imagen multiplexada, (ii) permita que cada arista de ese gráfico tenga una intensidad espacialmente variable que difiera en los microambientes del tejido y (iii) haga que la inclusión de aristas sea una decisión coherente de todo o nada a nivel del gráfico en lugar de un agregado de decisiones a nivel de coeficiente independientes. Los requisitos (i) y (ii) juntos son necesarios para recuperar redes de interacción condicionales biológicamente significativas a partir de datos del MAT espacialmente heterogéneos. El requisito (iii) es una necesidad estadística: cuando un campo de interacción espacialmente variable se representa mediante un conjunto de coeficientes de expansión de base, la aplicación de una contracción local independiente a cada coeficiente no produce un mecanismo para la decisión de arista de todo o nada que exige la selección de modelos gráficos, y la tasa de falsos descubrimientos resultante aumenta con el número de funciones de base utilizadas para representar el campo espacial.

Abordamos esta brecha a través de GP-GHS, un marco de regresión a nivel de nodo bayesiano que combina tres componentes. En primer lugar, modelamos cada coeficiente de regresión espacialmente variable como un proceso gaussiano utilizando la aproximación del proceso gaussiano del espacio de Hilbert (HSGP) de [Riutort-Mayol et al. (2023)], que reduce el costo de la inferencia exacta del GP a operaciones de matriz al representar el campo espacial en una base espectral. En segundo lugar, colocamos un prior de herradura de grupo en los coeficientes de la base espectral correspondientes a cada arista candidata, compartiendo un único parámetro de contracción local en todas las funciones de base dentro de un bloque. Esta estructura de grupo impone la decisión binaria a nivel de arista: cuando el parámetro de contracción del grupo se reduce a cero, todo el campo de interacción para ese par de tipos de células se reduce simultáneamente a cero en todas las ubicaciones espaciales, lo que corresponde a la ausencia de aristas; cuando es grande, el prior espectral que gobierna la suavidad espacial dentro del grupo es libre de operar. En tercer lugar, los p problemas de regresión a nivel de nodo se paralelizan y se combinan mediante la regla AND para recuperar un gráfico no dirigido simétrico, siguiendo el marco de selección de vecindarios de [Meinshausen y Bühlmann (2006)].

Aplicamos GP-GHS al conjunto de datos de imagen de tejido multiplexado CODEX de [Schürch et al. (2020)], que comprende 140 secciones de tejido de 35 pacientes con cáncer colorrectal (CCR) en etapa avanzada, estratificados por dos subtipos de microambiente tumoral patológicamente distintos: reacción similar a la enfermedad de Crohn (CLR) e infiltrado inflamatorio difuso (DII). El cáncer colorrectal exhibe una heterogeneidad espacial bien caracterizada en su arquitectura inmunitaria, y los subtipos CLR y DII están asociados con pronósticos sustancialmente diferentes ([Schürch et al., 2020]) y patrones de infiltración inmunitaria distintos. El conjunto de datos proporciona, por lo tanto, un caso de prueba ideal para determinar si GP-GHS puede recuperar redes de interacción espacial que difieran de forma coherente entre los grupos patológicos definidos. Evaluamos GP-GHS frente a cinco competidores mediante un estudio de simulación que abarca dos tamaños de problema y cuatro niveles de dispersión, utilizando gráficos sin escala generados por la unión preferencial para reflejar la topología dominada por nodos característica de las redes inmunitarias tumorales ([Barabási y Albert, 1999]; [Chen y Mellman, 2017]). Luego, aplicamos el método a los datos de CCR para identificar aristas diferencialmente activas entre CLR y DII, utilizando un modelo de efectos mixtos lineales sobre la puntuación de contracción posterior continua para tener en cuenta la correlación intra-paciente en las cuatro imágenes por paciente.

Las principales contribuciones de este trabajo son las siguientes. Desarrollamos el primer marco de modelo gráfico bayesiano para redes de interacción entre células espacialmente variables en imágenes de tejido multiplexadas, que combina una aproximación escalable del espacio de Hilbert del proceso gaussiano con un prior de herradura de grupo que impone la dispersión a nivel de gráfico al tiempo que aprovecha la fuerza espacial en todo el dominio del tejido. Demostramos mediante la simulación que la estructura de grupo del prior es el ingrediente fundamental para recuperar gráficos espacialmente estructurados: reemplazar la herradura de grupo con un prior de herradura escalar estándar hace que el rendimiento se desplome en todos los niveles de dispersión y tamaños de problema, independientemente de todas las demás opciones de modelado, lo que subraya que la contracción a nivel de coeficiente es fundamentalmente insuficiente para la agregación de señales espaciales. Aplicamos el método a un conjunto de datos de imagen de tejido de cáncer colorrectal que comprende 140 imágenes de 35 pacientes e identificamos una red inmunosupresora centrada en Tregs que se amplifica significativamente en el subtipo de infiltrado inflamatorio difuso en relación con el subtipo de respuesta similar a la enfermedad de Crohn, un hallazgo que se basa en la biología espacial conocida de estos dos subtipos y se recupera con alta confianza estadística después de la corrección de múltiples pruebas a nivel de paciente. Modelar la puntuación de contracción posterior continua κ a través de un modelo de efectos mixtos lineales con efectos aleatorios del paciente proporciona una potencia sustancialmente mayor para las pruebas diferenciales, un resultado con implicaciones prácticas directas para los estudios que tienen como objetivo comparar redes de interacción espacial entre poblaciones de pacientes o condiciones experimentales.

Métodos

Estructura de datos y notación

Sea el conjunto de ubicaciones espaciales observadas (puntos de tejido) dentro de una única sección de tejido, donde cada ubicación representa una coordenada bidimensional en el espacio tisular físico. En cada ubicación , observamos valores de expresión normalizados para tipos de células, recopilados en la matriz , donde denota la expresión del tipo de célula en la ubicación . El objetivo es recuperar un grafo no dirigido disperso sobre los tipos de células, donde denota el conjunto de vértices y una arista indica una dependencia condicional espacialmente estructurada entre los tipos de células y después de tener en cuenta todos los tipos de células restantes. De manera crucial, permitimos que la intensidad de esta dependencia varíe continuamente a lo largo del tejido, una característica que distingue el marco propuesto de los modelos gráficos clásicos.

Marco de Regresión Espacial a Nivel de Nodo

Adoptamos la estrategia de selección de vecindad de [Meinshausen y Bühlmann (2006)], extendida para dar cabida a coeficientes de regresión que varían espacialmente. Para cada tipo de célula , formulamos la regresión a nivel de nodo:

donde es la expresión del tipo de célula en la ubicación , es la expresión de su vecino putativo , y es el coeficiente de regresión que varía espacialmente y que codifica la asociación condicional entre los tipos de células y en la ubicación . Se declara presente una arista entre y si la función no es idénticamente cero en el dominio del tejido. Esta formulación difiere de una regresión lasso estándar a nivel de nodo en que los coeficientes son funciones del espacio en lugar de escalares, lo que permite que la intensidad y la dirección de las interacciones entre células varíen en los microambientes del tejido, como el núcleo del tumor, el margen invasivo y el estroma.

Las regresiones a nivel de nodo son independientes y se paralelizan en varios núcleos durante la computación. Se recupera un grafo no dirigido final aplicando la regla AND: una arista () se incluye en si y solo si tanto la regresión de sobre como la regresión de sobre dan como resultado una función de coeficiente activa. La regla AND es conservadora en relación con la regla OR, pero controla las falsas detecciones de manera más confiable en el régimen moderado típico de los datos de imagen multiplexada ([Meinshausen y Bühlmann, 2006]).

Coeficientes que Varían Espacialmente a través de Procesos Gaussianos en el Espacio de Hilbert

Cada coeficiente que varía espacialmente se modela como una realización de un proceso gaussiano (GP) sobre el dominio del tejido. Específicamente, colocamos la distribución a priori:

donde es el núcleo de covarianza de Matérn con el parámetro de suavidad y la longitud de escala . La familia de Matérn se prefiere sobre el núcleo exponencial al cuadrado porque permite un control explícito sobre la diferenciabilidad del campo espacial: para , el proceso resultante es una vez diferenciable en el sentido del cuadrado medio, lo cual es apropiado para las interacciones espaciales entre células que son suaves dentro de los compartimentos del tejido, pero que pueden presentar transiciones bruscas en los límites de los compartimentos ([Stein, 1999]). Para dos ubicaciones espaciales y con distancia euclidiana , el núcleo de Matérn 3/2 toma la forma:

donde es la varianza marginal y controla el rango espacial de la dependencia. La longitud de escala se establece en la mediana de las distancias por pares entre las coordenadas espaciales normalizadas por defecto, aunque puede ser especificada por el analista cuando se dispone de conocimientos previos sobre el rango de interacción.

La inferencia exacta de GP requiere computación debido a la factorización de Cholesky de la matriz de covarianza n × n, lo cual es prohibitivo para las secciones de tejido con n > 500 puntos. Por lo tanto, empleamos la aproximación del proceso gaussiano en el espacio de Hilbert (HSGP) de [Riutort-Mayol et al. (2023)], que aproxima el GP utilizando las autofunciones del operador Laplaciano en un dominio acotado. Específicamente, para el dominio bidimensional obtenido al normalizar las coordenadas espaciales, las autofunciones del Laplaciano se factorizan como productos tensoriales de funciones base sinusoidales unidimensionales:

lo que da como resultado funciones base en total. La aproximación HSGP representa entonces el campo espacial como:

donde los pesos de la base heredan una distribución a priori gaussiana de media cero con varianza dada por la densidad espectral del núcleo de Matérn evaluada en la eigenfrecuencia correspondiente , con . Para el núcleo de Matérn 3/2 en d = 2 dimensiones, la densidad espectral es:

de modo que . Esta distribución a priori sobre los pesos de la base es lo que codifica la suavidad espacial: las funciones base de baja frecuencia (pequeña ) reciben una gran varianza a priori y son libres de capturar tendencias espaciales amplias, mientras que las funciones base de alta frecuencia (grandes ) se reducen fuertemente hacia cero, suprimiendo la variación espacial brusca. En todo momento, fijamos en , ya que la varianza marginal de cada campo de interacción se controla mediante el producto en la jerarquía de herradura de grupo; incluir un libre en la distribución a priori crearía una redundancia no identificada entre la varianza marginal del GP y la escala de reducción local del grupo. La aproximación HSGP reduce la complejidad computacional a para una regresión de nodo dada, lo que la hace factible para las dimensiones de datos encontradas en la imagen multiplexada.

Distribución a Priori de Herradura de Grupo para la Selección de Aristas

La aproximación HSGP transforma el modelo de coeficiente que varía espacialmente (1) en una regresión lineal sobre la matriz de diseño expandida. Para el nodo (tipo de célula) con vecinos (otros tipos de células), defina la matriz de diseño en bloques cuyo -ésimo bloque de columnas es:

donde es el vector de expresión del -ésimo tipo de célula vecino y es la matriz de funciones base evaluadas. El -ésimo bloque de los coeficientes de regresión recopila los pesos de la base espectral para el vecino , de modo que el campo de interacción espacial se recupera como .

Una elección de modelado fundamental concierne a la ubicación de la distribución a priori de dispersión dentro de esta regresión -dimensional. El enfoque ingenuo de colocar parámetros de reducción local independientes en cada uno de los predictores, como en una herradura escalar estándar ([Carvalho et al., 2010]), es inapropiado por dos razones que son específicas de esta estructura del problema. Primero, las funciones base dentro de un bloque no son científicamente significativas individualmente; son un dispositivo numérico para aproximar una única función espacial . Reducir los coeficientes de la base individualmente de forma independiente rompe esta conexión y puede producir estimaciones incoherentes donde un campo de interacción espacial se reduce parcialmente a través de los componentes de frecuencia sin ninguna interpretación geométrica; por ejemplo, reducir los componentes de alta frecuencia mientras se conservan los de baja frecuencia no corresponde a ninguna declaración biológica significativa sobre si el tipo de célula influye en en algunas regiones del tejido, pero no en otras. En segundo lugar, la selección de aristas en un modelo gráfico es inherentemente una decisión binaria a nivel de grupo: o toda la función de interacción es cero en todas partes del tejido (arista ausente) o es distinta de cero en al menos alguna región (arista presente). Una distribución a priori escalar independiente no proporciona ningún mecanismo para esta decisión de nivel de bloque y, en su lugar, acumula oportunidades independientes para declarar falsamente la actividad, lo que infla las tasas de falsos descubrimientos en proporción al número de funciones base, teniendo en cuenta que si la varianza de no es proporcional a para todos dentro de , entonces no se tiene la correspondencia del proceso gaussiano.

La resolución adecuada es compartir un único parámetro de reducción local en todos los coeficientes de la base dentro del bloque , siguiendo el marco de regularización agrupada de [Xu et al. (2016)]. Cuando , todo el bloque se reduce simultáneamente a cero, de modo que el campo espacial reconstruido es idénticamente cero en todo el dominio del tejido; la arista está ausente globalmente. Cuando es grande, todos los pesos se liberan simultáneamente y el patrón espacial dentro del grupo está gobernado enteramente por la distribución a priori espectral . Esta separación de roles - controla la existencia de la arista, controla la suavidad espacial condicional a la existencia - es la propiedad estructural clave de la distribución a priori propuesta y no se puede lograr mediante ninguna distribución a priori que asigne una reducción local independiente a las funciones base individuales.

Concretamente, el -ésimo elemento del -ésimo bloque recibe la distribución a priori jerárquica:

donde denota la distribución de media Cauchy ([Polson y Scott, 2012]). El parámetro global controla la dispersión general del grafo; un pequeño empuja todas las aristas hacia la ausencia simultáneamente, lo que es consistente con la expectativa de que los grafos de interacción entre células en el tejido sean dispersos. El parámetro local del grupo controla la existencia de la arista individual mediante la modulación de la varianza de todo el -ésimo bloque a través de un escalar común. La densidad espectral opera entonces dentro de cada grupo superviviente para hacer cumplir la regularidad espacial, dando una mayor varianza a priori a los componentes de baja frecuencia (suaves) de y una menor varianza a los componentes de alta frecuencia (bruscos). Los tres niveles de la jerarquía operan así en escalas estructuralmente ortogonales: dispersión a nivel de grafo, existencia de la arista y regularidad espacial dentro de la arista. Esta parametrización preserva la correspondencia del GP: la varianza a priori de es , lo que se reduce a la distribución a priori espectral de HSGP cuando . La estructura del grupo modula la escala general de todo el bloque a través del escalar , dejando que la ponderación de la frecuencia dentro del grupo sea determinada únicamente por . En consecuencia, cuando todo el bloque se colapsa proporcionalmente en todas las frecuencias, manteniendo la estructura del GP del campo superviviente en lugar de distorsionarla.

La estructura del grupo también tiene una consecuencia directa sobre cómo se acumula la información en las ubicaciones espaciales. Al evaluar si la arista está activa, la distribución a posteriori para integra la evidencia de todas las observaciones y todos los coeficientes de la base de forma conjunta a través de la forma cuadrática ponderada por el GP , que aparece en la distribución completa condicional para (Sección 2.5). Una interacción espacialmente suave pero globalmente débil -una en la que es distinta de cero pero pequeña en todas partes- acumulará suficiente señal en todos los componentes de la base para evitar que se colapse, mientras que una señal que solo esté presente en coeficientes de base individuales aislados se reducirá según la distribución a priori del grupo, ya que agrega sobre todo el bloque. Esta propiedad de agregación es lo que distingue una distribución a priori de grupo de las distribuciones a priori escalares independientes y es el mecanismo que le da al método propuesto su ventaja de sensibilidad sobre los métodos que ignoran por completo la estructura espacial.

El factor de reducción del grupo para la arista se define como , que se encuentra en (0, 1). Cuando el grupo se reduce completamente a cero y la arista está ausente; cuando el grupo no se reduce y la arista está presente. Declaramos que la arista está activa si la media a posteriori , el umbral en el que la distribución a posteriori asigna una probabilidad igual a la señal y al ruido para ese grupo ([Carvalho et al., 2010]).

La función de verosimilitud para el nodo es:

y recibe la distribución a priori dispersa-IG, que proporciona una regularización suave al tiempo que se mantiene difusa en un amplio rango de varianzas residuales.

Inferencia a Posteriori a través del Muestreo de Gibbs

Todas las distribuciones completas condicionales están disponibles en forma cerrada, lo que da como resultado un muestreador de Gibbs en bloques eficiente. Introduciendo las variables auxiliares y para las representaciones de media Cauchy y , el muestreador recorre las siguientes actualizaciones.

Actualización
. La distribución completa condicional es gaussiana multivariada:

donde donde es diagonal en bloques con bloque siendo la covarianza a priori diagonal en bloques, con las densidades espectrales que determinan la escala dentro del grupo.

Actualización
. Para cada , la distribución completa condicional es:

donde es la forma cuadrática ajustada por el GP. Esta norma ponderada tiene en cuenta la estructura espectral: los coeficientes asociados con las funciones base de baja frecuencia (suaves) se penalizan menos que los de alta frecuencia.

Actualización
. El parámetro de reducción global satisface:

Actualización
. La actualización de la varianza residual es:

Las variables auxiliares y se actualizan como y . Ejecutamos el muestreador durante 3.000 iteraciones con un período de adaptación de 1.000 y un adelgazamiento de 5, conservando 400 muestras a posteriori por regresión de nodo. Los p muestreadores de Gibbs a nivel de nodo son fácilmente paralelizables y se distribuyen en los núcleos de CPU disponibles.

Mapas de Aristas Espaciales y Resúmenes a Posteriori

Para cada arista activa declarada , el campo de interacción espacial medio a posteriori se recupera como:

donde es la media a posteriori del -ésimo bloque de coeficientes. Esto produce un mapa continuo sobre el dominio del tejido, asignando una fuerza de interacción con signo en cada ubicación observada. Las regiones donde es grande y positivo indican una co-enriquecimiento espacial de los tipos de células y , mientras que los valores negativos indican exclusión espacial. Los intervalos creíbles a posteriori del 95% por puntos se calculan a partir de las muestras a posteriori, proporcionando una cuantificación de la incertidumbre para los mapas espaciales. La colección completa de mapas de aristas , junto con la matriz de adyacencia global y los resúmenes de reducción a nivel de grupo, constituyen la salida del método propuesto.

Construcción del grafo.

Las p regresiones a nivel de nodo se ajustan de forma independiente, cada una produciendo una lista de aristas dirigidas: el nodo declara que el nodo vecino está activo si , donde denota el factor de reducción a posteriori a nivel de grupo de la regresión de sobre . Debido a que las regresiones a nivel de nodo no están restringidas conjuntamente, los conjuntos de aristas dirigidas no tienen por qué ser simétricos: es posible que la regresión de sobre declare que la arista está activa, mientras que la regresión de sobre no lo hace. Resolvemos esta asimetría mediante la regla AND: una arista no dirigida se incluye en el grafo estimado si y solo si ambas regresiones la declaran activa.

La regla AND es conservadora en relación con la regla OR (que incluiría una arista si alguna de las regresiones la declara activa) y se sabe que proporciona un mejor control de los falsos descubrimientos en el régimen moderado ([Meinshausen y Bühlmann, 2006]). La matriz de adyacencia resultante es simétrica por construcción, con para todos los pares. La fuerza de la arista para una arista retenida se resume mediante , tomando el posterior direccional más confiable como la puntuación de reducción representativa. cercano a cero indica una señal activa y , cercano a uno, indica una reducción hacia el valor nulo. Tomar el mínimo selecciona el posterior direccional que expresa la evidencia más sólida de interacción, que es el resumen apropiado para una arista no dirigida simétrica. La colección completa de mapas de interacción espacial para las aristas activas, la matriz de adyacencia simétrica y la matriz de puntuaciones de reducción a posteriori constituyen la salida completa del método. Presentamos una visión general completa del método en la Figura 1.

Estructura de datos y notación

Sea denotar el conjunto de ubicaciones espaciales observadas (puntos de tejido) dentro de una sola sección de tejido, donde cada ubicación representa una coordenada bidimensional en el espacio físico del tejido. En cada ubicación , observamos valores de expresión normalizados para tipos de células, recopilados en la matriz , donde denota la expresión del tipo de célula en la ubicación . El objetivo es recuperar un grafo no dirigido disperso sobre los tipos de células, donde denota el conjunto de vértices y una arista indica una dependencia condicional espacialmente estructurada entre los tipos de células y después de tener en cuenta todos los tipos de células restantes. Fundamentalmente, permitimos que la fuerza de esta dependencia varíe continuamente en todo el tejido, una característica que distingue el marco propuesto de los modelos gráficos clásicos.

Marco de regresión espacial a nivel de nodo

Adoptamos la estrategia de selección de vecindad de [Meinshausen y Bühlmann (2006)], extendida para dar cabida a los coeficientes de regresión espacialmente variables. Para cada tipo de célula , formulamos la regresión a nivel de nodo

donde es la expresión del tipo de célula en la ubicación , es la expresión de su posible vecino y es el coeficiente de regresión espacialmente variable que codifica la asociación condicional entre los tipos de células y en la ubicación . Se declara que existe una arista entre y si la función no es idénticamente cero en todo el dominio del tejido. Esta formulación difiere de una regresión lasso a nivel de nodo estándar en que los coeficientes son funciones del espacio en lugar de escalares, lo que permite que la fuerza y la dirección de las interacciones entre células varíen en todo el tejido.

Las regresiones a nivel de nodo son independientes y se paralelizan en varios núcleos durante el cálculo. Se recupera un grafo no dirigido final aplicando la regla AND: una arista () se incluye en si y solo si tanto la regresión de sobre como la regresión de sobre dan como resultado una función de coeficiente activa. La regla AND es conservadora en relación con la regla OR, pero controla los falsos descubrimientos de forma más fiable en el régimen moderado típico de los datos de imagen multiplexada ([Meinshausen y Bühlmann, 2006]).

Coeficientes espacialmente variables a través de procesos gaussianos del espacio de Hilbert

Cada coeficiente espacialmente variable se modela como una realización de un proceso gaussiano (GP) sobre el dominio del tejido. Específicamente, colocamos la distribución a priori

donde es el núcleo de covarianza de Matérn con el parámetro de suavidad y la longitud de escala . La familia Matérn se prefiere al núcleo exponencial al cuadrado porque permite un control explícito sobre la diferenciabilidad del campo espacial: para , el proceso resultante es una vez diferenciable en el sentido del cuadrado medio, lo que es apropiado para las interacciones espaciales entre células que son suaves dentro de los compartimentos del tejido, pero que pueden presentar transiciones bruscas en los límites de los compartimentos ([Stein, 1999]). Para dos ubicaciones espaciales y con distancia euclidiana , el núcleo Matérn 3/2 toma la forma

donde es la varianza marginal y controla el rango espacial de la dependencia. La longitud de escala se establece en la mediana de la distancia por pares entre las coordenadas espaciales normalizadas por defecto, aunque puede ser especificada por el analista cuando se dispone de conocimientos previos sobre el rango de interacción.

La inferencia exacta de GP requiere computación debido a la factorización de Cholesky de la matriz de covarianza n × n, lo que es prohibitivo para las secciones de tejido con n > 500 puntos. Por lo tanto, empleamos la aproximación del proceso gaussiano del espacio de Hilbert (HSGP) de [Riutort-Mayol et al. (2023)], que aproxima el GP utilizando las autofunciones del operador Laplaciano en un dominio acotado. Específicamente, para el dominio bidimensional obtenido al normalizar las coordenadas espaciales, las autofunciones del Laplaciano se factorizan como productos tensoriales de funciones base sinusoidales unidimensionales,

lo que da un total de funciones base. La aproximación HSGP representa entonces el campo espacial como

donde los pesos de la base heredan una distribución a priori gaussiana de media cero con varianza dada por la densidad espectral del núcleo de Matérn evaluada en la frecuencia propia correspondiente , con . Para el núcleo Matérn 3/2 en d = 2 dimensiones, la densidad espectral es

de modo que . Esta distribución a priori sobre los pesos de la base es lo que codifica la suavidad espacial: las funciones base de baja frecuencia (pequeña ) reciben una gran varianza a priori y son libres de capturar tendencias espaciales amplias, mientras que las funciones base de alta frecuencia (grandes ) se reducen fuertemente hacia cero, suprimiendo la variación espacialmente rugosa. En todo momento, fijamos en , ya que la varianza marginal de cada campo de interacción está controlada por el producto en la jerarquía de herradura de grupo; absorber un libre en la distribución a priori crearía una redundancia no identificada entre la varianza marginal del GP y la escala de reducción local del grupo. La aproximación HSGP reduce la complejidad computacional a para una regresión de nodo dada, lo que la hace factible para las dimensiones de datos encontradas en la imagen multiplexada.

Distribución a priori de herradura de grupo para la selección de aristas

La aproximación HSGP transforma el modelo de coeficiente espacialmente variable (1) en una regresión lineal sobre la matriz de diseño ampliada. Para el nodo (tipo de célula) con vecinos (otros tipos de células), defina la matriz de diseño bloqueada cuyo -ésimo bloque de columnas es

donde es el vector de expresión del -ésimo tipo de célula vecino y es la matriz de funciones base evaluadas. El -ésimo bloque de coeficientes de regresión recopila los pesos de la base espectral para el vecino , de modo que el campo de interacción espacial se recupera como .

Una elección de modelado fundamental concierne a la colocación de la distribución a priori de dispersión dentro de esta regresión -dimensional. El enfoque ingenuo de colocar parámetros de reducción locales independientes en cada uno de los predictores, como en una herradura escalar estándar ([Carvalho et al., 2010]), es inapropiado por dos razones que son específicas de esta estructura del problema. Primero, las funciones base dentro de un bloque no tienen un significado científico individual; son un dispositivo numérico para aproximar una sola función espacial . Reducir los coeficientes de la base individualmente de forma independiente rompe esta conexión y puede producir estimaciones incoherentes donde un campo de interacción espacial se reduce parcialmente en los componentes de frecuencia sin ninguna interpretación geométrica; por ejemplo, reducir los componentes de alta frecuencia mientras se retienen los componentes de baja frecuencia no corresponde a ninguna declaración biológica significativa sobre si el tipo de célula influye en en algunas regiones del tejido, pero no en otras. En segundo lugar, la selección de aristas en un modelo gráfico es inherentemente una decisión binaria a nivel de grupo: o toda la función de interacción es cero en todas partes del tejido (arista ausente) o es distinta de cero en al menos alguna región (arista presente). Una distribución a priori escalar independiente no proporciona ningún mecanismo para esta decisión de todo o nada a nivel de bloque y, en su lugar, acumula oportunidades independientes para declarar falsamente la actividad, lo que infla las tasas de falsos descubrimientos en proporción al número de funciones base, teniendo en cuenta que si la varianza de no es proporcional a para todos dentro del , entonces no se tiene la correspondencia del proceso gaussiano.

La resolución adecuada es compartir un único parámetro de reducción local en todos los coeficientes de la base dentro del bloque , siguiendo el marco de regularización agrupada de [Xu et al. (2016)]. Cuando , todo el bloque se reduce simultáneamente a cero, de modo que el campo espacial reconstruido es idénticamente cero en todo el dominio del tejido: la arista está ausente globalmente. Cuando es grande, todos los pesos se liberan simultáneamente y el patrón espacial dentro del grupo está gobernado enteramente por la distribución a priori espectral . Esta separación de roles - controla la existencia de la arista, controla la suavidad espacial condicional a la existencia - es la propiedad estructural clave de la distribución a priori propuesta y no se puede lograr mediante ninguna distribución a priori que asigne una reducción local independiente a las funciones base individuales.

Concretamente, el -ésimo elemento del -ésimo bloque recibe la distribución a priori jerárquica

donde denota la distribución de media de Cauchy ([Polson y Scott, 2012]). El parámetro global controla la dispersión del grafo general: un pequeño empuja todas las aristas hacia la ausencia simultáneamente, lo que es consistente con la expectativa de que los grafos de interacción entre células en el tejido sean dispersos. El parámetro local del grupo controla la existencia de la arista individual mediante la modulación de la varianza de todo el -ésimo bloque a través de un escalar común. La densidad espectral opera entonces dentro de cada grupo superviviente para hacer cumplir la regularidad espacial, dando una mayor varianza a priori a los componentes de baja frecuencia (suaves) de y una menor varianza a los componentes de alta frecuencia (rugosos). Los tres niveles de la jerarquía operan así en escalas estructuralmente ortogonales: dispersión a nivel de grafo, existencia de aristas a nivel de arista y regularidad espacial dentro de la arista. Esta parametrización preserva la correspondencia del GP: la varianza a priori de es , que se reduce a la distribución a priori espectral del HSGP cuando . La estructura del grupo modula la escala general de todo el bloque a través del escalar , dejando que la ponderación de la frecuencia dentro del grupo esté determinada únicamente por . En consecuencia, cuando todo el bloque se colapsa proporcionalmente en todas las frecuencias, manteniendo la estructura del GP del campo superviviente en lugar de distorsionarla.

La estructura del grupo también tiene una consecuencia directa en la forma en que se acumula la información en diferentes ubicaciones espaciales. Al evaluar si un borde está activo, la distribución a posteriori de integra la evidencia de todas las observaciones y todos los coeficientes de base de forma conjunta a través de la forma cuadrática ponderada por el proceso gaussiano (GP), que aparece en la distribución condicional completa de (Sección 2.5). Una interacción espacialmente suave pero globalmente débil (es decir, una en la que es diferente de cero pero pequeña en todas partes) acumulará suficiente señal en todos los componentes de la base para evitar que colapse, mientras que una señal que solo está presente en coeficientes de base individuales aislados se reducirá mediante la restricción del grupo, ya que se agrega sobre todo el bloque. Esta propiedad de agregación es precisamente lo que distingue una restricción de grupo de las restricciones escalares independientes y es el mecanismo que le da al método propuesto su ventaja en sensibilidad sobre los métodos que ignoran por completo la estructura espacial.

El factor de reducción del grupo para el borde se define como , que se encuentra en el intervalo (0, 1). Cuando , el grupo se reduce por completo a cero y el borde está ausente; cuando , el grupo no se reduce y el borde está presente. Declaramos que el borde está activo si la media a posteriori , el umbral en el que la distribución a posteriori asigna una probabilidad igual a la señal y al ruido para ese grupo ([Carvalho et al., 2010]).

La función de verosimilitud para el nodo es

y recibe la restricción IG dispersa, que proporciona una regularización suave al tiempo que se mantiene difusa en un amplio rango de varianzas residuales.

Inferencia a posteriori mediante muestreo de Gibbs

Todas las distribuciones condicionales completas están disponibles en forma cerrada, lo que da como resultado un muestreador de Gibbs de bloques eficiente. Introduciendo las variables auxiliares y para las representaciones de la media de Cauchy, el muestreador itera a través de las siguientes actualizaciones.

Actualización
. La distribución condicional completa es gaussiana multivariada,

donde , donde es diagonal en bloques, siendo el bloque la matriz de covarianza a priori diagonal en bloques, con las densidades espectrales que determinan la escala dentro del grupo.

Actualización
. Para cada , la distribución condicional completa es

donde es la forma cuadrática ajustada por el GP. Esta norma ponderada tiene en cuenta la estructura espectral: los coeficientes asociados con las funciones de base de baja frecuencia (suaves) se penalizan menos que los de alta frecuencia.

Actualización
. El parámetro de reducción global satisface

Actualización
. La actualización de la varianza residual es

Las variables auxiliares y se actualizan como y . Ejecutamos el muestreador durante 3.000 iteraciones con un período de "burn-in" de 1.000 y un adelgazamiento de 5, conservando 400 muestras a posteriori por regresión de nodo. Los muestreadores nodales p son fácilmente paralelizables y se distribuyen entre los núcleos de CPU disponibles.

Mapas de bordes espaciales y resúmenes a posteriori

Para cada borde declarado activo , el campo de interacción espacial a posteriori se recupera como

donde es la media a posteriori del -ésimo bloque de coeficientes. Esto produce un mapa continuo sobre el dominio del tejido, asignando una fuerza de interacción con signo en cada ubicación observada. Las regiones donde es grande y positivo indican una co-enriquecimiento espacial de los tipos de células y , mientras que los valores negativos indican exclusión espacial. Los intervalos de credibilidad a posteriori del 95% por punto se calculan a partir de las muestras a posteriori, proporcionando una cuantificación de la incertidumbre para los mapas espaciales. La colección completa de mapas de bordes, junto con la matriz de adyacencia global y los resúmenes de reducción a nivel de grupo, constituyen la salida del método propuesto.

Construcción del grafo.

Las regresiones nodales p se ajustan de forma independiente, cada una produciendo una lista de bordes dirigida: el nodo declara que el vecino está activo si , donde denota la media a posteriori del factor de reducción del grupo de la regresión de sobre . Debido a que las regresiones nodales no están restringidas conjuntamente, los conjuntos de bordes dirigidos no tienen por qué ser simétricos: es posible que la regresión de sobre declare que el borde está activo, mientras que la regresión de sobre no lo hace. Resolvemos esta asimetría mediante la regla AND: un borde no dirigido se incluye en el grafo estimado si y solo si ambas regresiones lo declaran activo.

La regla AND es conservadora en relación con la regla OR (que incluiría un borde si alguna de las regresiones lo declara activo) y se sabe que proporciona un mejor control del falso descubrimiento en el régimen moderado ([Meinshausen y Bühlmann, 2006]). La matriz de adyacencia resultante es simétrica por construcción, con para todos los pares. La fuerza del borde para un borde retenido se resume mediante , tomando el de los dos valores a posteriori direccionales que sea más confiable como la puntuación de reducción representativa. cerca de cero indica una señal activa y , cerca de uno, indica una reducción hacia el valor nulo. Tomar el mínimo selecciona el valor a posteriori direccional que expresa la evidencia más sólida de interacción, que es el resumen apropiado para un borde no dirigido simétrico. La colección completa de mapas de interacción espacial para los bordes activos, la matriz de adyacencia simétrica y la matriz de puntuaciones de reducción a posteriori constituyen la salida completa del método. Presentamos una descripción esquemática completa del método en la Figura 1.

Construcción del grafo.

Las regresiones nodales p se ajustan de forma independiente, cada una produciendo una lista de bordes dirigida: el nodo declara que el vecino está activo si , donde denota la media a posteriori del factor de reducción del grupo de la regresión de sobre . Debido a que las regresiones nodales no están restringidas conjuntamente, los conjuntos de bordes dirigidos no tienen por qué ser simétricos: es posible que la regresión de sobre declare que el borde está activo, mientras que la regresión de sobre no lo hace. Resolvemos esta asimetría mediante la regla AND: un borde no dirigido se incluye en el grafo estimado si y solo si ambas regresiones lo declaran activo.

La regla AND es conservadora en relación con la regla OR (que incluiría un borde si alguna de las regresiones lo declara activo) y se sabe que proporciona un mejor control del falso descubrimiento en el régimen moderado ([Meinshausen y Bühlmann, 2006]). La matriz de adyacencia resultante es simétrica por construcción, con para todos los pares. La fuerza del borde para un borde retenido se resume mediante , tomando el de los dos valores a posteriori direccionales que sea más confiable como la puntuación de reducción representativa. cerca de cero indica una señal activa y , cerca de uno, indica una reducción hacia el valor nulo. Tomar el mínimo selecciona el valor a posteriori direccional que expresa la evidencia más sólida de interacción, que es el resumen apropiado para un borde no dirigido simétrico. La colección completa de mapas de interacción espacial para los bordes activos, la matriz de adyacencia simétrica y la matriz de puntuaciones de reducción a posteriori constituyen la salida completa del método. Presentamos una descripción esquemática completa del método en la Figura 1.

Estudio de simulación

Mecanismo de generación de datos

Evaluamos el método propuesto en comparación con cuatro competidores a través de un estudio de simulación estructurado diseñado para reflejar las características clave de los datos de imagen de tejidos multiplexados. Consideramos un conjunto fijo de ubicaciones espaciales n = 600 extraídas uniformemente del cuadrado unitario [0, 10]2, con p = 15, 25 tipos de células. El grafo de interacción célula-célula verdadero se genera utilizando el modelo de apego preferencial de Barabási-Albert ([Barabási y Albert, 1999]), que produce una topología sin escala con un pequeño número de tipos de células "hub" que tienen muchas conexiones y muchos tipos de células que tienen pocas conexiones. Esta arquitectura está biológicamente motivada: en los microambientes tumorales, un puñado de tipos de células centrales, como los macrófagos y las células tumorales, tienden a coordinar las interacciones con muchos otros, mientras que las poblaciones más periféricas, como las células plasmáticas o los mastocitos, participan en menos interacciones ([Chen y Mellman, 2017]). Examinamos cuatro niveles de dispersión, apuntando a probabilidades de borde de 0,7, 0,5, 0,3 y 0,1, correspondientes a grafos verdaderos densos, moderados, dispersos y muy dispersos.

Condicionalmente a la matriz de adyacencia verdadera , los datos se generan de acuerdo con la ecuación estructural

donde y los verdaderos coeficientes espacialmente variables se extraen de un GP de media cero con kernel de Matérn 3/2, longitud de escala apropiada y desviación estándar marginal . Esta elección de en relación con el dominio normalizado produce patrones espaciales que varían en una fracción moderada de la extensión del tejido, lo que es consistente con la escala de los patrones de infiltración inmune observados en los datos de imagen multiplexados. Cada réplica de simulación extrae nuevas realizaciones de GP para todos los bordes verdaderos y nuevo ruido, lo que garantiza que los resultados reflejen la variabilidad en ambas configuraciones espaciales y realizaciones de ruido. Se generan diez réplicas independientes por nivel de dispersión, lo que da como resultado un total de 40 conjuntos de datos simulados.

Métodos competidores

Comparamos GP-GHS con cinco métodos competidores que abarcan la verosimilitud penalizada, la reducción escalar bayesiana, la regresión GP exacta y el umbral de correlación ingenuo.

Graphical Lasso (GLasso).

El graphical lasso ([Friedman et al., 2008]) estima una matriz de precisión dispersa resolviendo , donde es la matriz de covarianza de la muestra de () y es un parámetro de penalización seleccionado mediante BIC en una cuadrícula de 15 valores. GLasso asume la estacionariedad e ignora las coordenadas espaciales, produciendo una única matriz de precisión global.

Nodewise Lasso (NL).

Siguiendo a [Meinshausen y Bühlmann (2006)], ajustamos una regresión penalizada - para cada nodo frente a todos los demás nodos, con la penalización seleccionada mediante validación cruzada de cinco pliegues. El grafo no dirigido se recupera mediante la regla AND. Al igual que GLasso, NL trata todos los coeficientes de regresión como escalares y no tiene en cuenta la estructura espacial.

Standard Horseshoe (SHS).

Este método aplica la restricción de herradura escalar de [Carvalho et al. (2010)] a cada regresión nodal, asignando un parámetro de reducción local independiente a cada predictor escalar en lugar de uno por bloque de vecinos. Debido a que la reducción opera a nivel de coeficiente en lugar de a nivel de grupo, SHS no tiene un mecanismo para tomar decisiones coherentes de borde "todo o nada": los coeficientes de base individuales dentro de un bloque de vecinos pueden reducirse selectivamente, lo que da como resultado estimaciones del campo de interacción que no son uniformemente cero ni uniformemente diferentes de cero. La inclusión de bordes se evalúa mediante el criterio de pseudo-p-valor aplicado a cada coeficiente escalar. SHS conserva el marco bayesiano nodal de GP-GHS, pero descarta la estructura de grupo y la restricción espectral del GP, aislando así su contribución conjunta al rendimiento.

Pairwise GP Regression (Exact GP –EGP).

Este competidor ajusta directamente una regresión independiente del proceso gaussiano para cada par de nodos ordenados , utilizando un kernel de Matérn-3/2 con una longitud de escala que coincide con la auto-seleccionada de HSGP. Para cada nodo de respuesta , el residuo parcial con respecto al vecino se forma regresionando todos los demás vecinos mediante OLS, y el residuo resultante se trata como una observación ruidosa de . Llamamos a este competidor "pairwise" porque no puede ajustar simultáneamente las restricciones del GP para todos los vecinos de un nodo determinado. Hacerlo requeriría invertir una matriz de covarianza del GP conjunta de dimensión , lo que es computacionalmente prohibitivo en los tamaños de problema considerados. El enfoque pairwise utiliza en su lugar los residuos OLS de un solo paso para eliminar los vecinos restantes, reemplazando un ajuste del GP conjunto con una secuencia de ajustes del GP marginal a un costo significativamente menor. La inclusión de bordes se determina mediante votación por mayoría: un borde se declara activo si el intervalo de credibilidad a posteriori del 95% excluye cero en más del 50% de las ubicaciones espaciales, y el grafo no dirigido final se simetriza mediante la regla AND. EGP es la línea de base exacta natural para GP-GHS: comparte la misma restricción espacial y la lógica de selección de bordes, pero omite la aproximación HSGP, la reducción de la herradura de grupo y la estimación conjunta de múltiples vecinos, a un costo computacional significativamente mayor.

Correlation Threshold (CT).

Una línea de base ingenua que declara un borde siempre que en la matriz de correlación de la muestra. CT no hace suposiciones estructurales y sirve como un límite de rendimiento inferior.

GLasso y NL representan métodos estándar de verosimilitud penalizada para la selección de modelos gráficos. SHS aísla el valor añadido del componente espacial GP manteniendo fijo el marco de herradura, al tiempo que elimina la estructura espacial. EGP aísla el valor de la aproximación HSGP y la reducción de grupos, manteniendo la estructura espacial pero eliminándolos. CT proporciona un límite inferior ingenuo. Todos los métodos se aplican a las mismas matrices de expresión escaladas. Los métodos basados en MCMC (GP-GHS y SHS) utilizan 2.000 iteraciones con 500 iteraciones de "burn-in" y un adelgazamiento de 5.

Métricas de evaluación

Sea la matriz de adyacencia estimada y la matriz de adyacencia verdadera. La evaluación se realiza en la mitad superior de ambas matrices, tratando la recuperación de aristas como un problema de clasificación binaria sobre las posibles aristas para , 25. Calculamos las siguientes métricas.

Sean TP, FP, TN, FN los verdaderos positivos, falsos positivos, verdaderos negativos y falsos negativos, respectivamente. Las métricas son:

El coeficiente de correlación de Matthews (MCC) es particularmente informativo en casos de desequilibrio de clases ([Chicco y Jurman, 2020]), lo que ocurre en el entorno muy disperso donde las aristas verdaderas constituyen solo aproximadamente el 10% de todos los pares. En la misma línea, la puntuación F1 favorece los métodos que seleccionan más verdaderos positivos sin cometer errores de tipo I o tipo II. Todas las métricas se promedian en las 10 réplicas por nivel de dispersión, y presentamos las medias con las desviaciones estándar. También se registra el tiempo de ejecución real por conjunto de datos.

Resultados: n = 600, p = 15

Los resultados para la configuración n = 600, p = 15 se presentan en la Figura 2. Organizamos la discusión en torno a cinco hallazgos sustanciales que, en conjunto, caracterizan las propiedades operativas de GP-GHS y el competidor Exact GP recién añadido en relación con los métodos restantes.

Hallazgo 1: GP-GHS supera a todos los competidores en F1 y MCC en todo el rango de dispersión, y Exact GP solo es competitivo en gráficos densos.

GP-GHS alcanza puntuaciones F1 de aproximadamente 0,78, 0,88, 0,88 y 0,83 en las probabilidades de arista π ∈ {0,1, 0,3, 0,5, 0,7}, respectivamente, con valores de MCC correspondientes de 0,83, 0,80, 0,78 y 0,62. El competidor Exact GP muestra una trayectoria cualitativamente diferente: su F1 aumenta monótonamente desde cerca de cero en π = 0,1 hasta aproximadamente 0,83 en π = 0,7, igualando efectivamente a GP-GHS solo en la configuración más densa. En todos los entornos más dispersos, Exact GP tiene un rendimiento significativamente inferior a GP-GHS en F1 y MCC, a pesar de compartir el mismo prior espacial de Matérn. Todos los métodos restantes (GLasso, Nodewise Lasso, Standard Horseshoe, Correlation Threshold) se mantienen por debajo de F1 ≈ 0,20 en cada nivel de dispersión, lo que confirma que ni la verosimilitud penalizada ni la reducción bayesiana escalar son adecuadas para la recuperación de gráficos con estructura espacial.

El rendimiento de GP-GHS alcanza su punto máximo en el rango de disperso a moderado (π ∈ {0,3, 0,5}) y disminuye modestamente en el extremo denso, un patrón consistente con el comportamiento del parámetro global de herradura . Cuando el gráfico es moderadamente disperso, encuentra un punto operativo bien separado que discrimina entre los bloques de vecinos activos e inactivos. En π = 0,7, con aproximadamente 74 aristas verdaderas de 105, el prior difunde su masa de reducción en un gran número de bloques activos, lo que reduce la discriminación fina y produce la disminución observada en MCC de 0,83 a 0,62.

Hallazgo 2: Exact GP alcanza una alta sensibilidad en todo momento, pero sufre una grave inflación de falsos descubrimientos en gráficos dispersos.

El TPR de Exact GP es notablemente estable en los niveles de dispersión, permaneciendo cerca de 0,85 desde π = 0,1 hasta π = 0,7. Este es el TPR más alto de cualquier método en entornos moderados y densos, y coincide con GP-GHS en entornos muy dispersos y dispersos. Sin embargo, su FDR es catastróficamente alto en configuraciones dispersas: aproximadamente 0,90 en π = 0,1 y 0,70 en π = 0,3, lo que significa que la gran mayoría de las aristas declaradas son falsos positivos. El FDR de Exact GP mejora a medida que el gráfico se densifica, alcanzando aproximadamente 0,30 en π = 0,7, pero sigue siendo sustancialmente superior a GP-GHS en todos los entornos.

Este comportamiento refleja una propiedad estructural de la regla de selección de aristas de Exact GP. El criterio de votación por mayoría (declarar una arista activa si el intervalo de credibilidad puntual del 95% excluye cero en más del 50% de las ubicaciones) trata cada par espacial de forma independiente y no tiene una regularización global entre las aristas. En un gráfico disperso, el posterior marginal para cada par de nodos integra solo la señal residual local para ese par, sin aprovechar la información de que la mayoría de los pares deberían ser nulos. En consecuencia, los intervalos de credibilidad posteriores para los pares de ruido frecuente excluyen cero en una fracción no trivial de las ubicaciones debido al ruido espacialmente coherente, lo que produce una amplia inclusión de falsos positivos. GP-GHS evita esto al colocar un prior de herradura de grupo en toda la colección de bloques de vecinos de forma conjunta, de modo que el parámetro global de reducción calibre la dispersión a nivel de arista de una manera adaptativa a los datos que refleje la densidad general del gráfico. El contraste de FDR entre los dos métodos espaciales en π = 0,1 (aproximadamente 0,90 para Exact GP frente a 0,55 para GP-GHS) cuantifica directamente el valor de esta regularización conjunta de la dispersión.

Hallazgo 3: GP-GHS mantiene una alta sensibilidad y controla el FDR de forma progresiva en los niveles de dispersión.

GP-GHS alcanza un TPR de aproximadamente 0,85, 0,85, 0,57 y 0,48 en los cuatro niveles de dispersión, consistentemente el más alto entre los métodos que también mantienen un FDR razonable. El FDR de GP-GHS disminuye monótonamente de ≈ 0,55 en π = 0,1 a ≈ 0,05 en π = 0,7, lo que indica que la precisión mejora a medida que el gráfico se densifica y la jerarquía de herradura concentra su masa de forma más uniforme en las aristas activas. El FDR elevado en π = 0,1 es de esperar: en este nivel de dispersión, solo existen aproximadamente 11 aristas verdaderas, y el parámetro global de reducción debe suprimir simultáneamente aproximadamente 94 pares nulos al tiempo que protege un pequeño número de señales verdaderas, un régimen que inevitablemente sacrifica algo de precisión en aras de la exhaustividad. Las grandes desviaciones estándar en el FDR de GP-GHS en π = 0,1 (que abarcan aproximadamente de 0,30 a 0,80 en las réplicas) reflejan una verdadera variabilidad de réplica a réplica impulsada por el pequeño número de señales verdaderas en lugar de la inestabilidad del estimador.

Hallazgo 4: SHS y GLasso confirman que la estructura espacial y la reducción de grupos son conjuntamente necesarias.

Standard Horseshoe (SHS) comparte el marco bayesiano nodo a nodo con GP-GHS, pero descarta el prior de grupo y el prior espectral de GP en favor de una reducción escalar independiente en cada coeficiente de base. Su F1 se mantiene por debajo de 0,05 y su TPR cerca de cero en todos los niveles de dispersión, un rendimiento indistinguible del de GLasso y Correlation Threshold. El fracaso casi completo de SHS es informativo: establece que la estructura bayesiana nodo a nodo por sí sola no confiere ninguna ventaja sobre los métodos de penalización estándar cuando se ignora la estructura espacial. Las contribuciones fundamentales de GP-GHS son la jerarquía de herradura de grupo, que agrega evidencia de todos los coeficientes de base espectral para que cada inclusión de arista sea una decisión conjunta, y el prior espectral de GP, que codifica la expectativa de que los campos de interacción verdaderos sean espacialmente suaves y, por lo tanto, generen una acumulación coherente de señales en las frecuencias de base. Sin ambos componentes, el prior de herradura no tiene un mecanismo para distinguir una arista débil pero espacialmente estructurada del ruido independiente.

GLasso tiene un rendimiento uniformemente deficiente (F1 < 0,10 en todo momento) por una razón diferente: estima una única matriz de precisión estacionaria a partir de 600 observaciones espacialmente variables, confundiendo el promedio espacial de cada campo de interacción con su asociación marginal global. La penalización opera entonces sobre esta señal promediada en lugar de la interacción resuelta espacialmente, seleccionando aristas basándose en la estructura de correlación marginal que no está relacionada con la red de dependencia condicional verdadera. Esta no es una limitación del gráfico lasso como herramienta general, sino más bien una demostración de que los métodos diseñados para observaciones intercambiables no se pueden aplicar a datos espacialmente heterogéneos sin una modificación fundamental.

Hallazgo 5: GP-GHS alcanza su rendimiento a un costo computacional que sigue siendo manejable, mientras que Exact GP es el competidor más costoso.

GP-GHS tiene un tiempo de reloj promedio de aproximadamente 2 a 3 segundos por conjunto de datos en estos experimentos, lo que refleja la paralelización de 11 núcleos de p = 15 regresiones nodo a nodo con funciones de base por dimensión (componentes espectrales) y NMC = 2.000 iteraciones. La aproximación HSGP reduce el costo por iteración de – que sería necesario para la inferencia exacta de GP — a productos de matrices, lo que hace que el muestreador sea factible en n = 600. El tiempo de ejecución es constante en los niveles de dispersión, como se esperaba, ya que el cuello de botella computacional es la dimensión del sistema de base espectral en lugar de la densidad del gráfico verdadero.

Exact GP es el método más costoso, con aproximadamente 10 a 20 segundos por conjunto de datos, aproximadamente 5 a 8 veces más lento que GP-GHS, a pesar de ejecutarse de forma secuencial en lugar de en paralelo. Su costo crece como debido a la factorización de Cholesky requerida por cada par nodo-vecino, y el tiempo de ejecución casi constante en los niveles de dispersión refleja el hecho de que todos los pares de nodos p(p - 1) se evalúan independientemente de la arista verdadera. Nodewise Lasso y SHS ocupan un rango intermedio de aproximadamente 2 a 3 segundos, mientras que GLasso es el método más rápido, con aproximadamente 0,1 segundos, debido a sus actualizaciones de descenso de coordenadas de forma cerrada. Correlation Threshold es prácticamente instantáneo. La ventaja computacional de GP-GHS sobre Exact GP, combinada con su control de FDR sustancialmente mejor, establece que la aproximación HSGP y el prior de herradura de grupo juntos superan a la línea de base de Exact GP tanto estadística como computacionalmente.

Resultados: n = 600, p = 25

Los resultados para la configuración de mayor dimensión (n = 600, p = 25, 300 posibles aristas) se presentan en la Figura 3. La configuración de p = 25 expande la dimensión de la regresión nodo a nodo de q = 14 a q = 24 predictores por respuesta y aumenta el conjunto de candidatos de aristas de 105 a 300 pares. La pregunta central es si la ventaja de rendimiento de GP-GHS documentada en la Sección 3.4 sobrevive a la transición a un problema más difícil, y si el comportamiento relativo de Exact GP cambia cualitativamente con la dimensión.

Hallazgo 1: GP-GHS mantiene un rendimiento dominante en F1 y MCC en p = 25, y Exact GP solo es competitivo en gráficos densos.

GP-GHS lidera todos los métodos en F1 y MCC en todo el rango de dispersión, y la magnitud de su ventaja sobre los competidores no espaciales es comparable a la observada en p = 15. La curva F1 alcanza su punto máximo en una densidad dispersa y se mantiene fuerte en entornos moderados y densos, con solo una ligera disminución en relación con p = 15, que es pequeña dado el espacio de aristas candidatas sustancialmente más difícil. Ningún competidor supera una quinta parte del valor de F1 de GP-GHS en ningún nivel de dispersión. Exact GP muestra nuevamente una trayectoria F1 en aumento que converge con GP-GHS solo en el extremo denso, donde la prevalencia de aristas verdaderas es lo suficientemente alta como para que cualquier método de alta sensibilidad logre una precisión razonable por defecto. En todos los entornos más dispersos, Exact GP se queda significativamente por detrás de GP-GHS en F1 y MCC, a pesar de mantener una sensibilidad casi perfecta, lo que confirma que un TPR alto por sí solo no es suficiente para una recuperación competitiva de gráficos cuando el conjunto de aristas nulas es grande.

Hallazgo 2: Exact GP logra un TPR uniformemente alto, pero su inflación de FDR empeora con la dimensión.

El TPR de Exact GP es notablemente estable en los niveles de dispersión, permaneciendo cerca de 0,85 desde π = 0,1 hasta π = 0,7. Este es el TPR más alto de cualquier método en entornos moderados y densos, y coincide con GP-GHS en entornos muy dispersos y dispersos. Sin embargo, su FDR es catastróficamente alto en las configuraciones dispersas: aproximadamente 0,90 en π = 0,1 y 0,70 en π = 0,3, lo que significa que la gran mayoría de las aristas declaradas son falsos positivos. El FDR de Exact GP mejora a medida que el gráfico se densifica, alcanzando aproximadamente 0,30 en π = 0,7, pero sigue siendo sustancialmente superior a GP-GHS en todos los entornos.

Este comportamiento refleja una propiedad estructural de la regla de selección de aristas de Exact GP. El criterio de votación por mayoría (declarar una arista activa si el intervalo de credibilidad puntual del 95% excluye cero en más del 50% de las ubicaciones) trata cada par espacial de forma independiente y no tiene una regularización global entre las aristas. En un gráfico disperso, el posterior marginal para cada par de nodos integra solo la señal residual local para ese par, sin aprovechar la información de que la mayoría de los pares deberían ser nulos. En consecuencia, los intervalos de credibilidad posteriores para los pares de ruido frecuente excluyen cero en una fracción no trivial de las ubicaciones debido al ruido espacialmente coherente, lo que produce una amplia inclusión de falsos positivos. GP-GHS evita esto al colocar un prior de herradura de grupo en toda la colección de bloques de vecinos de forma conjunta, de modo que el parámetro global de reducción calibre la dispersión a nivel de arista de una manera adaptativa a los datos que refleje la densidad general del gráfico. El contraste de FDR entre los dos métodos espaciales en π = 0,1 (aproximadamente 0,90 para Exact GP frente a 0,55 para GP-GHS) cuantifica directamente el valor de esta regularización conjunta de la dispersión.

Hallazgo 3: GP-GHS mantiene una alta sensibilidad y controla el FDR de forma progresiva en los niveles de dispersión.

GP-GHS alcanza un TPR de aproximadamente 0,85, 0,85, 0,57 y 0,48 en los cuatro niveles de dispersión, consistentemente el más alto entre los métodos que también mantienen un FDR razonable. El FDR de GP-GHS disminuye monótonamente de ≈ 0,55 en π = 0,1 a ≈ 0,05 en π = 0,7, lo que indica que la precisión mejora a medida que el gráfico se densifica y la jerarquía de herradura concentra su masa de forma más uniforme en las aristas activas. El FDR elevado en π = 0,1 es de esperar: en este nivel de dispersión, solo existen aproximadamente 11 aristas verdaderas, y el parámetro global de reducción debe suprimir simultáneamente aproximadamente 94 pares nulos al tiempo que protege un pequeño número de señales verdaderas, un régimen que inevitablemente sacrifica algo de precisión en aras de la exhaustividad. Las grandes desviaciones estándar en el FDR de GP-GHS en π = 0,1 (que abarcan aproximadamente de 0,30 a 0,80 en las réplicas) reflejan una verdadera variabilidad de réplica a réplica impulsada por el pequeño número de señales verdaderas en lugar de la inestabilidad del estimador.

Hallazgo 4: SHS y GLasso confirman que la estructura espacial y la reducción de grupos son conjuntamente necesarias.

Standard Horseshoe (SHS) comparte el marco bayesiano nodo a nodo con GP-GHS, pero descarta el prior de grupo y el prior espectral de GP en favor de una reducción escalar independiente en cada coeficiente de base. Su F1 se mantiene por debajo de 0,05 y su TPR cerca de cero en todos los niveles de dispersión, un rendimiento indistinguible del de GLasso y Correlation Threshold. El fracaso casi completo de SHS es informativo: establece que la estructura bayesiana nodo a nodo por sí sola no confiere ninguna ventaja sobre los métodos de penalización estándar cuando se ignora la estructura espacial. Las contribuciones fundamentales de GP-GHS son la jerarquía de herradura de grupo, que agrega evidencia de todos los coeficientes de base espectral para que cada inclusión de arista sea una decisión conjunta, y el prior espectral de GP, que codifica la expectativa de que los campos de interacción verdaderos sean espacialmente suaves y, por lo tanto, generen una acumulación coherente de señales en las frecuencias de base. Sin ambos componentes, el prior de herradura no tiene un mecanismo para distinguir una arista débil pero espacialmente estructurada del ruido independiente.

GLasso tiene un rendimiento uniformemente deficiente (F1 < 0,10 en todo momento) por una razón diferente: estima una única matriz de precisión estacionaria a partir de 600 observaciones espacialmente variables, confundiendo el promedio espacial de cada campo de interacción con su asociación marginal global. La penalización opera entonces sobre esta señal promediada en lugar de la interacción resuelta espacialmente, seleccionando aristas basándose en la estructura de correlación marginal que no está relacionada con la red de dependencia condicional verdadera. Esta no es una limitación del gráfico lasso como herramienta general, sino más bien una demostración de que los métodos diseñados para observaciones intercambiables no se pueden aplicar a datos espacialmente heterogéneos sin una modificación fundamental.

Hallazgo 5: GP-GHS alcanza su rendimiento a un costo computacional que sigue siendo manejable, mientras que Exact GP es el competidor más costoso.

GP-GHS tiene un tiempo de reloj promedio de aproximadamente 2 a 3 segundos por conjunto de datos en estos experimentos, lo que refleja la paralelización de 11 núcleos de p = 15 regresiones nodo a nodo con funciones de base por dimensión (componentes espectrales) y NMC = 2.000 iteraciones. La aproximación HSGP reduce el costo por iteración de – que sería necesario para la inferencia exacta de GP — a productos de matrices, lo que hace que el muestreador sea factible en n = 600. El tiempo de ejecución es constante en los niveles de dispersión, como se esperaba, ya que el cuello de botella computacional es la dimensión del sistema de base espectral en lugar de la densidad del gráfico verdadero.

Exact GP es el método más costoso, con aproximadamente 10 a 20 segundos por conjunto de datos, aproximadamente 5 a 8 veces más lento que GP-GHS, a pesar de ejecutarse de forma secuencial en lugar de en paralelo. Su costo crece como debido a la factorización de Cholesky requerida por cada par nodo-vecino, y el tiempo de ejecución casi constante en los niveles de dispersión refleja el hecho de que todos los pares de nodos p(p - 1) se evalúan independientemente de la

Exact GP mantiene una TPR cercana a la constante y uniformemente alta en todos los niveles de dispersión en p = 25, reproduciendo exactamente su comportamiento en p = 15. Sin embargo, su FDR es notablemente peor en p = 25 que en p = 15 en los entornos dispersos y moderados, siendo la diferencia más visible en la densidad moderada, donde el FDR se duplica aproximadamente en relación con el problema más pequeño. Este empeoramiento es fácilmente comprensible: la regla de selección por mayoría evalúa cada par de nodos de forma independiente, sin ninguna calibración global, por lo que, a medida que el conjunto de aristas nulas crece con p, las inclusiones falsas esperadas se acumulan proporcionalmente. GP-GHS evita esto mediante el parámetro global de herradura, que se ajusta a la baja a medida que el espacio nulo se expande y suprime inclusiones falsas adicionales de forma adaptativa a los datos. Por lo tanto, el valor de esta regularización conjunta de dispersión se escala con p, y su ventaja sobre Exact GP se vuelve más pronunciada a mayor dimensión.### Hallazgo 3: El control del FDR de GP-GHS mejora con la densidad del grafo, pero el régimen muy disperso sigue siendo el más difícil.

El perfil del FDR de GP-GHS sigue la misma trayectoria cualitativa que en p = 15: alto en el extremo muy disperso y disminuyendo monótonamente a medida que el grafo se densifica, alcanzando una buena precisión en entornos densos. El FDR en el extremo muy disperso es ligeramente más alto que en p = 15, lo que refleja el mayor conjunto de aristas nulas en p = 25, pero la disminución en los niveles de dispersión es más pronunciada, de modo que el FDR en los entornos moderados y densos es comparable entre los dos tamaños de problema. La TPR sigue el patrón de disminución esperado con la densidad del grafo, por las mismas razones descritas en la Sección 3.4: a medida que el grafo se densifica, el parámetro de contracción global se estabiliza en un valor más alto, lo que difunde el poder discriminatorio en el límite. La gran variabilidad de réplica a réplica en todas las métricas de GP-GHS en el extremo muy disperso confirma que este régimen es genuinamente difícil independientemente de la dimensión, impulsado por la pequeña proporción de señales verdaderas en relación con las aristas candidatas, en lugar de la inestabilidad del muestreador.### Hallazgo 4: El argumento de ablación para la estructura de grupo se fortalece en p = 25.

Standard Horseshoe sigue fallando en la recuperación de aristas en todos los niveles de dispersión en p = 25, seleccionando un número moderado de aristas que no tienen una relación coherente con el grafo verdadero. El mecanismo se vuelve más claro a medida que p crece: cada regresión a nivel de nodo ahora implica un número sustancialmente mayor de coeficientes de base escalares, y la contracción escalar independiente no tiene un mecanismo para acoplar los coeficientes que pertenecen al mismo bloque de vecinos en una decisión de selección coherente. El prior de grupo en GP-GHS realiza una única selección binaria por bloque de vecinos, independientemente del tamaño del bloque, por lo que el número de decisiones de selección a nivel de grupo por nodo es q = 24 en lugar de qm2 = 384. A medida que q aumenta, el número absoluto de decisiones escalares espurias en SHS crece, mientras que la estructura de grupo de GP-GHS permanece anclada en q decisiones por nodo, lo que hace que la ventaja de la contracción de grupo sea más visible a mayor dimensión. GLasso y Nodewise Lasso se mantienen cerca de sus valores mínimos en p = 15, lo que confirma que ninguno de los métodos de verosimilitud penalizada tiene un mecanismo para explotar la estructura de interacción espacialmente heterogénea, independientemente del tamaño del problema.### Hallazgo 5: La escala del tiempo de ejecución en p = 25 es la principal limitación práctica, siendo Exact GP ahora el método más costoso.

El tiempo de ejecución de GP-GHS aumenta sustancialmente de p = 15 a p = 25, lo que es consistente con la escala de las operaciones de matriz dominantes dentro del muestreador de Gibbs. Exact GP es ahora el método más costoso en general, creciendo más rápido que GP-GHS porque su costo se escala como — la factorización de Cholesky en n = 600 se repite para cada par de nodos ordenados, y el número de pares casi se triplica de p = 15 a p = 25. Nodewise Lasso y SHS también aumentan notablemente, lo que es consistente con sus costos de regresión. GLasso y Correlation Threshold siguen siendo insignificantes. El tiempo de ejecución es constante en los niveles de dispersión para todos los métodos, lo que confirma que el cuello de botella es la dimensión del problema en lugar de la densidad del grafo. El tiempo de ejecución de GP-GHS en p = 25 con paralelización de 11 núcleos es operativamente aceptable para estudios de imagen de tejidos a escala moderada; para conjuntos de datos más grandes, distribuir las regresiones a nivel de nodo en un clúster de computación o reducir la dimensión de la base HSGP de a son las vías naturales para acelerar aún más el proceso.### Comparación multidimensional

En conjunto, los experimentos de p = 15 y p = 25 establecen varias conclusiones sobre el comportamiento de la escala de GP-GHS y sus competidores. El dominio cualitativo de GP-GHS se conserva en ambos tamaños de problema, y las pérdidas absolutas de rendimiento del aumento de p son pequeñas en relación con el crecimiento del espacio de aristas candidatas, lo que sugiere que el modelo estadístico está bien especificado para grafos espacialmente estructurados en este rango de dimensiones. El costo principal del aumento de p es computacional en lugar de estadístico, y el cuello de botella práctico para aplicaciones más grandes es el tiempo de ejecución en lugar de la calidad inferencial.

El comportamiento de Exact GP en ambos entornos revela una propiedad estructural de la selección independiente por pares: se conserva una TPR alta con la dimensión, pero el control del FDR se degrada a medida que crece el conjunto de aristas nulas, porque no existe un mecanismo global para recalibrar el umbral de selección. GP-GHS resuelve esto mediante la regularización conjunta de la herradura de grupo, y la brecha entre los dos métodos espaciales en el FDR se amplía con p. Su convergencia en F1 en grafos densos es un artefacto específico del régimen de alta prevalencia de aristas en lugar de evidencia de equivalencia estadística.

El equilibrio entre FDR y TPR en el extremo muy disperso es esencialmente estable en ambos entornos para GP-GHS, lo que refleja una propiedad estructural del prior de herradura cuando la fracción de señal verdadera es pequeña. La calibración adaptativa del parámetro de contracción global o el umbral basado en datos del criterio de inclusión de aristas son direcciones naturales para mejorar la precisión en este régimen sin sacrificar la alta sensibilidad que distingue a GP-GHS de todos los competidores.

Mecanismo de generación de datos

Evaluamos el método propuesto en comparación con cuatro competidores a través de un estudio de simulación estructurada diseñado para reflejar las características clave de los datos de imagen de tejidos multiplexados. Consideramos un conjunto fijo de n = 600 ubicaciones espaciales extraídas uniformemente del cuadrado unitario [0, 10]2, con p = 15, 25 tipos de células. El grafo de interacción célula-célula verdadero se genera utilizando el modelo de apego preferencial de Barabási-Albert ([Barabási y Albert, 1999]), que produce una topología sin escala con un pequeño número de tipos de células centrales que tienen muchas conexiones y muchos tipos de células que tienen pocas conexiones. Esta arquitectura está biológicamente motivada: en los microambientes tumorales, un puñado de tipos de células centrales, como los macrófagos y las células tumorales, tienden a coordinar las interacciones con muchas otras, mientras que las poblaciones más periféricas, como las células plasmáticas o los mastocitos, participan en menos interacciones ([Chen y Mellman, 2017]). Examinamos cuatro niveles de dispersión al apuntar a probabilidades de aristas de 0.7, 0.5, 0.3 y 0.1, correspondientes a grafos verdaderos densos, moderados, dispersos y muy dispersos.

Condicionalmente a la matriz de adyacencia verdadera , los datos se generan de acuerdo con la ecuación estructural

donde y los coeficientes espacialmente variables verdaderos se extraen de un proceso gaussiano de media cero con kernel de Matérn 3/2, longitud de escala apropiada y desviación estándar marginal . Esta elección de en relación con el dominio normalizado produce patrones espaciales que varían en una fracción moderada de la extensión del tejido, lo que es consistente con la escala de los patrones de infiltración inmune observados en los datos de imagen multiplexada. Cada réplica de simulación extrae nuevas realizaciones de procesos gaussianos para todas las aristas verdaderas y nuevo ruido, lo que garantiza que los resultados reflejen la variabilidad en ambas configuraciones espaciales y realizaciones de ruido. Se generan diez réplicas independientes por nivel de dispersión, lo que da un total de 40 conjuntos de datos simulados.

Métodos de la competencia

Comparamos GP-GHS con cinco métodos de la competencia que abarcan la verosimilitud penalizada, la contracción escalar bayesiana, la regresión exacta de GP y el umbral de correlación ingenuo.

Graphical Lasso (GLasso).

Graphical lasso ([Friedman et al., 2008]) estima una matriz de precisión dispersa resolviendo , donde es la matriz de covarianza de muestra de () y es un parámetro de penalización seleccionado mediante BIC en una cuadrícula de 15 valores. GLasso asume la estacionariedad e ignora por completo las coordenadas espaciales, produciendo una única matriz de precisión global.### Nodewise Lasso (NL).

Siguiendo a [Meinshausen y Bühlmann (2006)], ajustamos una regresión penalizada - para cada nodo frente a todos los nodos restantes, con la penalización seleccionada mediante validación cruzada de cinco pliegues. El grafo no dirigido se recupera mediante la regla AND. Al igual que GLasso, NL trata todos los coeficientes de regresión como escalares y no tiene en cuenta la estructura espacial.### Standard Horseshoe (SHS).

Este método aplica el prior de herradura escalar de [Carvalho et al. (2010)] a cada regresión a nivel de nodo, asignando un parámetro de contracción local independiente a cada predictor escalar en lugar de uno por bloque de vecinos. Debido a que la contracción opera a nivel de coeficiente en lugar de a nivel de grupo, SHS no tiene un mecanismo para tomar decisiones coherentes de aristas de tipo "todo o nada": los coeficientes de base individuales dentro de un bloque de vecinos pueden contraerse selectivamente, lo que da como resultado estimaciones de campo de interacción que no son uniformemente cero ni uniformemente diferentes de cero. La inclusión de aristas se evalúa mediante el criterio de pseudo-p-valor aplicado a cada coeficiente escalar. SHS conserva el marco bayesiano a nivel de nodo de GP-GHS, pero descarta la estructura de grupo y el prior espectral de GP, aislando así su contribución conjunta al rendimiento.### Pairwise GP Regression (Exact GP –EGP).

Este competidor ajusta directamente una regresión independiente del proceso gaussiano para cada par de nodos ordenados , utilizando un kernel de Matérn-3/2 con una longitud de escala que coincide con la auto-seleccionada de HSGP. Para cada nodo de respuesta , el residuo parcial con respecto al vecino se forma mediante la regresión de todos los demás vecinos mediante OLS, y el residuo resultante se trata como una observación ruidosa de . Denominamos a este competidor "pairwise" porque no puede ajustar simultáneamente los priors de GP para todos los vecinos de un nodo determinado. Hacerlo requeriría invertir una matriz de covarianza de GP conjunta de dimensión , lo que es computacionalmente prohibitivo en los tamaños de problema considerados. El enfoque pairwise utiliza en su lugar residuos OLS de un solo paso para eliminar los vecinos restantes, reemplazando un ajuste de GP conjunto con una secuencia de ajustes de GP marginales a un costo sustancialmente menor. La inclusión de aristas se determina mediante votación por mayoría: una arista se declara activa si el intervalo de credibilidad posterior puntual del 95% excluye cero en más del 50% de las ubicaciones espaciales, y el grafo no dirigido final se simetriza mediante la regla AND. EGP es la línea de base exacta natural para GP-GHS: comparte el mismo prior espacial y la lógica de selección de aristas, pero renuncia a la aproximación HSGP, la contracción de la herradura de grupo y la estimación conjunta de varios vecinos, a un costo computacional sustancialmente mayor.### Correlation Threshold (CT).

Una línea de base ingenua que declara una arista siempre que en la matriz de correlación de muestra. CT no hace suposiciones estructurales y sirve como un límite de rendimiento inferior.

GLasso y NL representan métodos estándar de verosimilitud penalizada para la selección de modelos gráficos. SHS aísla el valor añadido del componente espacial de GP manteniendo fijo el marco de la herradura mientras elimina la estructura espacial. EGP aísla el valor de la aproximación HSGP y la contracción de grupo al mantener la estructura espacial mientras las elimina. CT proporciona un límite inferior ingenuo. Todos los métodos se aplican a las mismas matrices de expresión escaladas. Los métodos basados en MCMC (GP-GHS y SHS) utilizan 2000 iteraciones con 500 iteraciones de "burn-in" y un adelgazamiento de 5.

Graphical Lasso (GLasso).

El lasso gráfico ([Friedman et al., 2008]) estima una matriz de precisión dispersa resolviendo , donde es la matriz de covarianza muestral de () y es un parámetro de penalización seleccionado mediante BIC en una cuadrícula de 15 valores. GLasso asume estacionariedad e ignora por completo las coordenadas espaciales, produciendo una única matriz de precisión global.

Lasso por nodos (NL).

Siguiendo a [Meinshausen y Bühlmann (2006)], ajustamos una regresión -penalizada para cada nodo contra todos los nodos restantes, con la penalización seleccionada mediante validación cruzada de cinco pliegues. El grafo no dirigido se recupera mediante la regla AND. Al igual que GLasso, NL trata todos los coeficientes de regresión como escalares y no tiene en cuenta la estructura espacial.

Horseshoe estándar (SHS).

Este método aplica la distribución a priori de horseshoe escalar de [Carvalho et al. (2010)] a cada regresión por nodos, asignando un parámetro de reducción local independiente a cada predictor escalar en lugar de uno por bloque de vecinos. Dado que la reducción opera a nivel de coeficiente en lugar de a nivel de grupo, SHS no tiene un mecanismo para tomar decisiones coherentes de inclusión o exclusión de aristas: los coeficientes de base individuales dentro de un bloque de vecinos pueden reducirse selectivamente, lo que da como resultado estimaciones del campo de interacción que no son uniformemente cero ni uniformemente distintos de cero. La inclusión de aristas se evalúa mediante el criterio del pseudo-p-valor aplicado a cada coeficiente escalar. SHS conserva el marco bayesiano por nodos de GP-GHS, pero descarta tanto la estructura de grupo como la distribución a priori espectral de GP, aislando así su contribución conjunta al rendimiento.

Regresión GP por pares (GP exacta – EGP).

Este competidor ajusta directamente una regresión de proceso gaussiano independiente para cada par de nodos ordenado , utilizando un kernel Matérn-3/2 con una longitud de escala que coincide con la auto-seleccionada de HSGP. Para cada nodo de respuesta , el residuo parcial con respecto al vecino se forma mediante la eliminación de todos los demás vecinos a través de OLS, y el residuo resultante se trata como una observación ruidosa de . Denominamos a este competidor "por pares" porque no puede ajustar simultáneamente las distribuciones a priori de GP para todos los vecinos de un nodo dado. Hacerlo requeriría la inversión de una matriz de covarianza GP conjunta de dimensión , lo cual es computacionalmente prohibitivo para los tamaños de problema considerados. El enfoque por pares utiliza en su lugar residuos OLS de un solo paso para eliminar los vecinos restantes, reemplazando un ajuste GP conjunto con una secuencia de ajustes GP marginales a un costo sustancialmente menor. La inclusión de aristas se determina mediante votación por mayoría: se declara que una arista es activa si el intervalo creíble posterior del 95% excluye cero en más del 50% de las ubicaciones espaciales, y el grafo no dirigido final se simetriza mediante la regla AND. EGP es la línea de base exacta natural para GP-GHS: comparte la misma distribución a priori espacial y la lógica de selección de aristas, pero omite la aproximación HSGP, la reducción de horseshoe de grupo y la estimación conjunta de múltiples vecinos, a un costo computacional sustancialmente mayor.

Umbral de correlación (CT).

Una línea de base ingenua que declara una arista siempre que en la matriz de correlación muestral. CT no hace suposiciones estructurales y sirve como un límite de rendimiento inferior.

GLasso y NL representan métodos estándar de verosimilitud penalizada para la selección de modelos gráficos. SHS aísla el valor añadido del componente espacial GP manteniendo fija la estructura de horseshoe mientras elimina la estructura espacial. EGP aísla el valor de la aproximación HSGP y la reducción de grupo manteniendo la estructura espacial mientras las elimina. CT proporciona un límite inferior ingenuo. Todos los métodos se aplican a las mismas matrices de expresión escaladas. Los métodos basados en MCMC (GP-GHS y SHS) utilizan 2.000 iteraciones con 500 iteraciones de "burn-in" y un adelgazamiento de 5.

Métricas de evaluación

Sea denota la matriz de adyacencia estimada y la matriz de adyacencia verdadera. La evaluación se realiza en la mitad superior de ambas matrices, tratando la recuperación de aristas como un problema de clasificación binaria sobre las posibles aristas para , 25. Calculamos las siguientes métricas.

Sea TP, FP, TN, FN denotan los verdaderos positivos, falsos positivos, verdaderos negativos y falsos negativos respectivamente. Las métricas son:

El coeficiente de correlación de Matthews (MCC) es particularmente informativo en caso de desequilibrio de clases ([Chicco y Jurman, 2020]), lo que ocurre en el entorno muy disperso donde las aristas verdaderas constituyen solo aproximadamente el 10% de todos los pares. En la misma línea, la puntuación F1 favorece los métodos que seleccionan más verdaderos positivos sin cometer errores de tipo I o tipo II. Todas las métricas se promedian en las 10 réplicas por nivel de dispersión, y presentamos las medias con las desviaciones estándar. También se registra el tiempo de ejecución real por conjunto de datos.

Resultados: n = 600, p = 15

Los resultados para la configuración n = 600, p = 15 se presentan en la Figura 2. Organizamos la discusión en torno a cinco hallazgos sustantivos que, en conjunto, caracterizan las propiedades operativas de GP-GHS y el competidor Exact GP recién añadido en relación con los métodos restantes.

Hallazgo 1: GP-GHS supera a todos los competidores en F1 y MCC en todo el rango de dispersión, con Exact GP competitivo solo en grafos densos.

GP-GHS alcanza puntuaciones F1 de aproximadamente 0,78, 0,88, 0,88 y 0,83 en las probabilidades de aristas π ∈ {0,1, 0,3, 0,5, 0,7} respectivamente, con valores MCC correspondientes de 0,83, 0,80, 0,78 y 0,62. El competidor Exact GP muestra una trayectoria cualitativamente diferente: su F1 aumenta monótonamente desde cerca de cero en π = 0,1 hasta aproximadamente 0,83 en π = 0,7, igualando efectivamente a GP-GHS solo en la configuración más densa. En todos los entornos más dispersos, Exact GP tiene un rendimiento significativamente inferior a GP-GHS en F1 y MCC, a pesar de compartir la misma distribución a priori espacial de Matérn. Todos los métodos restantes (GLasso, Lasso por nodos, Horseshoe estándar, Umbral de correlación) se mantienen por debajo de F1 ≈ 0,20 en cada nivel de dispersión, lo que confirma que ni la verosimilitud penalizada ni la reducción bayesiana escalar son adecuadas para la recuperación de grafos con estructura espacial.

El rendimiento de GP-GHS alcanza su punto máximo en el rango disperso a moderado (π ∈ {0,3, 0,5}) y disminuye modestamente en el extremo denso, un patrón consistente con el comportamiento del parámetro global de horseshoe . Cuando el grafo es moderadamente disperso, encuentra un punto de operación bien separado que discrimina entre los bloques de vecinos activos e inactivos. En π = 0,7, con aproximadamente 74 aristas verdaderas de 105, la distribución a priori difunde su masa de reducción en un gran número de bloques activos, lo que reduce la discriminación fina y produce la disminución observada en MCC de 0,83 a 0,62.### Hallazgo 2: Exact GP logra una alta sensibilidad en todo momento, pero sufre una grave inflación de falsos descubrimientos en grafos dispersos.

El TPR de Exact GP es notablemente estable en los niveles de dispersión, permaneciendo cerca de 0,85 desde π = 0,1 hasta π = 0,7. Este es el TPR más alto de cualquier método en entornos moderados y densos y coincide con GP-GHS en entornos muy dispersos y dispersos. Sin embargo, su FDR es catastróficamente alto en configuraciones dispersas: aproximadamente 0,90 en π = 0,1 y 0,70 en π = 0,3, lo que significa que la gran mayoría de las aristas declaradas son falsos positivos. El FDR de Exact GP mejora a medida que el grafo se densifica, alcanzando aproximadamente 0,30 en π = 0,7, pero sigue siendo sustancialmente superior a GP-GHS en todos los entornos.

Este comportamiento refleja una propiedad estructural de la regla de selección de aristas de Exact GP. El criterio de votación por mayoría (declarar una arista activa si el intervalo creíble posterior del 95% excluye cero en más del 50% de las ubicaciones) trata cada par espacial de forma independiente y no tiene una regularización global entre las aristas. En un grafo disperso, el posterior marginal para cada par de nodos integra solo la señal residual local para ese par, sin aprovechar la información de que la mayoría de los pares deberían ser nulos. En consecuencia, los intervalos creíbles posteriores para los pares que solo contienen ruido a menudo excluyen cero en una fracción no trivial de las ubicaciones debido al ruido espacialmente coherente, lo que produce una amplia inclusión de falsos positivos. GP-GHS evita esto al colocar una distribución a priori de horseshoe de grupo en toda la colección de bloques de vecinos de forma conjunta, de modo que el parámetro global de reducción calibra la dispersión a nivel de arista de una manera adaptativa a los datos que refleja la densidad general del grafo. El contraste de FDR entre los dos métodos espaciales en π = 0,1 (aproximadamente 0,90 para Exact GP frente a 0,55 para GP-GHS) cuantifica directamente el valor de esta regularización conjunta de la dispersión.### Hallazgo 3: GP-GHS mantiene una alta sensibilidad y controla el FDR de forma progresiva en los niveles de dispersión.

GP-GHS alcanza un TPR de aproximadamente 0,85, 0,85, 0,57 y 0,48 en los cuatro niveles de dispersión, consistentemente el más alto entre los métodos que también mantienen un FDR razonable. El FDR de GP-GHS disminuye monótonamente de ≈ 0,55 en π = 0,1 a ≈ 0,05 en π = 0,7, lo que indica que la precisión mejora a medida que el grafo se densifica y la jerarquía de horseshoe concentra su masa de forma más uniforme en las aristas activas. El FDR elevado en π = 0,1 es de esperar: en este nivel de dispersión, solo existen aproximadamente 11 aristas verdaderas, y el parámetro global de reducción debe suprimir simultáneamente aproximadamente 94 pares nulos al tiempo que protege un pequeño número de señales verdaderas, un régimen que inevitablemente intercambia algo de precisión por la sensibilidad. Las grandes desviaciones estándar en el FDR de GP-GHS en π = 0,1 (que abarcan aproximadamente de 0,30 a 0,80 en las réplicas) reflejan una verdadera variabilidad de réplica a réplica impulsada por el pequeño número de señales verdaderas en lugar de la inestabilidad del estimador.### Hallazgo 4: SHS y GLasso confirman que la estructura espacial y la reducción de grupo son conjuntamente necesarias.

Standard Horseshoe (SHS) comparte el marco bayesiano por nodos con GP-GHS, pero descarta la distribución a priori de grupo y la distribución a priori espectral de GP a favor de una reducción escalar independiente en cada coeficiente de base. Su F1 se mantiene por debajo de 0,05 y su TPR cerca de cero en todos los niveles de dispersión, un rendimiento indistinguible de GLasso y Umbral de correlación. El fracaso casi completo de SHS es informativo: establece que la estructura bayesiana por nodos por sí sola no confiere ninguna ventaja sobre los métodos de verosimilitud penalizada estándar cuando no se tiene en cuenta la estructura espacial. Las contribuciones críticas de GP-GHS son la jerarquía de horseshoe de grupo, que agrega evidencia de todos los coeficientes de base espectrales para que la inclusión de cada arista sea una decisión conjunta, y la distribución a priori espectral de GP, que codifica la expectativa de que los campos de interacción verdaderos sean espacialmente suaves y, por lo tanto, generen una acumulación coherente de señales en las frecuencias de base. Sin ambos componentes, la distribución a priori de horseshoe no tiene un mecanismo para distinguir una arista débil pero espacialmente estructurada de fluctuaciones de ruido independientes.

GLasso tiene un rendimiento uniformemente deficiente (F1 < 0,10 en todo momento) por una razón diferente: estima una única matriz de precisión estacionaria a partir de 600 observaciones espacialmente variables, confundiendo el promedio espacial de cada campo de interacción con su asociación marginal global. La penalización opera entonces sobre esta señal promediada en lugar de sobre la interacción resuelta espacialmente, seleccionando aristas en función de la estructura de correlación marginal que no está relacionada con la red de dependencia condicional verdadera. Esta no es una limitación del lasso gráfico como herramienta general, sino más bien una demostración de que los métodos diseñados para observaciones intercambiables no se pueden aplicar a datos espacialmente heterogéneos sin una modificación fundamental.### Hallazgo 5: GP-GHS logra su rendimiento a un costo computacional que sigue siendo manejable, mientras que Exact GP es el competidor más costoso.

GP-GHS tiene un tiempo de reloj promedio de aproximadamente 2-3 segundos por conjunto de datos en estos experimentos, lo que refleja la paralelización de 11 núcleos de p = 15 regresiones por nodos con funciones de base por dimensión (componentes espectrales) y NMC = 2.000 iteraciones. La aproximación HSGP reduce el costo por iteración de – que sería necesario para la inferencia GP exacta — a productos de matrices, lo que hace que el muestreador sea factible en n = 600. El tiempo de ejecución es plano en los niveles de dispersión, como se esperaba, ya que el cuello de botella computacional es la dimensión del sistema de base espectral en lugar de la densidad del grafo verdadero.

Exact GP es el método más costoso, con un tiempo de aproximadamente 10 a 20 segundos por conjunto de datos, lo que supone unas 5 a 8 veces más lento que GP-GHS, a pesar de que se ejecuta de forma secuencial en lugar de en paralelo. Su coste aumenta a medida que lo hace debido a la factorización de Cholesky requerida para cada par nodo-vecino, y el tiempo de ejecución casi constante en los diferentes niveles de dispersión refleja el hecho de que se evalúan todos los pares de nodos p(p - 1), independientemente del grafo real. Nodewise Lasso y SHS ocupan un rango intermedio de aproximadamente 2 a 3 segundos, mientras que GLasso es el método más rápido, con aproximadamente 0,1 segundos, debido a sus actualizaciones de descenso de coordenadas en forma cerrada. Correlation Threshold es prácticamente instantáneo. La ventaja computacional de GP-GHS sobre Exact GP, combinada con su control del FDR sustancialmente mejor, establece que la aproximación HSGP y la distribución previa de grupo horseshoe juntos superan a la distribución previa exacta GP tanto a nivel estadístico como computacional.

Hallazgo 1: GP-GHS supera a todos los competidores en F1 y MCC en todo el rango de dispersión, y Exact GP solo es competitivo en grafos densos.

GP-GHS alcanza puntuaciones F1 de aproximadamente 0,78, 0,88, 0,88 y 0,83 para probabilidades de arista π ∈ {0,1, 0,3, 0,5, 0,7}, respectivamente, con valores MCC correspondientes de 0,83, 0,80, 0,78 y 0,62. El competidor Exact GP muestra una trayectoria cualitativamente diferente: su F1 aumenta de forma monótona desde cerca de cero en π = 0,1 hasta aproximadamente 0,83 en π = 0,7, igualando efectivamente a GP-GHS solo en la configuración más densa. En todos los entornos más dispersos, Exact GP tiene un rendimiento sustancialmente inferior a GP-GHS en F1 y MCC, a pesar de compartir la misma distribución espacial de Matérn. Todos los métodos restantes (GLasso, Nodewise Lasso, Standard Horseshoe, Correlation Threshold) se mantienen por debajo de F1 ≈ 0,20 en todos los niveles de dispersión, lo que confirma que ni la verosimilitud penalizada ni la reducción bayesiana escalar son adecuados para la recuperación de grafos con estructura espacial.

El rendimiento de GP-GHS alcanza su punto máximo en el rango de dispersión moderada a baja (π ∈ {0,3, 0,5}) y disminuye modestamente en el extremo denso, un patrón consistente con el comportamiento del parámetro global horseshoe. Cuando el grafo es moderadamente disperso, encuentra un punto de operación bien separado que discrimina entre los bloques de vecinos activos e inactivos. En π = 0,7, con aproximadamente 74 aristas verdaderas de 105, la distribución previa difunde su masa de reducción en un gran número de bloques activos, lo que reduce la discriminación fina y produce la disminución observada en MCC de 0,83 a 0,62.

Hallazgo 2: Exact GP alcanza una alta sensibilidad en todo momento, pero sufre una grave inflación de falsos descubrimientos en grafos dispersos.

La TPR de Exact GP es notablemente estable en los diferentes niveles de dispersión, manteniéndose cerca de 0,85 desde π = 0,1 hasta π = 0,7. Esta es la TPR más alta de todos los métodos en entornos moderados y densos, y coincide con GP-GHS en entornos muy dispersos y dispersos. Sin embargo, su FDR es catastróficamente alto en las configuraciones dispersas: aproximadamente 0,90 en π = 0,1 y 0,70 en π = 0,3, lo que significa que la gran mayoría de las aristas declaradas son falsos positivos. El FDR de Exact GP mejora a medida que el grafo se densifica, alcanzando aproximadamente 0,30 en π = 0,7, pero sigue siendo sustancialmente superior a GP-GHS en todos los entornos.

Este comportamiento refleja una propiedad estructural de la regla de selección de aristas de Exact GP. El criterio de votación por mayoría (declarar una arista activa si el intervalo de credibilidad puntual del 95% excluye el cero en más del 50% de las ubicaciones) trata cada par espacial de forma independiente y no tiene una regularización global entre las aristas. En un grafo disperso, la distribución posterior marginal para cada par de nodos integra solo la señal residual local para ese par, sin aprovechar la información de que la mayoría de los pares deberían ser nulos. En consecuencia, los intervalos de credibilidad posteriores para los pares con solo ruido a menudo excluyen el cero en una fracción no trivial de las ubicaciones debido al ruido espacialmente coherente, lo que produce una amplia inclusión de falsos positivos. GP-GHS evita esto al colocar una distribución previa de grupo horseshoe en toda la colección de bloques de vecinos de forma conjunta, de modo que el parámetro global de reducción calibra la dispersión a nivel de arista de forma adaptativa a los datos, lo que refleja la densidad general del grafo. El contraste de FDR entre los dos métodos espaciales en π = 0,1 (aproximadamente 0,90 para Exact GP frente a 0,55 para GP-GHS) cuantifica directamente el valor de esta regularización conjunta de la dispersión.

Hallazgo 3: GP-GHS mantiene una alta sensibilidad y controla el FDR de forma progresiva en los diferentes niveles de dispersión.

GP-GHS alcanza una TPR de aproximadamente 0,85, 0,85, 0,57 y 0,48 en los cuatro niveles de dispersión, siendo consistentemente la más alta entre los métodos que también mantienen un FDR razonable. El FDR de GP-GHS disminuye de forma monótona desde ≈ 0,55 en π = 0,1 hasta ≈ 0,05 en π = 0,7, lo que indica que la precisión mejora a medida que el grafo se densifica y la jerarquía horseshoe concentra su masa de forma más uniforme en las aristas activas. El FDR elevado en π = 0,1 es de esperar: en este nivel de dispersión, solo existen aproximadamente 11 aristas verdaderas, y el parámetro global de reducción debe suprimir simultáneamente aproximadamente 94 pares nulos al tiempo que protege un pequeño número de señales verdaderas, un régimen que inevitablemente sacrifica algo de precisión en aras de la exhaustividad. Las grandes desviaciones estándar en el FDR de GP-GHS en π = 0,1 (que abarcan aproximadamente de 0,30 a 0,80 en las diferentes réplicas) reflejan una verdadera variabilidad de una réplica a otra, impulsada por el pequeño número de señales verdaderas y no por la inestabilidad del estimador.

Hallazgo 4: SHS y GLasso confirman que la estructura espacial y la reducción de grupos son conjuntamente necesarias.

Standard Horseshoe (SHS) comparte el marco bayesiano a nivel de nodo con GP-GHS, pero descarta la distribución previa de grupo y la distribución previa espectral GP en favor de una reducción escalar independiente en cada coeficiente de base. Su F1 se mantiene por debajo de 0,05 y su TPR cerca de cero en todos los niveles de dispersión, un rendimiento indistinguible del de GLasso y Correlation Threshold. El fracaso casi completo de SHS es revelador: establece que la estructura bayesiana a nivel de nodo por sí sola no confiere ninguna ventaja sobre los métodos de penalización estándar cuando se ignora la estructura espacial. Las contribuciones fundamentales de GP-GHS son la jerarquía de grupo horseshoe, que agrega pruebas en todos los coeficientes de base espectral para que cada inclusión de arista sea una decisión conjunta, y la distribución previa espectral GP, que codifica la expectativa de que los campos de interacción verdaderos sean espacialmente suaves y, por lo tanto, generen una acumulación coherente de la señal en las diferentes frecuencias de base. Sin ambos componentes, la distribución previa horseshoe no tiene un mecanismo para distinguir una arista débil pero espacialmente estructurada del ruido independiente.

GLasso tiene un rendimiento uniformemente deficiente (F1 < 0,10 en todo momento) por una razón diferente: estima una única matriz de precisión estacionaria a partir de 600 observaciones espacialmente variables, confundiendo el promedio espacial de cada campo de interacción con su asociación marginal global. La penalización opera entonces sobre esta señal promediada en lugar de sobre la interacción resuelta espacialmente, seleccionando aristas en función de la estructura de correlación marginal que no está relacionada con la red de dependencia condicional verdadera. Esta no es una limitación del graphical lasso como herramienta general, sino más bien una demostración de que los métodos diseñados para observaciones intercambiables no se pueden aplicar a datos espacialmente heterogéneos sin una modificación fundamental.

Hallazgo 5: GP-GHS alcanza su rendimiento a un coste computacional que sigue siendo manejable, mientras que Exact GP es el competidor más costoso.

GP-GHS tiene un tiempo medio de ejecución de aproximadamente 2 a 3 segundos por conjunto de datos en estos experimentos, lo que refleja la paralelización de 11 núcleos de p = 15 regresiones a nivel de nodo con funciones de base por dimensión (componentes espectrales) y NMC = 2000 iteraciones. La aproximación HSGP reduce el coste por iteración de – que sería necesario para la inferencia exacta de GP — a productos de matrices, lo que hace que el muestreador sea factible en n = 600. El tiempo de ejecución es constante en los diferentes niveles de dispersión, como se esperaba, ya que el cuello de botella computacional es la dimensión del sistema de base espectral y no la densidad del grafo real.

Exact GP es el método más costoso, con un tiempo de aproximadamente 10 a 20 segundos por conjunto de datos, lo que supone unas 5 a 8 veces más lento que GP-GHS, a pesar de que se ejecuta de forma secuencial en lugar de en paralelo. Su coste aumenta a medida que lo hace debido a la factorización de Cholesky requerida para cada par nodo-vecino, y el tiempo de ejecución casi constante en los diferentes niveles de dispersión refleja el hecho de que se evalúan todos los pares de nodos p(p - 1), independientemente del grafo real. Nodewise Lasso y SHS ocupan un rango intermedio de aproximadamente 2 a 3 segundos, mientras que GLasso es el método más rápido, con aproximadamente 0,1 segundos, debido a sus actualizaciones de descenso de coordenadas en forma cerrada. Correlation Threshold es prácticamente instantáneo. La ventaja computacional de GP-GHS sobre Exact GP, combinada con su control del FDR sustancialmente mejor, establece que la aproximación HSGP y la distribución previa de grupo horseshoe juntos superan a la distribución previa exacta GP tanto a nivel estadístico como computacional.

Resultados: n = 600, p = 25

Los resultados para la configuración de mayor dimensión (n = 600, p = 25, 300 posibles aristas) se presentan en la Figura 3. El ajuste de p = 25 amplía la dimensión de la regresión a nivel de nodo de q = 14 a q = 24 predictores por respuesta y aumenta el conjunto de pares candidatos de 105 a 300. La pregunta clave es si la ventaja de rendimiento de GP-GHS documentada en la Sección 3.4 se mantiene en la transición a un problema más difícil, y si el comportamiento relativo de Exact GP cambia cualitativamente con la dimensión.

Hallazgo 1: GP-GHS mantiene un rendimiento dominante en F1 y MCC en p = 25, y Exact GP vuelve a ser competitivo solo en grafos densos.

GP-GHS lidera todos los métodos en F1 y MCC en todo el rango de dispersión, y la magnitud de su ventaja sobre los competidores no espaciales es comparable a la observada en p = 15. La curva F1 alcanza su punto máximo en una densidad dispersa y se mantiene fuerte en entornos moderados y densos, con solo una ligera disminución en relación con p = 15, lo que es pequeño dado el espacio de aristas candidatas sustancialmente más difícil. Ningún competidor supera a GP-GHS en más de una quinta parte del valor de F1 en ningún nivel de dispersión. Exact GP muestra de nuevo una trayectoria F1 en ascenso que converge con GP-GHS solo en el extremo denso, donde la prevalencia de aristas verdaderas es lo suficientemente alta como para que cualquier método de alta sensibilidad alcance una precisión razonable por defecto. En todos los entornos más dispersos, Exact GP se queda muy por detrás de GP-GHS tanto en F1 como en MCC, a pesar de mantener una sensibilidad casi perfecta, lo que confirma que una TPR alta por sí sola no es suficiente para una recuperación competitiva del grafo cuando el conjunto de aristas nulas es grande.

Hallazgo 2: Exact GP alcanza una TPR uniformemente alta, pero su inflación de FDR empeora con la dimensión.

Exact GP mantiene una TPR casi constante y uniformemente alta en todos los niveles de dispersión en p = 25, reproduciendo exactamente su comportamiento en p = 15. Sin embargo, su FDR es notablemente peor en p = 25 que en p = 15 en los entornos dispersos y moderados, y la diferencia es más visible en la densidad moderada, donde el FDR se duplica aproximadamente en relación con el problema más pequeño. Este empeoramiento es mecánicamente transparente: la regla de selección de aristas de votación por mayoría evalúa cada par de nodos de forma independiente, sin ninguna calibración global, por lo que, a medida que el conjunto de aristas nulas crece con p, las inclusiones falsas esperadas se acumulan proporcionalmente. GP-GHS evita esto mediante el parámetro global horseshoe, que se ajusta a la baja a medida que el espacio nulo se expande y suprime inclusiones falsas adicionales de forma adaptativa a los datos. El valor de esta regularización conjunta de la dispersión, por lo tanto, se escala con p, y su ventaja sobre Exact GP se vuelve más pronunciada en una dimensión más alta.

Hallazgo 3: El control del FDR de GP-GHS mejora con la densidad del grafo, pero el régimen muy disperso sigue siendo el entorno más difícil.

GP-GHS alcanza una TPR de aproximadamente 0,85, 0,85, 0,57 y 0,48 en los cuatro niveles de dispersión, siendo consistentemente la más alta entre los métodos que también mantienen un FDR razonable. El FDR de GP-GHS disminuye de forma monótona desde ≈ 0,55 en π = 0,1 hasta ≈ 0,05 en π = 0,7, lo que indica que la precisión mejora a medida que el grafo se densifica y la jerarquía horseshoe concentra su masa de forma más uniforme en las aristas activas. El FDR elevado en π = 0,1 es de esperar: en este nivel de dispersión, solo existen aproximadamente 11 aristas verdaderas, y el parámetro global de reducción debe suprimir simultáneamente aproximadamente 94 pares nulos al tiempo que protege un pequeño número de señales verdaderas, un régimen que inevitablemente sacrifica algo de precisión en aras de la exhaustividad. Las grandes desviaciones estándar en el FDR de GP-GHS en π = 0,1 (que abarcan aproximadamente de 0,30 a 0,80 en las diferentes réplicas) reflejan una verdadera variabilidad de una réplica a otra, impulsada por el pequeño número de señales verdaderas y no por la inestabilidad del estimador.

Hallazgo 4: SHS y GLasso confirman que la estructura espacial y la reducción de grupos son conjuntamente necesarias.

Standard Horseshoe (SHS) comparte el marco bayesiano a nivel de nodo con GP-GHS, pero descarta la distribución previa de grupo y la distribución previa espectral GP en favor de una reducción escalar independiente en cada coeficiente de base. Su F1 se mantiene por debajo de 0,05 y su TPR cerca de cero en todos los niveles de dispersión, un rendimiento indistinguible del de GLasso y Correlation Threshold. El fracaso casi completo de SHS es revelador: establece que la estructura bayesiana a nivel de nodo por sí sola no confiere ninguna ventaja sobre los métodos de penalización estándar cuando se ignora la estructura espacial. Las contribuciones fundamentales de GP-GHS son la jerarquía de grupo horseshoe, que agrega pruebas en todos los coeficientes de base espectral para que cada inclusión de arista sea una decisión conjunta, y la distribución previa espectral GP, que codifica la expectativa de que los campos de interacción verdaderos sean espacialmente suaves y, por lo tanto, generen una acumulación coherente de la señal en las diferentes frecuencias de base. Sin ambos componentes, la distribución previa horseshoe no tiene un mecanismo para distinguir una arista débil pero espacialmente estructurada del ruido independiente.

GLasso tiene un rendimiento uniformemente deficiente (F1 < 0,10 en todo momento) por una razón diferente: estima una única matriz de precisión estacionaria a partir de 600 observaciones espacialmente variables, confundiendo el promedio espacial de cada campo de interacción con su asociación marginal global. La penalización opera entonces sobre esta señal promediada en lugar de sobre la interacción resuelta espacialmente, seleccionando aristas en función de la estructura de correlación marginal que no está relacionada con la red de dependencia condicional verdadera. Esta no es una limitación del graphical lasso como herramienta general, sino más bien una demostración de que los métodos diseñados para observaciones intercambiables no se pueden aplicar a datos espacialmente heterogéneos sin una modificación fundamental.

Hallazgo 5: GP-GHS alcanza su rendimiento a un coste computacional que sigue siendo manejable, mientras que Exact GP es el competidor más costoso.

GP-GHS tiene un tiempo medio de ejecución de aproximadamente 2 a 3 segundos por conjunto de datos en estos experimentos, lo que refleja la paralelización de 11 núcleos de p = 15 regresiones a nivel de nodo con funciones de base por dimensión (componentes espectrales) y NMC = 2000 iteraciones. La aproximación HSGP reduce el coste por iteración de – que sería necesario para la inferencia exacta de GP — a productos de matrices, lo que hace que el muestreador sea factible en n = 600. El tiempo de ejecución es constante en los diferentes niveles de dispersión, como se esperaba, ya que el cuello de botella computacional es la dimensión del sistema de base espectral y no la densidad del grafo real.

Exact GP es el método más costoso, con un tiempo de aproximadamente 10 a 20 segundos por conjunto de datos, lo que supone unas 5 a 8 veces más lento que GP-GHS, a pesar de que se ejecuta de forma secuencial en lugar de en paralelo. Su coste aumenta a medida que lo hace debido a la factorización de Cholesky requerida para cada par nodo-vecino, y el tiempo de ejecución casi constante en los diferentes niveles de dispersión refleja el hecho de que se evalúan todos los pares de nodos p(p - 1), independientemente del grafo real. Nodewise Lasso y SHS ocupan un rango intermedio de aproximadamente 2 a 3 segundos, mientras que GLasso es el método más rápido, con aproximadamente 0,1 segundos, debido a sus actualizaciones de descenso de coordenadas en forma cerrada. Correlation Threshold es prácticamente instantáneo. La ventaja computacional de GP-GHS sobre Exact GP, combinada con su control del FDR sustancialmente mejor, establece que la aproximación HSGP y la distribución previa de grupo horseshoe juntos superan a la distribución previa exacta GP tanto a nivel estadístico como computacional.

El perfil FDR de GP-GHS sigue la misma trayectoria cualitativa que en p = 15: alto en entornos muy dispersos y disminuyendo monótonamente a medida que el grafo se densifica, alcanzando una buena precisión en entornos densos. El FDR en el extremo muy disperso es ligeramente superior al de p = 15, lo que refleja el mayor conjunto de aristas nulas en p = 25, pero la disminución a través de los niveles de dispersión es más pronunciada, de modo que el FDR en entornos moderados y densos es comparable entre los dos tamaños de problema. El TPR sigue el patrón de disminución esperado con la densidad del grafo, por las mismas razones descritas en la Sección 3.4: a medida que el grafo se densifica, el parámetro de contracción global se estabiliza en un valor mayor que difunde el poder discriminatorio en el margen. La gran variabilidad de un replicado a otro en todas las métricas de GP-GHS en entornos muy dispersos confirma que este régimen es genuinamente difícil independientemente de la dimensión, debido a la pequeña proporción de señales verdaderas en relación con las aristas candidatas, en lugar de la inestabilidad del muestreador.

Hallazgo 4: El argumento de ablación para la estructura de grupo se fortalece en p = 25.

Standard Horseshoe continúa fallando en la recuperación de aristas en todos los niveles de dispersión en p = 25, seleccionando un número moderado de aristas que no guardan una relación coherente con el grafo verdadero. El mecanismo se vuelve más claro a medida que p aumenta: cada regresión a nivel de nodo ahora implica sustancialmente más coeficientes de base escalares, y la contracción escalar independiente no tiene un mecanismo para acoplar los coeficientes que pertenecen al mismo bloque de vecinos en una decisión de selección coherente. El prior de grupo en GP-GHS realiza una única selección binaria por bloque de vecinos independientemente del tamaño del bloque, por lo que el número de decisiones de selección a nivel de grupo por nodo es q = 24 en lugar de qm2 = 384. A medida que q aumenta, el número absoluto de decisiones espurias a nivel escalar en SHS crece, mientras que la estructura de grupo de GP-GHS permanece anclada en q decisiones por nodo, lo que hace que la ventaja de la contracción de grupo sea más visible a mayor dimensión. GLasso y Nodewise Lasso se mantienen cerca de sus valores mínimos de rendimiento en p = 15, lo que confirma que ninguno de los dos métodos de máxima verosimilitud penalizada tiene un mecanismo para explotar una estructura de interacción espacialmente heterogénea independientemente del tamaño del problema.

Hallazgo 5: La escala del tiempo de ejecución en p = 25 es la principal limitación práctica, siendo Exact GP ahora el método más costoso.

El tiempo de ejecución de GP-GHS aumenta sustancialmente de p = 15 a p = 25, lo que es consistente con la escala de las operaciones de matriz dominantes dentro del muestreador de Gibbs. Exact GP es ahora el método más costoso en general, y su crecimiento es más rápido que el de GP-GHS porque su costo se escala como — la factorización de Cholesky en n = 600 se repite para cada par de nodos ordenados, y el número de pares casi se triplica de p = 15 a p = 25. Nodewise Lasso y SHS también aumentan notablemente, lo que es consistente con sus costos de regresión. GLasso y Correlation Threshold siguen siendo insignificantes. El tiempo de ejecución es constante en todos los niveles de dispersión para todos los métodos, lo que confirma que el cuello de botella es la dimensión del problema en lugar de la densidad del grafo. El tiempo de ejecución de GP-GHS en p = 25 con una paralelización de 11 núcleos es operativamente aceptable para estudios de imagen de tejidos a escala moderada; para conjuntos de datos más grandes, distribuir las regresiones a nivel de nodo en un clúster de computación o reducir la dimensión de la base HSGP de a son las vías naturales para acelerar aún más el proceso.

Hallazgo 1: GP-GHS mantiene un rendimiento dominante en F1 y MCC en p = 25, y Exact GP vuelve a ser competitivo solo en grafos densos.

GP-GHS lidera todos los métodos en F1 y MCC en todo el rango de dispersión, y la magnitud de su ventaja sobre los competidores no espaciales es comparable a la observada en p = 15. La curva de F1 alcanza su punto máximo en una densidad dispersa y se mantiene fuerte en entornos moderados y densos, con solo una ligera disminución en relación con p = 15, que es pequeña dado el espacio de aristas candidatas sustancialmente más difícil. Ningún competidor supera un quinto del valor de F1 de GP-GHS en ningún nivel de dispersión. Exact GP vuelve a mostrar una trayectoria de F1 en aumento que converge con GP-GHS solo en el extremo denso, donde la prevalencia de aristas verdaderas es lo suficientemente alta como para que cualquier método de alta sensibilidad logre una precisión razonable de forma predeterminada. En todos los entornos más dispersos, Exact GP supera sustancialmente a GP-GHS tanto en F1 como en MCC, a pesar de mantener una sensibilidad casi perfecta, lo que confirma que un TPR alto por sí solo no es suficiente para una recuperación de grafos competitiva cuando el conjunto de aristas nulas es grande.

Hallazgo 2: Exact GP logra un TPR uniformemente alto, pero su inflación de FDR empeora con la dimensión.

Exact GP mantiene un TPR casi constante y uniformemente alto en todos los niveles de dispersión en p = 25, reproduciendo exactamente su comportamiento en p = 15. Sin embargo, su FDR es notablemente peor en p = 25 que en p = 15 en los entornos dispersos y moderados, y la diferencia es más visible en una densidad moderada, donde el FDR se duplica aproximadamente en relación con el problema más pequeño. Este empeoramiento es mecánicamente transparente: la regla de selección por votación mayoritaria evalúa cada par de nodos de forma independiente sin ninguna calibración global, por lo que a medida que el conjunto de aristas nulas crece con p, las inclusiones falsas esperadas se acumulan proporcionalmente. GP-GHS evita esto a través del parámetro de contracción global del herradura, que se ajusta hacia abajo a medida que el espacio nulo se expande y suprime inclusiones falsas adicionales de forma adaptativa a los datos. El valor de esta regularización de la dispersión conjunta, por lo tanto, se escala con p, y su ventaja sobre Exact GP se vuelve más pronunciada a mayor dimensión.

Hallazgo 3: El control de FDR de GP-GHS mejora con la densidad del grafo, pero el régimen muy disperso sigue siendo el entorno más difícil.

El perfil FDR de GP-GHS sigue la misma trayectoria cualitativa que en p = 15: alto en entornos muy dispersos y disminuyendo monótonamente a medida que el grafo se densifica, alcanzando una buena precisión en entornos densos. El FDR en el extremo muy disperso es ligeramente superior al de p = 15, lo que refleja el mayor conjunto de aristas nulas en p = 25, pero la disminución a través de los niveles de dispersión es más pronunciada, de modo que el FDR en entornos moderados y densos es comparable entre los dos tamaños de problema. El TPR sigue el patrón de disminución esperado con la densidad del grafo, por las mismas razones descritas en la Sección 3.4: a medida que el grafo se densifica, el parámetro de contracción global se estabiliza en un valor mayor que difunde el poder discriminatorio en el margen. La gran variabilidad de un replicado a otro en todas las métricas de GP-GHS en entornos muy dispersos confirma que este régimen es genuinamente difícil independientemente de la dimensión, debido a la pequeña proporción de señales verdaderas en relación con las aristas candidatas, en lugar de la inestabilidad del muestreador.

Hallazgo 4: El argumento de ablación para la estructura de grupo se fortalece en p = 25.

Standard Horseshoe continúa fallando en la recuperación de aristas en todos los niveles de dispersión en p = 25, seleccionando un número moderado de aristas que no guardan una relación coherente con el grafo verdadero. El mecanismo se vuelve más claro a medida que p aumenta: cada regresión a nivel de nodo ahora implica sustancialmente más coeficientes de base escalares, y la contracción escalar independiente no tiene un mecanismo para acoplar los coeficientes que pertenecen al mismo bloque de vecinos en una decisión de selección coherente. El prior de grupo en GP-GHS realiza una única selección binaria por bloque de vecinos independientemente del tamaño del bloque, por lo que el número de decisiones de selección a nivel de grupo por nodo es q = 24 en lugar de qm2 = 384. A medida que q aumenta, el número absoluto de decisiones espurias a nivel escalar en SHS crece, mientras que la estructura de grupo de GP-GHS permanece anclada en q decisiones por nodo, lo que hace que la ventaja de la contracción de grupo sea más visible a mayor dimensión. GLasso y Nodewise Lasso se mantienen cerca de sus valores mínimos de rendimiento en p = 15, lo que confirma que ninguno de los dos métodos de máxima verosimilitud penalizada tiene un mecanismo para explotar una estructura de interacción espacialmente heterogénea independientemente del tamaño del problema.

Hallazgo 5: La escala del tiempo de ejecución en p = 25 es la principal limitación práctica, siendo Exact GP ahora el método más costoso.

El tiempo de ejecución de GP-GHS aumenta sustancialmente de p = 15 a p = 25, lo que es consistente con la escala de las operaciones de matriz dominantes dentro del muestreador de Gibbs. Exact GP es ahora el método más costoso en general, y su crecimiento es más rápido que el de GP-GHS porque su costo se escala como — la factorización de Cholesky en n = 600 se repite para cada par de nodos ordenados, y el número de pares casi se triplica de p = 15 a p = 25. Nodewise Lasso y SHS también aumentan notablemente, lo que es consistente con sus costos de regresión. GLasso y Correlation Threshold siguen siendo insignificantes. El tiempo de ejecución es constante en todos los niveles de dispersión para todos los métodos, lo que confirma que el cuello de botella es la dimensión del problema en lugar de la densidad del grafo. El tiempo de ejecución de GP-GHS en p = 25 con una paralelización de 11 núcleos es operativamente aceptable para estudios de imagen de tejidos a escala moderada; para conjuntos de datos más grandes, distribuir las regresiones a nivel de nodo en un clúster de computación o reducir la dimensión de la base HSGP de a son las vías naturales para acelerar aún más el proceso.

Comparación Interdimensional

En conjunto, los experimentos de p = 15 y p = 25 establecen varias conclusiones sobre el comportamiento de la escala de GP-GHS y sus competidores. El dominio cualitativo de GP-GHS se conserva en ambos tamaños de problema, y las pérdidas absolutas de rendimiento debido al aumento de p son pequeñas en relación con el crecimiento del espacio de aristas candidatas, lo que sugiere que el modelo estadístico está bien especificado para grafos espacialmente estructurados en este rango de dimensiones. El principal costo de aumentar p es computacional en lugar de estadístico, y el cuello de botella práctico para aplicaciones más grandes es el tiempo de ejecución en lugar de la calidad inferencial.

El comportamiento de Exact GP en ambos entornos revela una propiedad estructural de la selección independiente por pares: se conserva un TPR alto con la dimensión, pero el control de FDR se deteriora a medida que crece el conjunto de aristas nulas, porque no existe un mecanismo global para recalibrar el umbral de selección. GP-GHS resuelve esto a través de la regularización conjunta del herradura de grupo, y la brecha entre los dos métodos espaciales en FDR se amplía con p. Su convergencia en F1 en grafos densos es un artefacto específico del régimen de alta prevalencia de aristas en lugar de evidencia de equivalencia estadística.

El equilibrio entre FDR y TPR en el extremo muy disperso es esencialmente estable en ambos entornos para GP-GHS, lo que refleja una propiedad estructural del prior del herradura cuando la fracción de señales verdaderas es pequeña. La calibración adaptativa del parámetro de contracción global o el umbral basado en datos del criterio de inclusión de aristas son direcciones naturales para mejorar la precisión en este régimen sin sacrificar la alta sensibilidad que distingue a GP-GHS de todos los competidores.

Análisis de Datos: Redes de Interacción Celular-Celular Espacial en el Cáncer Colorrectal

Datos y Preprocesamiento

Aplicamos el marco de trabajo de GP-GHS al conjunto de datos CODEX disponible públicamente de [Schürch et al. (2020)], que comprende datos de imagen de una sola célula multiplexada de 35 pacientes con cáncer colorrectal en 140 imágenes de tejido que abarcan dos microambientes tumorales patológicamente distintos: reacción similar a la enfermedad de Crohn (CLR, n = 68 imágenes) e infiltrado inflamatorio difuso (DII, n = 72 imágenes). Siguiendo la fenotipificación previa proporcionada con el conjunto de datos, cada célula se asigna a uno de los K tipos de células discretos. Retenemos los tipos de células presentes en al menos el 0,5% de todas las células en toda la cohorte, lo que da como resultado K = 15 tipos: células B, células T CD4+, células T CD8+, macrófagos CD163+, macrófagos CD68+, células T reguladoras (Tregs), células T CD4+ de memoria, células inmunitarias genéricas, granulocitos, células plasmáticas, músculo liso, estroma, células tumorales, vasos sanguíneos y adipocitos. Cada paciente contribuye con exactamente cuatro imágenes, y los pacientes están anidados dentro del grupo de patología, con 17 pacientes asignados a CLR y 18 a DII.

Construcción de Características Espaciales mediante Estimación de la Densidad del Núcleo

El modelo GP-GHS requiere una matriz de características continua y referenciada espacialmente como entrada. Las coordenadas celulares sin procesar constituyen un proceso de puntos espaciales marcado con una dispersión significativa dentro de la imagen, lo que no es adecuado para la regresión espacial directa. Para convertir las observaciones de patrones de puntos en un formato adecuado para GP-GHS, estimamos una superficie de abundancia espacial suave para cada tipo de célula en cada imagen utilizando la estimación de densidad de kernel bidimensional (KDE) para capturar el gradiente de heterogeneidad de cada tipo de célula a través de las densidades marginales. Para cada imagen, definimos una cuadrícula regular de 30 × 30 sobre el cuadro delimitador de la imagen, expandida en un 5% en cada lado para mitigar los efectos de borde, lo que da como resultado 900 ubicaciones de cuadrícula por imagen. La KDE para cada tipo de célula se evalúa en esta cuadrícula común utilizando MASS::kde2d, con un ancho de banda seleccionado de forma independiente por cada eje espacial mediante la regla de referencia normal. Los tipos de células con menos de cinco células observadas en una imagen determinada se les asigna una densidad constante pequeña (10−6) para mantener la estabilidad numérica. La matriz de densidad resultante de 900 × K se transforma logarítmicamente después de escalarla por 104 para evitar el desbordamiento numérico, y cada columna se estandariza posteriormente a una media de cero y una varianza unitaria. Las coordenadas de la cuadrícula sirven directamente como ubicaciones espaciales para la construcción de la base HSGP. Con este fin, utilizamos estos perfiles KDE marginales para tipos de células individuales como expresiones celulares en nuestro modelo.

Inferencia de la red GP-GHS

Ajustamos el modelo GP-GHS de forma independiente a cada una de las 140 imágenes. Para cada imagen, se ejecutan en paralelo K = 15 regresiones a nivel de nodo, cada una de las cuales trata la superficie KDE de un tipo de célula como la variable dependiente y las K - 1 superficies restantes como predictores estructurados espacialmente. La aproximación HSGP utiliza funciones base por dimensión espacial ( funciones base en total) con un kernel de covarianza Matérn-3/2. El muestreador MCMC se ejecuta durante 3000 iteraciones con un período de "burn-in" de 1000 y un factor de adelgazamiento de 5, conservando 400 muestras posteriores por nodo. Los nodos cuya varianza de respuesta cae por debajo de 10−8 después de la estandarización se excluyen de los modelos de regresión y se les asignan sin bordes, ya que una superficie KDE espacialmente constante no proporciona información para la regresión. Se declara que un borde entre los tipos de células y es activo si la media posterior de la estadística de contracción , donde con valores cercanos a cero que indican una fuerte evidencia de interacción y valores cercanos a uno que indican una contracción casi completa hacia el valor nulo. La matriz de adyacencia se simetriza mediante una regla AND, que requiere una inclusión mutua de ambas regresiones direccionales, y la matriz se simetriza tomando el mínimo por elementos, lo que constituye la opción más conservadora al requerir que ambos nodos proporcionen evidencia de una interacción compartida. Todo el proceso se ejecutó hasta su finalización en las 140 imágenes, con un tiempo total de ejecución de aproximadamente 251 minutos en un servidor de múltiples núcleos.

Datos y preprocesamiento

Aplicamos el marco GP-GHS al conjunto de datos CODEX disponible públicamente de [Schürch et al. (2020)], que comprende datos de imagen de una sola célula multiplexados de 35 pacientes con cáncer colorrectal en 140 imágenes de tejido que abarcan dos microambientes tumorales patológicamente distintos: reacción similar a la enfermedad de Crohn (CLR, n = 68 imágenes) e infiltrado inflamatorio difuso (DII, n = 72 imágenes). Siguiendo la fenotipificación previa proporcionada con el conjunto de datos, cada célula se asigna a uno de los K tipos de células discretos. Conservamos los tipos de células presentes en al menos el 0,5% de todas las células en toda la cohorte, lo que da como resultado K = 15 tipos: células B, células T CD4+, células T CD8+, macrófagos CD163+, macrófagos CD68+, células T reguladoras (Tregs), células T CD4+ de memoria, células inmunitarias genéricas, granulocitos, células plasmáticas, músculo liso, estroma, células tumorales, vasos sanguíneos y adipocitos. Cada paciente contribuye con exactamente cuatro imágenes, y los pacientes están anidados dentro del grupo de patología, con 17 pacientes asignados a CLR y 18 a DII.

Construcción de características espaciales mediante la estimación de la densidad del kernel

El modelo GP-GHS requiere una matriz de características continua y referenciada espacialmente como entrada. Las coordenadas celulares sin procesar constituyen un proceso de puntos espaciales marcado con una dispersión significativa dentro de la imagen, lo que no es adecuado para la regresión espacial directa. Para convertir las observaciones de patrones de puntos en un formato adecuado para GP-GHS, estimamos una superficie de abundancia espacial suave para cada tipo de célula en cada imagen utilizando la estimación de densidad de kernel bidimensional (KDE) para capturar el gradiente de heterogeneidad de cada tipo de célula a través de las densidades marginales. Para cada imagen, definimos una cuadrícula regular de 30 × 30 sobre el cuadro delimitador de la imagen, expandida en un 5% en cada lado para mitigar los efectos de borde, lo que da como resultado 900 ubicaciones de cuadrícula por imagen. La KDE para cada tipo de célula se evalúa en esta cuadrícula común utilizando MASS::kde2d, con un ancho de banda seleccionado de forma independiente por cada eje espacial mediante la regla de referencia normal. Los tipos de células con menos de cinco células observadas en una imagen determinada se les asigna una densidad constante pequeña (10−6) para mantener la estabilidad numérica. La matriz de densidad resultante de 900 × K se transforma logarítmicamente después de escalarla por 104 para evitar el desbordamiento numérico, y cada columna se estandariza posteriormente a una media de cero y una varianza unitaria. Las coordenadas de la cuadrícula sirven directamente como ubicaciones espaciales para la construcción de la base HSGP. Con este fin, utilizamos estos perfiles KDE marginales para tipos de células individuales como expresiones celulares en nuestro modelo.

Inferencia de la red GP-GHS

Ajustamos el modelo GP-GHS de forma independiente a cada una de las 140 imágenes. Para cada imagen, se ejecutan en paralelo K = 15 regresiones a nivel de nodo, cada una de las cuales trata la superficie KDE de un tipo de célula como la variable dependiente y las K - 1 superficies restantes como predictores estructurados espacialmente. La aproximación HSGP utiliza funciones base por dimensión espacial ( funciones base en total) con un kernel de covarianza Matérn-3/2. El muestreador MCMC se ejecuta durante 3000 iteraciones con un período de "burn-in" de 1000 y un factor de adelgazamiento de 5, conservando 400 muestras posteriores por nodo. Los nodos cuya varianza de respuesta cae por debajo de 10−8 después de la estandarización se excluyen de los modelos de regresión y se les asignan sin bordes, ya que una superficie KDE espacialmente constante no proporciona información para la regresión. Se declara que un borde entre los tipos de células y es activo si la media posterior de la estadística de contracción , donde con valores cercanos a cero que indican una fuerte evidencia de interacción y valores cercanos a uno que indican una contracción casi completa hacia el valor nulo. La matriz de adyacencia se simetriza mediante una regla AND, que requiere una inclusión mutua de ambas regresiones direccionales, y la matriz se simetriza tomando el mínimo por elementos, lo que constituye la opción más conservadora al requerir que ambos nodos proporcionen evidencia de una interacción compartida. Todo el proceso se ejecutó hasta su finalización en las 140 imágenes, con un tiempo total de ejecución de aproximadamente 251 minutos en un servidor de múltiples núcleos.

Resultados

Gráficos de interacción espacial de imágenes de tejido de CRC

Aplicamos GP-GHS a un conjunto de datos de imágenes de tejido de cáncer colorrectal (CRC) que comprende 140 imágenes de 35 pacientes, con cuatro imágenes obtenidas por paciente. La clasificación de la patología asignó 17 pacientes (68 imágenes) al grupo CLR y 18 pacientes (72 imágenes) al grupo DII. Cada imagen contenía coordenadas espaciales y etiquetas de tipo de célula en 16 tipos de células: adipocitos, células B, macrófagos CD163+, células T CD4+, macrófagos CD68+, células T reguladoras (Tregs), células T CD4+ de memoria, células inmunitarias genéricas, granulocitos, células plasmáticas, músculo liso, estroma, células tumorales, vasos sanguíneos y adipocitos.

Para cada imagen, primero estimamos las superficies de abundancia de tipos de células utilizando la estimación de la densidad del kernel bidimensional en una cuadrícula espacial de 30 × 30, produciendo una matriz de características de 900 × 16 con una columna por tipo de célula. Las densidades referenciadas en la cuadrícula se transformaron logarítmicamente y se estandarizaron antes del modelado. GP-GHS se ajustó luego de forma independiente a cada imagen, lo que dio como resultado un gráfico de interacción espacial no dirigido con una matriz de adyacencia binaria y una matriz de contracción simétrica . Bajo la convención del prior de herradura, indica una interacción espacialmente covariante entre los tipos de células y , mientras que indica una contracción hacia el valor nulo. La simetrización de las estimaciones de borde direccionales siguió una regla AND para la matriz de adyacencia y una regla de mínimo para , esta última seleccionando la posterior direccional con mayor evidencia de interacción.

Pruebas diferenciales de interacciones espaciales entre CLR y DII

Para identificar los bordes con una fuerza de interacción diferencial entre los grupos de patología, probamos todos los posibles pares de tipos de células; después de restringirnos a los bordes observados en las imágenes, se conservaron 105 bordes para la prueba. Para cada borde , modelamos la estadística de contracción log-transformada a nivel de imagen utilizando un modelo lineal mixto (LMM):

donde es un intercepto aleatorio a nivel de paciente que tiene en cuenta la correlación entre las cuatro imágenes por paciente, y . El nivel de referencia fue CLR. El efecto fijo es la cantidad inferencial primaria: representa la diferencia de grupo en después de ajustar la correlación dentro del paciente. Las estimaciones kappa marginales específicas del grupo y se transformaron de nuevo para facilitar la interpretación. Una prueba de razón de verosimilitud que compara el modelo completo con un modelo nulo sin el término de grupo produjo valores p a nivel de borde. Se aplicó una corrección para pruebas múltiples utilizando el procedimiento de Benjamini-Hochberg a un umbral de tasa de descubrimiento falso (FDR) del 5%.

Interacciones espaciales diferencialmente activas

Trece de los 105 bordes probados alcanzaron significación estadística a FDR < 0,05 (Tabla 1). Los 13 bordes significativos mostraron una covariación espacial más fuerte en DII en relación con CLR, como lo indican los valores uniformemente negativos (Figura 5). Ningún borde mostró una interacción significativamente más débil en DII.

Las señales diferenciales más fuertes y más significativas se concentraron en torno a las Tregs. Diez de los 13 bordes significativos involucraron a las Tregs como uno de los socios, abarcando interacciones con macrófagos CD68+ (, FDR = 0,0021), células T CD4+ de memoria (, FDR = 0,0021), vasos sanguíneos (, FDR = 0,0022), células T CD8+ (, FDR = 0,0022), estroma (, FDR = 0,0022), granulocitos (, FDR = 0,0022), células tumorales (, FDR = 0,0022), macrófagos CD163+ (, FDR = 0,0033), células B (, FDR = 0,0057) y músculo liso (, FDR = 0,013). Los tres bordes significativos restantes involucraron a los macrófagos CD68+ que interactúan con las células T CD4+ de memoria (, FDR = 0,037), el estroma (, FDR = 0,037) y el músculo liso (, FDR = 0,047).

Los valores kappa marginales estimados por el LMM ilustran la magnitud de estas diferencias grupales en la escala kappa (Figura 6). En CLR, los bordes que involucran a las Tregs tenían valores estimados cercanos a 1,0, lo que indica que estas interacciones se redujeron en gran medida hacia el valor nulo, lo que es consistente con una localización espacial dispersa o poco frecuente. En DII, los mismos bordes tenían valores que oscilaban entre aproximadamente 0,56 y 0,90, lo que indica una menor contracción y, por lo tanto, una interacción espacial más fuerte. Por ejemplo, el borde de macrófagos CD68+ y Tregs tenía versus , y el borde de células T CD4+ de memoria y Tregs tenía versus .

El mapa de calor beta-hat (Figura 7) confirma la asimetría direccional en todos los 105 pares simultáneamente. Los mapas de calor de la prevalencia de bordes (Figura 4) corroboran estos hallazgos a nivel de la topología del gráfico. El grupo DII mostró una prevalencia de bordes visiblemente mayor para las interacciones centradas en las Tregs y para los bordes que conectan los macrófagos CD68+ con las poblaciones estromales e inmunitarias, lo que es consistente con los resultados de las pruebas diferenciales basadas en LMM.

Interpretación biológica

Los resultados de las pruebas diferenciales apuntan a una red inmunosupresora centrada en las células T reguladoras (Tregs) que se amplifica selectivamente en DII en relación con CLR. En el tejido DII, las Tregs exhiben una colocalización espacial significativamente más fuerte con prácticamente todos los principales compartimentos inmunitarios y estromales representados en los datos, incluidos los macrófagos CD68+, los macrófagos CD163+, las células T CD8+, las células T CD4+ de memoria, las células B, los granulocitos, las células tumorales, los vasos sanguíneos, el estroma y el músculo liso. Los valores kappa estimados mediante LMM confirman que estas diferencias reflejan cambios genuinos en la estructura de dependencia espacial y no simplemente un aumento en la prevalencia de los bordes: en CLR, estos bordes que involucran a las Tregs se reducen casi por completo hacia el valor nulo (), mientras que en DII, los mismos bordes presentan una reducción posterior sustancialmente menor (que oscila entre 0,56 y 0,90), lo que indica que la colocalización de las Tregs con el microambiente inmunitario más amplio es una característica espacial definitoria de los tumores DII.

Este patrón es consistente con la biología conocida del subtipo de infiltrado inflamatorio difuso. Los tumores DII se caracterizan por un infiltrado inmunitario espacialmente disperso y heterogéneo, en el que se cree que las Tregs se acumulan en múltiples compartimentos tisulares y suprimen las respuestas inmunitarias efectoras a través de mecanismos dependientes del contacto y mediados por citocinas ([Schürch et al., 2020]). El acoplamiento espacial concurrente de las Tregs con las poblaciones mieloides (macrófagos CD68+ y CD163+) y las poblaciones linfoides (células T CD8+, células T CD4+ de memoria, células B) sugiere una supresión inmunitaria coordinada que opera en múltiples ramas de la respuesta adaptativa e innata. La fuerte interacción Treg↔vasculatura en DII () también puede reflejar nichos perivasculares de Tregs, una disposición espacial que se ha relacionado con una alteración del tráfico de las células inmunitarias hacia el parénquima tumoral.

Los tres bordes significativos de macrófagos CD68+—con las células T CD4+ de memoria, el estroma y el músculo liso—sugieren un eje secundario de reorganización espacial en DII que es independiente de las Tregs. Los macrófagos CD68+ en DII parecen estar incrustados en un contexto estromal y de músculo liso más denso, lo que puede reflejar un microambiente enriquecido en mieloides, lo que es consistente con la polarización de los macrófagos asociados al tumor hacia un fenotipo inmunosupresor. La colocalización de los macrófagos CD68+ con las células T CD4+ de memoria en DII podría indicar interacciones de presentación de antígenos que se facilitan espacialmente en este subtipo, pero que se interrumpen o se difunden espacialmente en CLR.

En contraste, el tejido CLR muestra una estructura de interacción espacial globalmente dispersa en todos los 105 bordes probados, con valores kappa estimados mediante LMM cercanos a 1 para la mayoría de los pares de tipos de células. Esto es consistente con la arquitectura inflamatoria similar a la enfermedad de Crohn de CLR, en la que la infiltración inmunitaria se organiza en agregados linfoides discretos con una distribución espacial más focal que difusa ([Schürch et al., 2020]). La relativa ausencia de fuertes dependencias espaciales en CLR puede reflejar una respuesta inmunitaria más compartimentada, donde las poblaciones efectoras están segregadas en lugar de estar espacialmente intermezcladas con los elementos reguladores y estromales. Es importante destacar que ningún borde mostró una interacción significativamente más fuerte en CLR en relación con DII, lo que indica que la remodelación de la red entre los subtipos es direccionalmente asimétrica y refleja predominantemente una ganancia de interacción en DII en lugar de una reestructuración bidireccional de la red de colocalización espacial.

En conjunto, estos hallazgos demuestran que GP-GHS recupera diferencias interpretables y biológicamente coherentes en las redes de interacción espacial a partir de datos de imagen multiplexada. El método no solo identifica qué bordes difieren, sino que también cuantifica la dirección y la magnitud de esas diferencias en una escala inferencial bien definida, lo que permite realizar comparaciones entre subgrupos de pacientes definidos patológicamente a una resolución que los análisis estándar de coocurrencia o enriquecimiento de vecindad no proporcionan.

Grafos de interacción espacial a partir de imágenes de tejido de CRC

Aplicamos GP-GHS a un conjunto de datos de imagen de tejido de cáncer colorrectal (CRC) que comprende 140 imágenes de 35 pacientes, con cuatro imágenes obtenidas por paciente. La clasificación patológica asignó 17 pacientes (68 imágenes) al grupo CLR y 18 pacientes (72 imágenes) al grupo DII. Cada imagen contenía coordenadas espaciales y etiquetas de tipo de célula en 16 tipos de células: adipocitos, células B, macrófagos CD163+, células T CD4+, macrófagos CD68+, células T CD8+, células inmunitarias genéricas, granulocitos, células T CD4+ de memoria, células plasmáticas, músculo liso, estroma, Tregs, células tumorales y vasos sanguíneos.

Para cada imagen, primero estimamos las superficies de abundancia de tipos de células utilizando la estimación de densidad del núcleo bidimensional en una cuadrícula espacial de 30 × 30, lo que produce una matriz de características de 900 × 16 con una columna por tipo de célula. Las densidades referenciadas a la cuadrícula se transformaron logarítmicamente y se estandarizaron antes del modelado. A continuación, se ajustó GP-GHS de forma independiente a cada imagen, lo que dio como resultado un grafo de interacción espacial no dirigido con una matriz de adyacencia binaria y una matriz de coeficientes de reducción simétrica . Bajo la convención del prior de herradura, indica una interacción espacialmente covariante entre los tipos de células y , mientras que indica una reducción hacia el valor nulo. La simetrización de las estimaciones de los bordes direccionales siguió una regla AND para la matriz de adyacencia y una regla mínima para , esta última seleccionando la dirección posterior con mayor evidencia de interacción.

Pruebas diferenciales de las interacciones espaciales entre CLR y DII

Para identificar los bordes con una fuerza de interacción diferencial entre los grupos patológicos, probamos todos los posibles pares de tipos de células; después de restringirnos a los bordes observados en las imágenes, se retuvieron 105 bordes para la prueba. Para cada borde , modelamos el coeficiente de reducción transformado logarítmicamente a nivel de imagen utilizando un modelo lineal mixto (LMM):

donde es un intercepto aleatorio a nivel de paciente que tiene en cuenta la correlación entre las cuatro imágenes por paciente, y . El nivel de referencia fue CLR. El efecto fijo es la cantidad inferencial primaria: representa la diferencia de grupo en después de ajustar por la correlación intra-paciente. Los valores kappa marginales específicos del grupo y se transformaron de nuevo para facilitar su interpretación. Una prueba de razón de verosimilitud que comparó el modelo completo con un modelo nulo sin el término de grupo produjo valores p a nivel de borde. Se aplicó una corrección para pruebas múltiples utilizando el procedimiento de Benjamini-Hochberg a un umbral de tasa de descubrimiento falso (FDR) del 5 %.

Interacciones espaciales diferencialmente activas

Trece de los 105 bordes probados alcanzaron significación estadística a FDR < 0,05 (Tabla 1). Todos los 13 bordes significativos mostraron una covariación espacial más fuerte en DII en relación con CLR, como lo indican los valores uniformemente negativos (Figura 5). Ningún borde mostró una interacción significativamente más débil en DII.

Las señales diferenciales más fuertes y estadísticamente significativas se concentraron en torno a las Tregs. Diez de los 13 bordes significativos involucraron a las Tregs como uno de los socios, abarcando interacciones con los macrófagos CD68+ (, FDR = 0,0021), las células T CD4+ de memoria (, FDR = 0,0021), los vasos sanguíneos (, FDR = 0,0022), las células T CD8+ (, FDR = 0,0022), el estroma (, FDR = 0,0022), los granulocitos (, FDR = 0,0022), las células tumorales (, FDR = 0,0022), los macrófagos CD163+ (, FDR = 0,0033), las células B (, FDR = 0,0057) y el músculo liso (, FDR = 0,013). Los tres bordes significativos restantes involucraron a los macrófagos CD68+ que interactuaban con las células T CD4+ de memoria (, FDR = 0,037), el estroma (, FDR = 0,037) y el músculo liso (, FDR = 0,047).

Los valores kappa marginales estimados mediante LMM ilustran la magnitud de estas diferencias de grupo en la escala kappa (Figura 6). En CLR, los bordes que involucran a las Tregs tenían valores kappa estimados cercanos a 1,0, lo que indica que estas interacciones se redujeron en gran medida hacia el valor nulo, lo que es consistente con una colocalización dispersa o espacialmente difusa. En DII, los mismos bordes tenían valores que oscilaban entre aproximadamente 0,56 y 0,90, lo que indica una reducción sustancialmente menor y, por lo tanto, una interacción espacial más fuerte. Por ejemplo, el borde macrófagos CD68+–Tregs tenía un valor de versus , y el borde células T CD4+ de memoria–Tregs tenía un valor de versus .

El mapa de calor beta-hat (Figura 7) confirma la asimetría direccional en todos los 105 pares simultáneamente. Los mapas de calor de la prevalencia de los bordes (Figura 4) corroboran estos hallazgos a nivel de la topología del grafo. El grupo DII mostró una prevalencia de bordes visiblemente mayor para las interacciones centradas en las Tregs y para los bordes que conectan los macrófagos CD68+ con las poblaciones estromales e inmunitarias, lo que es consistente con los resultados de las pruebas diferenciales basadas en LMM.

Interpretación biológica

Los resultados de las pruebas diferenciales apuntan a una red inmunosupresora centrada en las Tregs que se amplifica selectivamente en DII en relación con CLR. En el tejido DII, las Tregs exhiben una colocalización espacial significativamente más fuerte con prácticamente todos los principales compartimentos inmunitarios y estromales representados en los datos, incluidos los macrófagos CD68+, los macrófagos CD163+, las células T CD8+, las células T CD4+ de memoria, las células B, los granulocitos, las células tumorales, los vasos sanguíneos, el estroma y el músculo liso. Los valores kappa estimados mediante LMM confirman que estas diferencias reflejan cambios genuinos en la estructura de dependencia espacial y no simplemente un aumento en la prevalencia de los bordes: en CLR, estos bordes que involucran a las Tregs se reducen casi por completo hacia el valor nulo (), mientras que en DII, los mismos bordes presentan una reducción posterior sustancialmente menor (que oscila entre 0,56 y 0,90), lo que indica que la colocalización de las Tregs con el microambiente inmunitario más amplio es una característica espacial definitoria de los tumores DII.

Este patrón es consistente con la biología conocida del subtipo de infiltrado inflamatorio difuso. Los tumores DII se caracterizan por un infiltrado inmunitario espacialmente disperso y heterogéneo, en el que se cree que las Tregs se acumulan en múltiples compartimentos tisulares y suprimen las respuestas inmunitarias efectoras a través de mecanismos dependientes del contacto y mediados por citocinas ([Schürch et al., 2020]). El acoplamiento espacial concurrente de las Tregs con las poblaciones mieloides (macrófagos CD68+ y CD163+) y las poblaciones linfoides (células T CD8+, células T CD4+ de memoria, células B) sugiere una supresión inmunitaria coordinada que opera en múltiples ramas de la respuesta adaptativa e innata. La fuerte interacción Treg↔vasculatura en DII () también puede reflejar nichos perivasculares de Tregs, una disposición espacial que se ha relacionado con una alteración del tráfico de las células inmunitarias hacia el parénquima tumoral.

Los tres bordes significativos de macrófagos CD68+—con las células T CD4+ de memoria, el estroma y el músculo liso—sugieren un eje secundario de reorganización espacial en DII que es independiente de las Tregs. Los macrófagos CD68+ en DII parecen estar incrustados en un contexto estromal y de músculo liso más denso, lo que puede reflejar un microambiente enriquecido en mieloides, lo que es consistente con la polarización de los macrófagos asociados al tumor hacia un fenotipo inmunosupresor. La colocalización de los macrófagos CD68+ con las células T CD4+ de memoria en DII podría indicar interacciones de presentación de antígenos que se facilitan espacialmente en este subtipo, pero que se interrumpen o se difunden espacialmente en CLR.

En contraste, el tejido CLR muestra una estructura de interacción espacial globalmente dispersa en todos los 105 bordes probados, con valores kappa estimados mediante LMM cercanos a 1 para la mayoría de los pares de tipos de células. Esto es consistente con la arquitectura inflamatoria similar a la enfermedad de Crohn de CLR, en la que la infiltración inmunitaria se organiza en agregados linfoides discretos con una distribución espacial más focal que difusa ([Schürch et al., 2020]). La relativa ausencia de fuertes dependencias espaciales en CLR puede reflejar una respuesta inmunitaria más compartimentada, donde las poblaciones efectoras están segregadas en lugar de estar espacialmente intermezcladas con los elementos reguladores y estromales. Es importante destacar que ningún borde mostró una interacción significativamente más fuerte en CLR en relación con DII, lo que indica que la remodelación de la red entre los subtipos es direccionalmente asimétrica y refleja predominantemente una ganancia de interacción en DII en lugar de una reestructuración bidireccional de la red de colocalización espacial.

En conjunto, estos hallazgos demuestran que GP-GHS recupera diferencias interpretables y biológicamente coherentes en las redes de interacción espacial a partir de datos de imagen multiplexada. El método no solo identifica qué aristas difieren, sino que cuantifica la dirección y la magnitud de esas diferencias en una escala inferencial bien definida, lo que permite realizar comparaciones entre subgrupos de pacientes definidos patológicamente a una resolución que los análisis estándar de coocurrencia o enriquecimiento de vecindad no proporcionan.

Interpretación biológica

Los resultados de las pruebas diferenciales señalan una red inmunosupresora centrada en las células T reguladoras (Treg), que se amplifica selectivamente en DII en relación con CLR. En el tejido DII, las células Treg exhiben una colocalización espacial significativamente más fuerte con casi todos los principales compartimentos inmunitarios y estromales representados en los datos, incluidos los macrófagos CD68+, los macrófagos CD163+, las células T CD8+, las células T CD4+ de memoria, las células B, los granulocitos, las células tumorales, los vasos sanguíneos, el estroma y el músculo liso. Los valores kappa estimados por LMM confirman que estas diferencias reflejan cambios genuinos en la estructura de la dependencia espacial y no simplemente un aumento en la prevalencia de las aristas: en CLR, estas aristas que involucran a las células Treg se reducen casi por completo hacia el valor nulo (), mientras que en DII, las mismas aristas tienen una reducción posterior sustancialmente menor (que oscila entre 0,56 y 0,90), lo que indica que la colocalización de las células Treg con el microambiente inmunitario más amplio es una característica espacial definitoria de los tumores DII.

Este patrón es consistente con la biología conocida del subtipo de infiltrado inflamatorio difuso. Los tumores DII se caracterizan por un infiltrado inmunitario heterogéneo y espacialmente disperso, en el que se cree que las células Treg se acumulan en múltiples compartimentos tisulares y suprimen las respuestas inmunitarias efectoras a través de mecanismos dependientes del contacto y mediados por citocinas ([Schürch et al., 2020]). El acoplamiento espacial concurrente de las células Treg con las poblaciones mieloides (macrófagos CD68+ y CD163+) y las poblaciones linfoides (células T CD8+, células T CD4+ de memoria, células B) sugiere una supresión inmunitaria coordinada que opera en múltiples ramas de la respuesta adaptativa e innata. La fuerte interacción Treg↔vasculatura en DII () también puede reflejar nichos perivasculares de células Treg, una disposición espacial que se ha relacionado con una alteración del tráfico de las células inmunitarias hacia el parénquima tumoral.

Las tres aristas significativas de macrófagos CD68+—con células T CD4+ de memoria, estroma y músculo liso—sugieren un eje secundario de reorganización espacial en DII que es independiente de las células Treg. Los macrófagos CD68+ en DII parecen estar incrustados en un contexto estromal y de músculo liso más denso, lo que puede reflejar un microambiente enriquecido en mieloides, coherente con la polarización de los macrófagos asociados al tumor hacia un fenotipo inmunosupresor. La colocalización de los macrófagos CD68+ con las células T CD4+ de memoria en DII podría indicar interacciones de presentación de antígenos que se facilitan espacialmente en este subtipo, pero que se interrumpen o se difunden espacialmente en CLR.

En contraste, el tejido CLR muestra una estructura de interacción espacial globalmente dispersa en todas las 105 aristas probadas, con valores kappa estimados por LMM cercanos a 1 para la mayoría de los pares de tipos de células. Esto es consistente con la arquitectura inflamatoria similar a la enfermedad de Crohn de CLR, en la que la infiltración inmunitaria se organiza en agregados linfoides discretos con una distribución espacial más focal que difusa ([Schürch et al., 2020]). La relativa ausencia de fuertes dependencias espaciales en CLR puede reflejar una respuesta inmunitaria más compartimentada, donde las poblaciones efectoras están segregadas en lugar de estar espacialmente intermezcladas con los elementos reguladores y estromales. Es importante destacar que ninguna arista mostró una interacción significativamente más fuerte en CLR en relación con DII, lo que indica que la remodelación de la red entre los subtipos es direccionalmente asimétrica y refleja predominantemente una ganancia de interacción en DII en lugar de una reestructuración bidireccional de la red de colocalización espacial.

En conjunto, estos hallazgos demuestran que GP-GHS recupera diferencias interpretables y biológicamente coherentes en las redes de interacción espacial a partir de datos de imagen multiplexada. El método no solo identifica qué aristas difieren, sino que cuantifica la dirección y la magnitud de esas diferencias en una escala inferencial bien definida, lo que permite realizar comparaciones entre subgrupos de pacientes definidos patológicamente a una resolución que los análisis estándar de coocurrencia o enriquecimiento de vecindad no proporcionan.

Discusión

Hemos presentado GP-GHS, un marco de regresión bayesiana a nivel de nodo para inferir redes de interacción célula-célula espacialmente variables a partir de datos de imagen tisular multiplexada. El método combina tres componentes que están individualmente motivados y son conjuntamente necesarios: una aproximación de proceso gaussiano en el espacio de Hilbert para la estimación espacial escalable, un prior de herradura de grupo que impone la selección de aristas como una decisión binaria a nivel de grupo en todo el conjunto de coeficientes de base espectral y una estrategia de regresión a nivel de nodo que descompone el problema del modelo gráfico p-dimensional en p regresiones paralelizables. El estudio de simulación estableció que ningún competidor recupera gráficos espacialmente estructurados con una precisión significativa, y el análisis de 140 imágenes de tejido de cáncer colorrectal demostró que el marco recupera redes de interacción centradas en las células Treg que son biológicamente interpretables y que difieren de manera coherente entre dos subgrupos de pacientes definidos patológicamente.

Una pregunta natural es por qué GP-GHS caracteriza solo la presencia de una interacción célula-célula en lugar de su signo, magnitud o patrón espacial. Tres propiedades del entorno de múltiples imágenes dificultan la presentación de estos resúmenes. Primero, el signo y la magnitud de pueden variar continuamente en todo el tejido dentro de una sola imagen: una colocalización en el núcleo del tumor puede coincidir con una exclusión espacial en el borde invasor, por lo que ni un signo global ni una magnitud global están bien definidos. En segundo lugar, no existe un sistema de coordenadas espacial común en las imágenes de diferentes pacientes: el dominio tisular es específico de la imagen, lo que impide el promedio directo de los campos espaciales entre las imágenes. En tercer lugar, las regresiones a nivel de nodo no están restringidas conjuntamente, por lo que las magnitudes de y no están en la misma escala. Dadas estas limitaciones, la puntuación de reducción posterior, que resume la fuerza con la que el prior de grupo suprime todo el campo de interacción hacia cero, es el resumen más portátil e interpretable en todas las imágenes y los pacientes.

El prior de grupo es la contribución de modelado central.

La decisión de diseño más importante en GP-GHS es la colocación de un único parámetro de reducción local por bloque de vecinos en lugar de por coeficiente de base. El estudio de simulación hace que las consecuencias de esta elección sean concretas. Horseshoe estándar, que comparte todos los demás componentes del modelo pero aplica una reducción escalar independiente a cada uno de los coeficientes de base, no logra recuperar ninguna de las aristas verdaderas en todos los niveles de dispersión y ambos tamaños de problema. Este fracaso no es incidental: refleja la falta de correspondencia entre un prior que toma decisiones independientes a nivel de coeficiente y una pregunta científica que requiere una única decisión a nivel de arista. Cuando o 25 funciones de base representan un solo campo de interacción, la reducción independiente acumula oportunidades independientes para una actividad falsa, lo que en la práctica significa que el posterior para los coeficientes individuales nunca se compromete claramente con cero o distinto de cero. El prior de grupo resuelve esto al hacer que los coeficientes dentro de un bloque aumenten y disminuyan juntos bajo un único , de modo que el posterior se concentra fuertemente en el régimen de todo cero (arista ausente) o todo distinto de cero (arista presente). Esta es precisamente la propiedad que produce un valor de TPR superior a 0,85 en el entorno muy disperso donde los competidores logran menos de 0,15.

Limitaciones y direcciones para el trabajo futuro.

Varias limitaciones del marco actual merecen ser abordadas en el trabajo futuro. Primero, el preprocesamiento basado en KDE convierte un proceso de puntos marcado en una entrada de regresión espacial, lo que implica elecciones de selección de ancho de banda y resolución de cuadrícula que podrían influir en las redes inferidas. Un tratamiento más riguroso modelaría directamente el proceso de puntos del tipo de célula, por ejemplo, a través de una formulación de proceso de Cox log-gaussiano donde las superficies de intensidad latentes son las entradas de regresión. Esto propagaría la incertidumbre desde la estimación del proceso de puntos al paso de inferencia de la red en lugar de tratar las superficies KDE como cantidades fijas conocidas, que no lo son.

En segundo lugar, como se demostró con la aplicación de datos de CRC, cuando se tienen múltiples imágenes para comparar, el marco actual se implementa ajustando una red espacial particular para cada imagen y luego resumiendo y comparando las redes en todas las imágenes en un modelo lineal mixto post-hoc. Una formulación jerárquica alternativa modelaría todas las imágenes conjuntamente, colocando un prior sobre la distribución de las redes a nivel de imagen dentro de cada grupo de patología y estimando directamente las redes de consenso a nivel de grupo con cuantificación de la incertidumbre. Un modelo de este tipo aprovecharía la información de las imágenes dentro de un grupo, lo que podría permitir la recuperación de aristas más débiles que están presentes en la mayoría de las imágenes, pero no lo suficientemente fuertes como para sobrevivir al umbral a nivel de imagen. La estructura de correlación intra-paciente, que actualmente manejamos solo en la etapa de prueba a través del efecto aleatorio del paciente en el LMM, podría incorporarse en el modelo conjunto a través de un prior de red a nivel de paciente, pero esto aumenta la complejidad computacional, ya que tenemos GP específicos de la imagen.

En tercer lugar, el costo computacional de GP-GHS en p = 25 es de aproximadamente 31 minutos por imagen, lo que se convierte en un cuello de botella práctico para los conjuntos de datos de imágenes de toda la lámina con cientos de secciones de tejido. El costo dominante son las operaciones de matriz dentro del muestreador de Gibbs. Dos direcciones para la aceleración merecen ser perseguidas: aproximaciones de inferencia variacional al posterior de la herradura de grupo, que reemplazarían el muestreador MCMC con una optimización determinista que se escala de manera más favorable, y aproximaciones de GP dispersas más allá de HSGP que exploten la disposición espacial irregular de las ubicaciones de las células de manera más eficiente que una base de cuadrícula regular.

En cuarto lugar, si bien las funciones de base HSGP imponen la suavidad espacial a través del prior espectral, no incorporan ninguna información estructural a nivel de tejido, como los límites de los compartimentos tisulares o las anotaciones histológicas. En la práctica, las interacciones célula-célula pueden ser cualitativamente diferentes en el núcleo del tumor, el borde invasor y el estroma, y un modelo que pueda asignar redes de interacción separadas a regiones tisulares predefinidas, o que aprenda patrones de interacción espacialmente discontinuos, sería más biológicamente interpretable [Bhadury et al. (2026)]. Extender el marco a GP particionados por dominio o a kernels espacialmente no estacionarios es una dirección natural para el trabajo metodológico futuro.

Finalmente, el marco actual no distingue entre las interacciones de contacto directo célula-célula y las interacciones de señalización parácrina que operan en escalas espaciales más amplias. Las aristas inferidas reflejan las dependencias condicionales espaciales en la escala del ancho de banda KDE y la resolución de la cuadrícula, lo que podría confundir estos dos tipos de interacción mecánicamente distintos. La incorporación explícita de la escala espacial, por ejemplo, incluyendo múltiples escalas de ancho de banda en el KDE o utilizando una expansión de base de múltiples resoluciones, podría proporcionar una descomposición más detallada de los tipos de interacción.

Conclusión.

GP-GHS aborda una necesidad real en el conjunto de herramientas para el análisis de la omica espacial. Los métodos gráficos existentes ignoran por completo las coordenadas espaciales o tratan la estructura espacial como un factor de confusión que debe eliminarse, en lugar de como una fuente de señal que debe modelarse. El prior de grupo "horseshoe" que opera sobre los coeficientes de la base espectral de GP es, según nuestro conocimiento, el primer método que aplica simultáneamente una suavidad espacial a cada campo de interacción inferido y convierte la inclusión de aristas en una decisión coherente a nivel de grupo. El método resultante recupera gráficos con estructura espacial con una precisión sustancialmente mayor que todos los competidores en la simulación, e identifica una red inmunosupresora coherente biológicamente, centrada en las células Treg, en el cáncer colorrectal, que diferencia entre microambientes tumorales patológicamente distintos. A medida que las plataformas de imagen multiplexada se convierten en herramientas estándar en la oncología traslacional y el número de tipos de células perfiladas simultáneamente continúa creciendo, los métodos que puedan representar fielmente la heterogeneidad espacial del microambiente tumoral serán esenciales para traducir los datos de imagen en hipótesis mecanicistas y, en última instancia, en biomarcadores clínicos.

El prior de grupo es la contribución central del modelo.

La decisión de diseño más importante en GP-GHS es la colocación de un único parámetro de reducción local por bloque vecino en lugar de por coeficiente de base. El estudio de simulación hace que las consecuencias de esta elección sean evidentes. El "Horseshoe" estándar, que comparte todos los demás componentes del modelo pero aplica una reducción escalar independiente a cada uno de los coeficientes de base, no logra recuperar ninguna de las aristas verdaderas en todos los niveles de dispersión y en ambos tamaños de problema. Este fallo no es accidental: refleja la falta de correspondencia estructural entre un prior que toma decisiones independientes a nivel de coeficiente y una pregunta científica que requiere una única decisión a nivel de arista. Cuando uno o 25 funciones de base representan un único campo de interacción, la reducción independiente acumula oportunidades independientes de actividad falsa, lo que en la práctica significa que la distribución a posteriori de los coeficientes individuales nunca se compromete claramente con cero o con un valor distinto de cero. El prior de grupo resuelve esto al hacer que los coeficientes dentro de un bloque aumenten y disminuyan juntos bajo un único parámetro, de modo que la distribución a posteriori se concentra fuertemente en el régimen de "todo cero" (arista ausente) o "todo distinto de cero" (arista presente). Esta es precisamente la propiedad que produce un valor de TPR superior a 0,85 en el entorno muy disperso donde los competidores alcanzan un valor inferior a 0,15.

Limitaciones y líneas de trabajo futuro.

Varias limitaciones del marco actual merecen ser abordadas en el trabajo futuro. En primer lugar, el preprocesamiento basado en KDE convierte un proceso de puntos marcado en una entrada de regresión espacial, lo que implica la selección del ancho de banda y las opciones de resolución de la cuadrícula, lo que podría influir en las redes inferidas. Un tratamiento más riguroso modelaría directamente el proceso de puntos del tipo de célula, por ejemplo, mediante una formulación de proceso de Cox log-gaussiano en la que las superficies de intensidad latentes sean las entradas de regresión. Esto propagaría la incertidumbre de la estimación del proceso de puntos a la etapa de inferencia de la red en lugar de tratar las superficies de KDE como cantidades fijas y conocidas, que no lo son.

En segundo lugar, como se demuestra con la aplicación de datos de CRC, cuando se tienen varias imágenes para comparar, el marco actual se implementa ajustando una red espacial particular para cada imagen y, a continuación, resumiendo y comparando las redes entre las imágenes en un modelo lineal mixto post-hoc. Una formulación jerárquica alternativa modelaría todas las imágenes conjuntamente, colocando un prior sobre la distribución de las redes a nivel de imagen dentro de cada grupo de patología y estimando directamente las redes de consenso a nivel de grupo con cuantificación de la incertidumbre. Un modelo de este tipo aprovecharía la información de las imágenes dentro de un grupo, lo que podría permitir recuperar aristas más débiles que están presentes en la mayoría de las imágenes, pero no lo suficientemente fuertes como para sobrevivir al umbral a nivel de imagen. La estructura de correlación intra-paciente, que actualmente solo se trata en la etapa de prueba a través del efecto aleatorio del paciente en el LMM, podría incorporarse al modelo conjunto a través de un prior de red a nivel de paciente, pero esto aumenta la complejidad computacional, ya que tenemos GPs específicos de la imagen.

En tercer lugar, el coste computacional de GP-GHS en p = 25 es de aproximadamente 31 minutos por imagen, lo que se convierte en un cuello de botella práctico para los conjuntos de datos de imagen de toda la lámina con cientos de secciones de tejido. El coste dominante son las operaciones de matriz dentro del muestreador de Gibbs. Merece la pena explorar dos direcciones para la aceleración: aproximaciones de inferencia variacional al prior de grupo "horseshoe", que reemplazarían el muestreador MCMC por una optimización determinista que se escala de forma más favorable, y aproximaciones de GP dispersas más allá de HSGP que exploten la disposición espacial irregular de las ubicaciones de las células de forma más eficiente que una base de cuadrícula regular.

En cuarto lugar, si bien las funciones de base HSGP aplican una suavidad espacial a través del prior espectral, no incorporan ninguna información estructural a nivel de tejido, como los límites de los compartimentos tisulares o las anotaciones histológicas. En la práctica, las interacciones entre células pueden ser cualitativamente diferentes en el núcleo del tumor, el borde invasivo y el estroma, y un modelo que pueda asignar redes de interacción separadas a regiones de tejido predefinidas, o que aprenda patrones de interacción espacialmente discontinuos, sería más interpretable biológicamente [Bhadury et al. (2026)]. Extender el marco a GPs de dominio particionado o a kernels espacialmente no estacionarios es una dirección natural para el trabajo metodológico futuro.

Por último, el marco actual no distingue entre las interacciones de contacto directo entre células y las interacciones de señalización parácrina que operan a escalas espaciales más amplias. Las aristas inferidas reflejan las dependencias condicionales espaciales a la escala del ancho de banda de KDE y la resolución de la cuadrícula, lo que podría confundir estos dos tipos de interacción mecanicísticamente distintos. La incorporación explícita de la escala espacial, por ejemplo, incluyendo múltiples escalas de ancho de banda en el KDE o utilizando una expansión de base de múltiples resoluciones, podría proporcionar una descomposición más detallada de los tipos de interacción.

Conclusión.

GP-GHS aborda una necesidad real en el conjunto de herramientas para el análisis de la omica espacial. Los métodos gráficos existentes ignoran por completo las coordenadas espaciales o tratan la estructura espacial como un factor de confusión que debe eliminarse, en lugar de como una fuente de señal que debe modelarse. El prior de grupo "horseshoe" que opera sobre los coeficientes de la base espectral de GP es, según nuestro conocimiento, el primer método que aplica simultáneamente una suavidad espacial a cada campo de interacción inferido y convierte la inclusión de aristas en una decisión coherente a nivel de grupo. El método resultante recupera gráficos con estructura espacial con una precisión sustancialmente mayor que todos los competidores en la simulación, e identifica una red inmunosupresora coherente biológicamente, centrada en las células Treg, en el cáncer colorrectal, que diferencia entre microambientes tumorales patológicamente distintos. A medida que las plataformas de imagen multiplexada se convierten en herramientas estándar en la oncología traslacional y el número de tipos de células perfiladas simultáneamente continúa creciendo, los métodos que puedan representar fielmente la heterogeneidad espacial del microambiente tumoral serán esenciales para traducir los datos de imagen en hipótesis mecanicistas y, en última instancia, en biomarcadores clínicos.