IntegrateRigor: optimización de la integración sin necesidad de anotaciones para la recuperación de la identidad celular, que revela nichos de interacción entre el cáncer y el sistema inmunitario.

El secuenciamiento de ARN a nivel de célula única a gran escala (scRNA-seq) y la transcriptómica espacial (ST) están transformando la biología al permitir mediciones sistemáticas de la expresión génica a nivel de célula única y con resolución espacial [^[1]–^[3]]. La creciente disponibilidad de grandes atlas de referencia [^[4]–^[6]], junto con la rápida acumulación de conjuntos de datos generados a partir de diversas fuentes experimentales y biológicas, ha ampliado enormemente las oportunidades para estudiar la heterogeneidad celular, la organización de los tejidos y los procesos de la enfermedad a gran escala. Sin embargo, para aprovechar este potencial, cada vez es más necesario integrar datos de múltiples lotes para recuperar identidades celulares comparables, un paso que sigue siendo técnicamente y conceptualmente desafiante.

Aquí, un lote se refiere de manera general a una muestra o grupo de células recolectadas como un replicado o bajo una condición distinta, como de diferentes donantes, puntos de tiempo, laboratorios, plataformas de secuenciación o incluso especies [^[7]]. A lo largo de este manuscrito, identidad celular denota las características transcriptómicas de una célula, incluido su tipo, subtipo o estado [^[8]]. Las diferencias entre lotes pueden surgir de dos fuentes distintas: diferencias en la composición de la identidad celular entre los lotes y diferencias en la expresión génica medida para la misma identidad celular entre los lotes. La primera no impide la identificación de la identidad celular, ya que las células con la misma identidad siguen siendo comparables; la segunda, a la que denominamos variación entre lotes, oscurece la verdadera identidad celular y es el principal desafío para la integración de datos. Es importante destacar que la variación entre lotes no se limita a los artefactos técnicos, sino que también puede reflejar diferencias biológicas entre los lotes, como la expresión diferencial específica del tipo celular. Simplemente agrupar las células de diferentes lotes antes de identificar las identidades celulares no es una opción válida, ya que no tiene en cuenta la variación entre lotes y, por lo tanto, compromete la recuperación de la identidad celular, así como todos los análisis posteriores que dependen de ella, como las pruebas de expresión diferencial, la inferencia de trayectorias y linajes, y el análisis de la comunicación célula-célula [^[9], ^[10]].

Por lo tanto, la integración es esencial para construir un espacio compartido en el que las identidades celulares puedan recuperarse y compararse de manera consistente entre los lotes [^[11], ^[12]]. Sin una integración adecuada, las células con la misma identidad pero de diferentes lotes pueden permanecer separadas debido a la variación entre lotes [^[13], ^[14]]. En esencia, la recuperación de la identidad celular es un requisito previo para cualquier análisis posterior basado en la identidad celular; estos análisis solo deben realizarse después de que se hayan establecido identidades celulares confiables entre los lotes, y no deben estar entrelazados con el paso de integración en sí. Una integración eficaz debe equilibrar, por lo tanto, dos objetivos contrapuestos: preservar la estructura de la identidad celular al tiempo que se minimiza la variación entre lotes. No lograr este equilibrio conduce a dos modos de fallo opuestos. La subintegración deja una variación residual entre lotes, lo que hace que las células se agrupen por lote en lugar de por su verdadera identidad. La sobreintegración, por el contrario, elimina la variación biológicamente significativa entre las identidades celulares y fusiona artificialmente identidades celulares distintas [^[7], ^[15]–^[17]]. Ambos distorsionan la recuperación de la identidad celular y socavan el rigor y la reproducibilidad de cualquier análisis posterior. Se ha desarrollado una amplia gama de métodos computacionales para la integración de datos de células únicas, que incluyen Harmony [^[18]], Seurat-RPCA [^[19]], scVI [^[20]], FastMNN [^[21]] y LIGER [^[22]], entre muchos otros [^[10], ^[23]–^[25]]. Estos métodos están diseñados principalmente para alinear las células con la misma identidad biológica entre los lotes en una incrustación de baja dimensión compartida. Otros enfoques, como Crescendo [^[10]] y CellANOVA [^[26]], toman la incrustación integrada como entrada para realizar la corrección de lotes a nivel de recuento y el análisis de expresión diferencial específico del tipo celular, respectivamente; por lo tanto, sus resultados dependen directamente de la calidad de la incrustación integrada. Sin embargo, a pesar del importante progreso metodológico, la práctica actual de integración sigue siendo limitada en dos aspectos clave, especialmente en el caso común en el que no están disponibles anotaciones precisas de la identidad celular antes de la integración [^[7]].

En primer lugar, la elección de los genes de entrada es un determinante clave, pero a menudo subestimado, de la calidad de la integración para la recuperación de la identidad celular [^[7], ^[27]]. La mayoría de los métodos de integración comienzan con una selección genérica de genes, normalmente utilizando genes altamente variables (HVG) seleccionados ya sea dentro de cada lote y luego agrupados, o de células agrupadas entre los lotes [^[27], ^[28]]. Sin embargo, estos HVG pueden ser sensibles a la variación entre lotes; es decir, sus niveles de expresión para la misma identidad celular pueden diferir sustancialmente entre los lotes. Por ejemplo, muchos HVG reflejan principalmente la actividad o los programas de respuesta al estrés en lugar de la identidad celular intrínseca [^[29], ^[30]], y, por lo tanto, son muy sensibles a la variación entre lotes (por ejemplo, los genes mitocondriales y ribosomales) [^[27], ^[31]]. El uso de tales HVG introduce la variación entre lotes en la entrada de la integración, lo que dificulta la alineación confiable de la identidad celular entre los lotes. En contraste, los genes vinculados a programas conservados de identidad celular muestran patrones de expresión relativamente estables entre donantes, tejidos e incluso especies, lo que proporciona un conjunto de genes candidatos con una base biológica [^[32]–^[35]]. Un estudio de referencia reciente [^[27]] demostró que la selección de genes puede ser tan importante como la elección del método para la calidad de la integración. Aunque se ha propuesto una métrica dependiente de la anotación del tipo celular, los efectos técnicos del grupo (GTE), para cuantificar la variación entre lotes a nivel de gen [^[36]], las anotaciones confiables a menudo no están disponibles antes de la integración, y su uso para la selección de genes introduce circularidad en el flujo de trabajo.

En segundo lugar, las métricas existentes para evaluar el equilibrio entre la mezcla de lotes y la preservación de la identidad celular se dividen en dos categorías, cada una con limitaciones importantes. En primer lugar, las métricas dependientes de la anotación, como el índice local inverso de Simpson del tipo celular (cLISI) [^[13], ^[37]] y el índice de Rand ajustado (ARI) [^[37], ^[38]] para la preservación de la identidad celular, y el LISI de integración consciente del tipo celular (CiLISI) [^[13]] para la mezcla de lotes, son informativas cuando están disponibles anotaciones precisas de la identidad celular. Sin embargo, crean una dependencia circular: la integración se realiza precisamente para recuperar las identidades celulares entre los lotes, pero estas métricas presuponen que ya existen anotaciones confiables. En segundo lugar, las métricas independientes de la anotación, como el LISI de integración (iLISI), la prueba del efecto de lote de los k vecinos más cercanos (kBET), la regresión de componentes principales (PCR) y el ancho de silueta promedio del lote (batch ASW) [^[13], ^[17], ^[18], ^[37]], evitan esta dependencia, pero solo miden la mezcla de lotes sin evaluar la preservación de la identidad celular, y, por lo tanto, no pueden detectar la sobreintegración. Aunque la puntuación de silueta es una métrica independiente de la anotación que puede evaluar la separación de los grupos de células y, por lo tanto, la preservación de la identidad celular, depende del rendimiento del algoritmo de agrupación y se ha demostrado que es poco confiable en el análisis comparativo de la integración de células únicas [^[37]]. Como resultado, ningún enfoque independiente existente evalúa conjuntamente la mezcla de lotes y la preservación de la identidad celular para guiar las decisiones de integración para un conjunto de datos específico. En la práctica, los investigadores a menudo recurren a la inspección visual subjetiva de las parcelas UMAP [^[39]] y a la optimización de hiperparámetros por prueba y error (por ejemplo, la penalización de agrupación de diversidad en Harmony), lo que dificulta la reproducibilidad y la optimización de la integración para un conjunto de datos determinado.

En conjunto, estas limitaciones revelan una brecha crítica: ningún enfoque existente aborda conjuntamente la selección de genes, la optimización de hiperparámetros y la selección de métodos de una manera estadísticamente fundamentada, independiente de la anotación y específica del conjunto de datos, con el objetivo explícito de una recuperación confiable de la identidad celular. Para abordar esta brecha, proponemos IntegrateRigor, un marco basado en datos y agnóstico del método para optimizar la integración de la transcriptómica de células únicas y espaciales para la recuperación de la identidad celular. IntegrateRigor no es un método de integración independiente; más bien, proporciona un envoltorio fundamentado para evaluar y ajustar los métodos de integración existentes. Primero, identifica los genes estables entre lotes (BSG): genes cuyos patrones de expresión entre las identidades celulares permanecen consistentes entre los lotes, y selecciona los HVG dentro de este conjunto estable entre lotes como características de entrada informativas de la identidad celular para la integración. La selección de genes se logra a través de una puntuación de estabilidad del lote (BSS) por gen basada en la probabilidad y sin anotación, que cuantifica la estabilidad del gen entre los lotes sin requerir anotaciones de la identidad celular. Después de la integración, IntegrateRigor evalúa la incrustación de células resultante utilizando dos métricas independientes de la anotación: la puntuación de alineación del lote (BAS), que mide la variación residual entre lotes, y la puntuación de identidad celular (CIS), que mide la preservación de la estructura de la identidad celular. Estas dos puntuaciones se agregan en una puntuación de integración a nivel de conjunto de datos que identifica la configuración óptima (es decir, método y hiperparámetros) para una recuperación confiable de la identidad celular en un conjunto de datos determinado. Esto cambia el enfoque de las pruebas comparativas globales de métodos a la optimización específica del conjunto de datos: aunque los estudios comparativos proporcionan clasificaciones generales valiosas de los métodos de integración [^[13], ^[14], ^[31]], normalmente utilizan configuraciones de parámetros fijas en los conjuntos de datos y, por lo tanto, pueden no reflejar el mejor rendimiento posible de cada método en un conjunto de datos determinado. Al abordar conjuntamente la selección de genes, la optimización de hiperparámetros y la selección de métodos, IntegrateRigor transforma la integración de un flujo de trabajo heurístico de prueba y error en un procedimiento riguroso, independiente de la anotación y adaptable a los datos para la recuperación de la identidad celular.

En una aplicación a la transcriptómica de células únicas y espaciales del cáncer colorrectal [^[10]], IntegrateRigor permitió una recuperación confiable de la identidad celular que reveló nichos de interfaz inmune-cáncer previamente no anotados en el microambiente tumoral, una estructura que estaba oscurecida tanto por la configuración de integración predeterminada como por la incrustación sobreintegrada informada en la literatura anterior. En seis tareas de integración que abarcan diversas fuentes de variación entre lotes (es decir, plataformas de secuenciación, condiciones de muestra, etapa de desarrollo, especies y secciones de tejido espacial), IntegrateRigor mejoró constantemente la recuperación de la identidad celular de cinco métodos ampliamente utilizados (Seurat-RPCA, Harmony, scVI, FastMNN y LIGER). En particular, Harmony, que mostró solo un rendimiento moderado con la configuración predeterminada, logró uno de los mejores resultados generales cuando se combinó con IntegrateRigor, lo que demuestra que la optimización específica del conjunto de datos puede liberar todo el potencial de los métodos existentes para la recuperación de la identidad celular.

Resultados

Descripción general de IntegrateRigor: un marco independiente de la anotación para la selección de genes estables entre lotes y la optimización de la integración para la recuperación de la identidad celular

IntegrateRigor proporciona un envoltorio adaptable a los datos y sin anotación que optimiza la integración de datos para una recuperación confiable de la identidad celular al tiempo que minimiza la subjetividad humana en la optimización de hiperparámetros. Al combinar la selección de BSG con la optimización de hiperparámetros basada en datos, IntegrateRigor mejora la recuperación de la identidad celular al mitigar tanto la subintegración como la sobreintegración, sin requerir anotaciones previas de la identidad celular. Este diseño hace que IntegrateRigor sea ampliamente aplicable a diferentes conjuntos de datos y fácil de usar como un complemento dentro de los flujos de trabajo de integración existentes. Los detalles del método IntegrateRigor se describen en la sección de Métodos.

Como se ilustra en la Figura 1, IntegrateRigor introduce métricas estadísticas independientes de la anotación en tres niveles: (1) una puntuación de estabilidad del lote (BSS) por gen basada en la probabilidad que selecciona genes estables entre lotes e informativos de la identidad celular como características de entrada; (2) métricas de nivel de incrustación basadas en la probabilidad, BAS y CIS, que cuantifican conjuntamente la mezcla de lotes y la preservación de la identidad celular; y (3) una puntuación de integración a nivel de conjunto de datos que identifica el método y la configuración de hiperparámetros óptimos para un conjunto de datos determinado.

En la Etapa 1, IntegrateRigor identifica un conjunto de genes reguladores de lote (BSG) para utilizarlos como características de entrada para la integración posterior (Figura 1a). Aquí, la variación entre lotes se refiere a las diferencias en la expresión génica medida para la misma identidad celular en diferentes lotes, a diferencia de las diferencias en la composición de la identidad celular. Asumimos que un subconjunto de genes contiene información sobre la identidad celular y muestra una variación mínima entre lotes; es decir, su expresión dentro de cada identidad celular no cambia sustancialmente entre lotes (por ejemplo, no son genes expresados diferencialmente específicos de la identidad celular entre lotes), e IntegrateRigor tiene como objetivo identificar y seleccionar estos genes. Intuitivamente, un gen de este tipo debería mostrar los mismos patrones de expresión relativa en todas las identidades celulares en cada lote, incluso si las proporciones de identidades celulares difieren entre los lotes.

Para formalizar esta intuición, construimos una estadística basada en la probabilidad para cuantificar el grado en que la distribución de la expresión de un gen se conserva entre los lotes. La idea es aprovechar un lote de referencia, ya sea definido por el usuario (por ejemplo, un atlas o un conjunto de datos de referencia bien anotado) o, si no se especifica, seleccionado automáticamente como el lote con la mayor diversidad de identidades celulares, aproximada por el mayor número de grupos de células dentro de un lote, para guiar la recuperación de la estructura de identidad celular compartida entre los lotes. Nos referimos a los lotes restantes como lotes de consulta. Específicamente, para cada gen y cada lote de consulta, ajustamos dos modelos de mezcla: un modelo restringido por el lote de referencia que fija los componentes de la mezcla (correspondientes a las identidades celulares) según lo estimado en el lote de referencia, al tiempo que permite que solo las proporciones de los componentes difieran, y un modelo no restringido en el que tanto los componentes como las proporciones son libres de variar con respecto al lote de referencia. La diferencia en el ajuste del modelo entre estos dos modelos, medida por la probabilidad logarítmica normalizada, cuantifica el grado de variación entre lotes para ese gen en ese lote de consulta; el promedio de esta diferencia en todos los lotes de consulta produce la puntuación BSS por gen, que captura la consistencia con la que se mantiene el patrón de expresión de un gen en todos los lotes (Figura 1b). Los genes con puntuaciones más altas se seleccionan como BSG utilizando un umbral adaptativo de datos determinado mediante la detección del punto de inflexión, lo que produce características de entrada que conservan preferentemente las señales de identidad celular y, al mismo tiempo, reducen la variación entre lotes. Los HVG se seleccionan entonces dentro de los BSG como características de entrada finales, asegurando que sean tanto informativos para la identidad celular como estables entre los lotes.

En la Etapa 2, IntegrateRigor identifica el resultado de integración óptimo entre las incrustaciones candidatas generadas bajo diferentes configuraciones de hiperparámetros (Figura 1a). Cada incrustación se evalúa utilizando dos métricas complementarias sin anotación que corresponden a los dos objetivos de la integración para la recuperación de la identidad celular. La puntuación BAS cuantifica el grado de alineación entre lotes en cada dimensión de la incrustación, donde los valores más altos reflejan una mejor mezcla de lotes. La puntuación CIS mide la fuerza con la que cada dimensión de la incrustación conserva la estructura de identidad celular biológicamente significativa, cuantificada por el grado en que la distribución de los valores de la incrustación en esa dimensión refleja distintas subpoblaciones celulares en lugar de un fondo unimodal. Específicamente, para cada dimensión de la incrustación, calculamos la diferencia de probabilidad logarítmica entre un modelo restringido y un modelo no restringido por separado para el lote de referencia y cada lote de consulta, y luego hacemos el promedio en todos los lotes, con cada puntuación por lote normalizada por el número de células para garantizar la comparabilidad. Estas puntuaciones por dimensión se agregan en todas las dimensiones de la incrustación para obtener las puntuaciones BAS y CIS (Figura 1b–c). La puntuación de integración a nivel de conjunto de datos se define como la suma de BAS y CIS, lo que captura directamente el equilibrio entre la mezcla de lotes y la preservación de la identidad celular: BAS penaliza la subintegración, mientras que CIS penaliza la sobreintegración. Dado que ambas puntuaciones se definen como diferencias de probabilidad logarítmica normalizadas, están en una escala comparable y se pueden sumar directamente, sin necesidad de un parámetro de ponderación adicional. La configuración de hiperparámetros que maximiza la puntuación de integración se selecciona entonces como la configuración óptima para la recuperación de la identidad celular en el conjunto de datos dado.

La validación empírica de los modelos estadísticos subyacentes a IntegrateRigor se encuentra en la Sección de Métodos.### Las puntuaciones de estabilidad del lote identifican los genes estables entre lotes y mejoran la recuperación de la identidad celular al excluir los genes afectados por la variación de la muestra y la variación técnica

La puntuación BSS, calculada en la Etapa 1 de IntegrateRigor, cuantifica la estabilidad con la que se mantiene el patrón de expresión de un gen en todos los lotes dentro de la misma identidad celular. Los genes con una alta puntuación BSS muestran patrones de expresión consistentes en todos los lotes dentro de la misma identidad celular; los genes con una baja puntuación BSS muestran cambios importantes entre lotes, lo que indica que su expresión está fuertemente influenciada por la variación entre lotes y puede oscurecer la estructura de identidad celular compartida, lo que dificulta una recuperación confiable de la identidad celular. Validamos la puntuación BSS y evaluamos el beneficio posterior de la selección de genes basada en la puntuación BSS utilizando el conjunto de datos de enfermedad hepática asociada con la insuficiencia intestinal (IFALD), que comprende tres donantes [^[40]].

Los genes clasificados en alto por la puntuación BSS mostraron patrones de expresión más concordantes en todos los lotes, mientras que los genes clasificados en bajo mostraron una variación sustancial entre lotes (Figura 2a). El análisis de componentes principales de las células agrupadas confirmó esta separación: los genes con grandes cargas en el componente principal 1 (PC1), que está dominado por el lote, tendieron a tener una puntuación BSS más baja y no se clasificaron como estables entre lotes, mientras que los genes con grandes cargas en el PC2, que está dominado por el tipo de célula, tendieron a tener una puntuación BSS más alta y se clasificaron como estables entre lotes (Figura S4). Para comparar la puntuación BSS con una métrica dependiente de la anotación existente, comparamos la puntuación BSS con la puntuación GTE [^[36]], que requiere anotaciones de tipo de célula. La puntuación BSS se correlacionó linealmente con la puntuación GTE en la escala logarítmica en todos los conjuntos de datos (Pearson r = 0,50–0,69; todos los p < 2,2 × 10−16; Figura S5), lo que confirma que la puntuación BSS captura una noción similar de variación entre lotes que esta métrica dependiente de la anotación, sin requerir etiquetas de tipo de célula.

IntegrateRigor clasifica los genes por la puntuación BSS y utiliza un umbral de punto de inflexión adaptativo de datos [^[41]] para excluir los genes que claramente no son estables entre lotes, conservando el resto como BSG para la selección posterior de HVG. Este enfoque es deliberadamente conservador: en lugar de seleccionar agresivamente solo los genes más estables, elimina los que son claramente problemáticos y, al mismo tiempo, conserva un conjunto amplio para la selección de HVG. En el conjunto de datos IFALD, solo 331 de 14 349 genes se excluyeron como inestables entre lotes, dejando 14 018 BSG. Tanto en la configuración predeterminada como en la configuración restringida por BSG, se utilizó el mismo número de 2000 HVG como características de entrada; la única diferencia fue si esos HVG se extraían del conjunto de genes completo o se restringían a los BSG. Cuando se aplicó Harmony a los 2000 HVG predeterminados, la incrustación integrada mostró una clara subintegración, con las células del mismo tipo que permanecían separadas por lote. Restringir la selección de HVG a los BSG mejoró sustancialmente la recuperación de la identidad celular, aumentando la puntuación CiLISI de 1,33 a 1,47 y la puntuación ARI de 0,77 a 0,89 (Figura 2b). Se observaron mejoras consistentes en todos los cinco métodos de integración evaluados (FastMNN, Harmony, LIGER, scVI y Seurat-RPCA), lo que demuestra que la selección de BSG es una mejora modular y agnóstica del método que se puede integrar en las canalizaciones existentes.

Para comprender por qué los genes inestables entre lotes dificultan la recuperación de la identidad celular, examinamos su contenido biológico y técnico. Los tres donantes de IFALD difieren en edad, método de preservación de tejidos y condición de la enfermedad (Figura 2c). Tanto en el PCA no integrado como en el UMAP de Harmony predeterminado, el donante C54 (un adulto sano, congelado) se separó claramente de los otros dos donantes, lo que sugiere que la variación dominante entre lotes refleja las características a nivel de muestra en lugar del estado de la enfermedad en sí (Figura 2b; Figura S4). En consonancia con esto, los genes inestables entre lotes se enriquecieron con términos de ontología génica relacionados con la coagulación y la fibrinolisis, las respuestas inflamatorias y agudas, y los procesos lipídicos y metabólicos (Figura 2e); estos son procesos que probablemente están asociados con la edad y la preservación en lugar de las identidades celulares compartidas [^[42]–^[44]]. Los genes inestables entre lotes también se enriquecieron con genes mitocondriales y ribosomales, que se sabe que reflejan la calidad de la muestra, el estrés y la actividad transcripcional global en lugar de la identidad celular; por el contrario, estos genes estaban poco representados entre los BSG (Tabla S2). Cabe destacar que no todos los genes ribosomales eran uniformemente inestables: un subconjunto con funciones reguladoras especializadas se clasificó constantemente como estable entre lotes en todos los conjuntos de datos, lo que es consistente con los análisis gen por gen anteriores [^[36]] (Material suplementario; Figura S6 suplementaria).

En conjunto, estos resultados demuestran que la selección de BSG mejora la recuperación de la identidad celular al filtrar los genes dominados por los efectos específicos de la muestra o el ruido técnico, al tiempo que conserva los genes que codifican la estructura de identidad celular compartida en todos los lotes.### Las puntuaciones de integración guían la optimización de hiperparámetros para equilibrar la subintegración y la sobreintegración para la recuperación de la identidad celular

Los métodos de integración de última generación (por ejemplo, Seurat-RPCA, Harmony, scVI y FastMNN) producen una incrustación celular de baja dimensión como resultado de su integración, e IntegrateRigor evalúa cada dimensión de esta incrustación utilizando dos métricas complementarias sin anotación: BAS, que cuantifica el grado de alineación entre lotes, y CIS, que mide la fuerza con la que la dimensión conserva la estructura de identidad celular biológicamente distinta. Para calcular BAS y CIS, IntegrateRigor utiliza un modelo de mezcla gaussiana para evaluar la distribución de los valores en cada dimensión de la incrustación. Sin embargo, las incrustaciones LIGER no se capturaron bien con este modelo (véase la Sección de Métodos), por lo que LIGER se excluyó del ajuste de hiperparámetros y los análisis de selección de métodos.

En el conjunto de datos de PBMC estimuladas con IFN-β [^[45]] (13 999 células, dos lotes), BAS y CIS fueron directamente interpretables a nivel de dimensión en la incrustación de Harmony: las dimensiones con una alta puntuación BAS mostraron una fuerte mezcla de lotes, mientras que las dimensiones con una alta puntuación CIS resolvieron claramente los principales tipos de células inmunitarias (Figura 2f). Esto confirma que las dos puntuaciones capturan aspectos complementarios de la calidad de la integración y pueden respaldar el análisis exploratorio de las dimensiones individuales de la incrustación integrada.

Basándose en esto, IntegrateRigor calcula el promedio de las puntuaciones BAS y CIS por dimensión en todas las dimensiones de la incrustación y define la puntuación de integración como su suma. En el conjunto de datos de PBMC estimuladas con IFN-β, Harmony predeterminada mostró una clara subintegración, con las células que permanecieron parcialmente separadas por lote, lo que comprometió la recuperación de la identidad celular (Figura 2g). En contraste, el valor de hiperparámetro seleccionado por la puntuación de integración, incluso sin la selección de BSG, mejoró la mezcla de lotes (la puntuación CiLISI aumentó de 1,54 a 1,65) y, al mismo tiempo, mejoró sustancialmente la recuperación de la identidad celular (la puntuación ARI aumentó de 0,73 a 0,94) (Figura 2g). En esencia, esta optimización no requirió anotaciones de identidad celular y se basó únicamente en la estructura intrínseca de la incrustación integrada.

Una ventaja clave de la evaluación conjunta de la mezcla de lotes y la preservación de la identidad celular es la capacidad de proteger contra la sobreintegración, un modo de fallo que las métricas de mezcla de lotes sin anotación no pueden detectar por sí solas. Las métricas existentes sin anotación, como iLISI, PCR, batch ASW y kBET, recompensan una mayor mezcla de lotes, pero no evalúan si se preserva la estructura de la identidad celular. Para ilustrar esta limitación, analizamos un conjunto de datos de células inmunitarias humanas con dos lotes (uno derivado de tejido de médula ósea y el otro de PBMC) [^[13]]. Luego, evaluamos las incrustaciones de Harmony en un rango de valores de hiperparámetros. Las cuatro métricas de mezcla de lotes mostraron tendencias similares en los hiperparámetros y, aunque kBET seleccionó un valor de hiperparámetro ligeramente diferente al de iLISI, PCR y batch ASW, todos estos ajustes favorecieron una mayor mezcla de lotes y se desviaron del hiperparámetro en el que la puntuación de integración de IntegrateRigor alcanzó su punto máximo. Específicamente, a lo largo del camino del hiperparámetro, BAS aumentó con el valor del hiperparámetro, mientras que CIS disminuyó (Figura S10), lo que demuestra el equilibrio entre la mezcla de lotes y la preservación de la identidad celular, una consideración ausente en las métricas de mezcla de lotes existentes (Figura S10). En el hiperparámetro seleccionado mediante la maximización de la puntuación de integración, las células de diferentes lotes estaban bien alineadas, mientras que las principales poblaciones inmunitarias permanecieron claramente separadas (CiLISI = 1,43, cLISI = 1,16, ARI = 0,65). En contraste, el hiperparámetro seleccionado por iLISI, PCR, ASW y kBET logró una mezcla de lotes similar (CiLISI = 1,40), pero exhibió una sobreintegración: varias identidades celulares distintas se fusionaron (cLISI aumentó a 1,26 y ARI disminuyó a 0,61) (Figura 2h). Aquí, un valor más alto de cLISI indica una peor preservación de la identidad celular.

En conjunto, estos resultados demuestran que la puntuación de integración identifica los ajustes de hiperparámetros que logran un equilibrio equilibrado entre la mezcla de lotes y la preservación de la identidad celular para una recuperación fiable de la identidad celular, mitigando tanto la subintegración como la sobreintegración sin requerir anotaciones previas de la identidad celular.### IntegrateRigor mejora constantemente la recuperación de la identidad celular en diversos conjuntos de datos y métodos

Para evaluar la generalizabilidad de IntegrateRigor, lo evaluamos en una colección diversa de tareas de integración que abarcan múltiples fuentes de variación entre lotes: condiciones de muestra (PBMC estimuladas con IFN-β frente a controles [^[45]]; donantes de IFALD con distinta edad, preservación y contexto de la enfermedad [^[40]]), plataformas de secuenciación (células inmunitarias humanas [^[13]]), etapas de desarrollo (corazón de pollo [^[46]]), cortes de tejido (transcriptómica espacial del cerebro de ratón [^[47]]) y especies (páncreas humano-ratón [^[48]]). Esta diversidad nos permite probar si IntegrateRigor puede mejorar la recuperación de identidades celulares comparables entre lotes en entornos que van más allá de la corrección de lotes de una sola célula de rutina, incluidos los casos en los que las diferencias técnicas y biológicas hacen que la integración sea especialmente desafiante.

En todas estas tareas, IntegrateRigor mejoró constantemente la recuperación de la identidad celular en relación con la canalización predeterminada, como lo demuestra una mejor mezcla de lotes (CiLISI) junto con una preservación de la identidad celular mantenida o mejorada (ARI) (Figura 3a). En el conjunto de datos de PBMC estimuladas con IFN-β, la aplicación de IntegrateRigor con scVI mejoró la mezcla de lotes (CiLISI aumentó de 1,69 a 1,78) al tiempo que se mantuvo una alta recuperación de la identidad celular (ARI de 0,82 a 0,83), lo que refleja una mejor mezcla de células estimuladas y de control sin pérdida de la estructura del tipo celular. En IFALD, donde Harmony predeterminado dejó las células del mismo tipo parcialmente separadas por donante, IntegrateRigor mejoró sustancialmente tanto la mezcla de lotes como la recuperación de la identidad celular (CiLISI: 1,36 → 1,47; ARI: 0,77 → 0,88). En el conjunto de datos del corazón de pollo, FastMNN con IntegrateRigor mejoró de CiLISI = 2,49 a 2,63 y ARI de 0,62 a 0,63. En el conjunto de datos del páncreas humano-ratón, la mezcla de lotes de Seurat-RPCA aumentó de CiLISI = 1,88 a 2,07, mientras que ARI se mantuvo alto en 0,95. Los valores métricos completos y las visualizaciones UMAP para todos los conjuntos de datos se proporcionan en las Figuras Suplementarias S15-S21. En conjunto, estos resultados demuestran que IntegrateRigor mejora la recuperación de la identidad celular en una amplia gama de escenarios de conjuntos de datos, desde casos en los que la incrustación predeterminada ya es razonable hasta ejemplos con una subintegración pronunciada.

Las ganancias se resumen en la Figura 3b, que muestra la mejora relativa de IntegrateRigor sobre la integración predeterminada en los diferentes métodos y conjuntos de datos. La mezcla de lotes (CiLISI) mejoró en casi todos los entornos, con las mayores ganancias en los conjuntos de datos con una fuerte variación entre lotes, como PBMC estimuladas con IFN-β, IFALD y páncreas humano-ratón. La fuerte mejora en el conjunto de datos del páncreas humano-ratón muestra que IntegrateRigor también puede optimizar la integración entre especies, donde las diferencias entre los lotes no se limitan al ruido técnico, sino que incluyen diferencias transcriptómicas a nivel de especie. Este resultado respalda la aplicabilidad más amplia de IntegrateRigor: al seleccionar genes estables entre lotes y ajustar la fuerza de la integración, el marco puede recuperar identidades celulares comparables incluso cuando los conjuntos de datos están separados por grandes diferencias biológicas y técnicas. En esencia, estas mejoras en la alineación de lotes se acompañaron de una preservación de ARI mantenida o mejorada, lo que confirma que el aumento de la mezcla no se produjo a expensas de la pérdida de la estructura de la identidad celular (Figura 3c). En particular, Harmony mostró la mayor mejora y se convirtió en el método con mejor rendimiento en general cuando fue optimizado por IntegrateRigor, a pesar de su rendimiento predeterminado moderado, lo que demuestra cómo la optimización específica del conjunto de datos puede permitir una mejor recuperación de la identidad celular.

También examinamos si la puntuación de integración, además de guiar la selección de hiperparámetros dentro de cada método, puede ayudar a distinguir entre los métodos de integración. Esto es inherentemente difícil porque la integración no tiene un estándar de oro absoluto y diferentes métricas capturan diferentes aspectos del rendimiento [^[37], ^[49]]. Por lo tanto, no afirmamos que IntegrateRigor seleccione un método óptimo único, sino que preguntamos si su clasificación coincide con las métricas establecidas dependientes de la anotación del alineamiento de lotes y la preservación de la identidad celular.

En cinco de los seis conjuntos de datos, Seurat-RPCA y Harmony fueron los dos métodos con las puntuaciones de integración más altas (Figura S22). La excepción fue el conjunto de datos del corazón de pollo, para el cual Seurat-RPCA y FastMNN obtuvieron las puntuaciones más altas, con FastMNN logrando la puntuación de integración más alta (Figura S22). Este patrón es generalmente consistente con la clasificación general en la Figura 3c, donde Seurat-RPCA y Harmony se agrupan en la región superior derecha (CiLISI alto y ARI alto), y con la fuerte separación de tipos celulares observada para FastMNN en el conjunto de datos del corazón de pollo (Figura S18). Estas observaciones respaldan la puntuación de integración como un proxy útil sin anotación para la recuperación de la identidad celular a nivel de método y el alineamiento de lotes.### IntegrateRigor descubre nichos de interfaz cáncer-inmune en el cáncer colorrectal mediante la integración de datos de transcriptómica de una sola célula y espacial

Una integración bien fundamentada para una recuperación fiable de la identidad celular puede permitir descubrimientos biológicos que permanecen ocultos bajo opciones de hiperparámetros subóptimas. Ilustramos esto en un conjunto de datos de cáncer colorrectal perfilado con transcriptómica de una sola célula y espacial [^[10]] (192 166 células; dos cortes de tejido MERSCOPE espacial y un lote de una sola célula 10x Genomics). El estudio original [^[10]] utilizó Harmony como un paso previo para obtener una incrustación integrada para Crescendo, un método de corrección de recuento posterior, y proporcionó anotaciones de tipo celular basadas en la incrustación de Harmony (Figura 4a). Sin embargo, ni el código ni la incrustación integrada utilizada en ese estudio se publicaron. Por lo tanto, reproducimos la incrustación de Harmony con la configuración predeterminada y la encontramos marcadamente diferente de la incrustación publicada: la incrustación original mostró una extensa mezcla de lotes (Figura 4a), mientras que nuestra incrustación de Harmony predeterminada mostró una clara separación residual de lotes, particularmente entre las células epiteliales cancerosas (Figura S23). Esta discrepancia destaca un desafío común: sin un criterio bien definido, no está claro qué incrustación preserva mejor las identidades celulares biológicamente distintas.

La aplicación de IntegrateRigor a este conjunto de datos produjo una incrustación de Harmony optimizada que difería tanto de la original como del resultado predeterminado, logrando un mejor equilibrio entre la mezcla de lotes y la preservación de la identidad celular (Figura 4a). Cuando examinamos si las anotaciones de tipo celular del estudio original seguían siendo consistentes con la incrustación optimizada, encontramos una discordancia sustancial: un subconjunto de células originalmente anotadas como células epiteliales cancerosas se agrupó entre las poblaciones inmunitarias en la incrustación de IntegrateRigor. Esto nos motivó a realizar una anotación independiente de tipo celular en la incrustación optimizada utilizando la agrupación de Louvain y la expresión de genes marcadores.

En nuestra reanotación, las células que se desplazaron de las células epiteliales cancerosas a las células inmunitarias se mapearon de nuevo a sus ubicaciones físicas en los datos de transcriptómica espacial (Figura 4b). Estas células se localizaron específicamente en el límite entre las células epiteliales cancerosas canónicas (anotadas como epiteliales cancerosas tanto en las anotaciones originales como en las revisadas) y las células inmunitarias canónicas (anotadas como inmunitarias en ambas anotaciones). Nos referimos a estas poblaciones límite localizadas espacialmente como nichos de interfaz cáncer-inmune.

Confirmamos su identidad como poblaciones de interfaz mediante el análisis de expresión diferencial utilizando la expresión génica observada de un solo lote antes de la integración, para evitar artefactos de integración (Figura 4c). En comparación con las células epiteliales cancerosas canónicas, cada nicho se enriqueció con marcadores de su tipo celular inmunitario correspondiente; en comparación con las células inmunitarias canónicas, cada nicho retuvo una fuerte expresión de marcadores epiteliales (KRT8, MUC13, PERP). Esta firma dual, similar a la inmune en relación con el epitelio, y similar al epitelio en relación con la inmune, define la identidad del nicho de interfaz. Específicamente, el nicho cáncer-B expresó CD19, MS4A1 y CXCR5; el nicho cáncer-mast expresó CPA3, MS4A2 y KIT; el nicho cáncer-T expresó CCL5, CD3D y KLRB1; el nicho cáncer-plasma expresó MANF, CCR10 y MZB1; y el nicho cáncer-mieloide expresó CSF3R, FCGR3B y IL1A. Se esperan perfiles de transcripción mixtos en los datos de transcriptómica espacial, donde el contacto directo entre células, la segmentación celular imperfecta y la contaminación local de transcritos pueden contribuir.

El análisis de enriquecimiento funcional de los genes expresados diferencialmente entre cada nicho y su correspondiente población inmune canónica (excluyendo los marcadores epiteliales compartidos) reveló programas específicos del tipo de célula inmune, coherentes con funciones activas en la interfaz cáncer-inmune (Figura 4d; Figura S24) [^[50]–^[52]]. El nicho cáncer-mieloide se enriqueció en la “vía de señalización mediada por citocinas”, “regulación positiva de la producción de citocinas”, “adhesión célula-célula de leucocitos” y “regulación de la activación de las células T”, lo que sugiere un programa de comunicación inmune y coordinación local en el microambiente tumoral [^[53], ^[54]]. El nicho cáncer-B estuvo dominado por términos relacionados con el procesamiento de antígenos y la presentación de antígenos, lo que indica un estado especializado para el manejo de antígenos en la interfaz del cáncer epitelial [^[55]–^[57]]. El nicho cáncer-T se enriqueció en citotoxicidad mediada por leucocitos, producción de interleucina-10 y adhesión célula-célula, lo que sugiere un estado inmunorregulador de las células T modulado por el contacto epitelial [^[54], ^[58]]. El nicho cáncer-plasma se enriqueció en programas de presentación de antígenos y respuesta inmune activada, más allá de la función canónica de las células plasmáticas [^[59]]. El nicho cáncer-mastocitos se enriqueció en términos relacionados con la formación de la memoria inmunológica y la presentación de antígenos, lo que sugiere su participación en una regulación inmune adaptativa más amplia [^[60]–^[62]]. En todos los nichos, los programas enriquecidos reflejan las funciones características del tipo de célula inmune correspondiente, lo que es coherente con su posición en la interfaz cáncer-inmune.

Estos nichos de interfaz no se detectaron ni con el embedding original ni con el embedding predeterminado de Harmony. El embedding original los fusionó en un grupo amplio, enmascarando su identidad intermedia y conduciendo a su anotación como células epiteliales cancerosas únicamente. El embedding predeterminado de Harmony no logró alinearlos de manera consistente entre los lotes, lo que impidió su reconocimiento como tipos específicos de nichos. En contraste, el embedding optimizado de IntegrateRigor reveló estos nichos de interfaz cáncer-inmune y respaldó su interpretación a través de la localización espacial, las firmas duales de marcadores cáncer-inmune y el enriquecimiento funcional específico del tipo de célula inmune. Este ejemplo demuestra que la elección de los hiperparámetros de integración no es simplemente un detalle técnico: al optimizar la integración para la recuperación de la identidad celular, IntegrateRigor puede revelar poblaciones celulares de interés biológico que, de otro modo, permanecerían ocultas.

Descripción general de IntegrateRigor: un marco de trabajo sin anotación para la selección de genes estable entre lotes y la optimización de la integración para la recuperación de la identidad celular

IntegrateRigor proporciona un wrapper adaptable a los datos y sin anotación que optimiza la integración de datos para una recuperación fiable de la identidad celular, al tiempo que minimiza la subjetividad humana en el ajuste de los hiperparámetros. Al combinar la selección de genes estables entre lotes (BSG) con la optimización de hiperparámetros basada en los datos, IntegrateRigor mejora la recuperación de la identidad celular al mitigar tanto la subintegración como la sobreintegración, sin requerir anotaciones previas de la identidad celular. Este diseño hace que IntegrateRigor sea ampliamente aplicable en diferentes conjuntos de datos y fácil de usar como un complemento dentro de las canalizaciones de integración existentes. Los detalles del método IntegrateRigor se describen en la sección de Métodos.

Como se ilustra en la Figura 1, IntegrateRigor introduce métricas estadísticas sin anotación en tres niveles: (1) una estadística de similitud entre lotes (BSS) por gen, basada en la probabilidad, que selecciona genes estables entre lotes e informativos para la identidad celular como características de entrada; (2) métricas a nivel de embedding basadas en la probabilidad, BAS y CIS, que cuantifican conjuntamente la mezcla entre lotes y la preservación de la identidad celular; y (3) una puntuación de integración a nivel de conjunto de datos que identifica el método y la configuración de hiperparámetros óptimos para un conjunto de datos determinado.

En la Etapa 1, IntegrateRigor identifica un conjunto de BSG para usar como características de entrada para la integración posterior (Figura 1a). Aquí, la variación entre lotes se refiere a las diferencias en la expresión génica medida para la misma identidad celular en diferentes lotes, a diferencia de las diferencias en la composición de la identidad celular. Asumimos que un subconjunto de genes contiene información sobre la identidad celular y muestra una variación mínima entre lotes; es decir, su expresión dentro de cada identidad celular no cambia sustancialmente entre los lotes (por ejemplo, no son genes expresados diferencialmente entre lotes para una identidad celular específica), e IntegrateRigor tiene como objetivo identificar y seleccionar estos genes. Intuitivamente, un gen de este tipo debería mostrar los mismos patrones de expresión relativa entre las identidades celulares en cada lote, incluso si las proporciones de las identidades celulares difieren entre los lotes.

Para formalizar esta intuición, construimos una estadística basada en la probabilidad para cuantificar el grado en que se conserva la distribución de la expresión de un gen entre los lotes. La idea es aprovechar un lote de referencia, ya sea definido por el usuario (por ejemplo, un atlas o un conjunto de datos de referencia bien anotado) o, si no se especifica, seleccionado automáticamente como el lote con la mayor diversidad de identidades celulares, aproximada por el mayor número de grupos de células dentro de un lote, para guiar la recuperación de la estructura de identidad celular compartida entre los lotes. Nos referimos a los lotes restantes como lotes de consulta. Específicamente, para cada gen y cada lote de consulta, ajustamos dos modelos de mezcla: un modelo restringido por la referencia que fija los componentes de la mezcla (correspondientes a las identidades celulares) según lo estimado en el lote de referencia, al tiempo que solo permite que difieran las proporciones de los componentes, y un modelo no restringido en el que tanto los componentes como las proporciones son libres de variar con respecto al lote de referencia. La diferencia en el ajuste del modelo entre estos dos modelos, medida por la probabilidad logarítmica normalizada, cuantifica el grado de variación entre lotes para ese gen en ese lote de consulta; el promedio de esta diferencia en todos los lotes de consulta produce la BSS por gen, que captura la consistencia con la que se mantiene el patrón de expresión de un gen entre los lotes (Figura 1b). Los genes con puntuaciones más altas se seleccionan como BSG utilizando un umbral adaptable a los datos determinado por la detección del punto de inflexión, lo que produce características de entrada que conservan preferentemente las señales de identidad celular y reducen la variación entre lotes. Los HVG se seleccionan entonces dentro de los BSG como características de entrada finales, lo que garantiza que sean tanto informativos para la identidad celular como estables entre los lotes.

En la Etapa 2, IntegrateRigor identifica el resultado de integración óptimo entre los embeddings candidatos generados bajo diferentes configuraciones de hiperparámetros (Figura 1a). Cada embedding se evalúa utilizando dos métricas complementarias sin anotación que corresponden a los dos objetivos de la integración para la recuperación de la identidad celular. La BAS cuantifica el grado de alineación entre lotes en cada dimensión del embedding, donde los valores más altos reflejan una mejor mezcla entre lotes. La CIS mide la fuerza con la que cada dimensión del embedding preserva la estructura biológicamente significativa de la identidad celular, cuantificada por el grado en que la distribución de los valores del embedding en esa dimensión refleja distintas subpoblaciones celulares en lugar de un fondo unimodal. Específicamente, para cada dimensión del embedding, calculamos la diferencia de probabilidad logarítmica entre un modelo restringido y un modelo no restringido por separado para el lote de referencia y cada lote de consulta, y luego hacemos el promedio entre los lotes, con cada puntuación por lote normalizada por el número de células para garantizar la comparabilidad. Estas puntuaciones por dimensión se agregan en todas las dimensiones del embedding para obtener la BAS y la CIS (Figura 1b–c). La puntuación de integración a nivel de conjunto de datos se define como la suma de la BAS y la CIS, lo que captura directamente el equilibrio entre la mezcla entre lotes y la preservación de la identidad celular: la BAS penaliza la subintegración, mientras que la CIS penaliza la sobreintegración. Dado que ambas puntuaciones se definen como diferencias de probabilidad logarítmica normalizadas, están en una escala comparable y se pueden sumar directamente, sin requerir un parámetro de ponderación adicional. La configuración de hiperparámetros que maximiza la puntuación de integración se selecciona entonces como la configuración óptima para la recuperación de la identidad celular en el conjunto de datos dado.

La validación empírica de los modelos estadísticos subyacentes a IntegrateRigor se encuentra en la sección de Métodos.

Las puntuaciones de estabilidad entre lotes identifican genes estables entre lotes y mejoran la recuperación de la identidad celular al excluir los genes afectados por la variación de la muestra y la variación técnica

La BSS, calculada en la Etapa 1 de IntegrateRigor, cuantifica la estabilidad con la que se mantiene el patrón de expresión de un gen entre los lotes dentro de la misma identidad celular. Los genes con una BSS alta muestran patrones de expresión consistentes entre los lotes dentro de la misma identidad celular; los genes con una BSS baja muestran cambios sustanciales entre lotes, lo que indica que su expresión está fuertemente influenciada por la variación entre lotes y puede enmascarar la estructura de identidad celular compartida, lo que dificulta una recuperación fiable de la identidad celular. Validamos la BSS y evaluamos el beneficio posterior de la selección de genes basada en la BSS utilizando el conjunto de datos de la enfermedad hepática asociada con la insuficiencia intestinal (IFALD), que comprende tres donantes [^[40]].

Los genes clasificados en alto por la BSS mostraron patrones de expresión más concordantes entre los lotes, mientras que los genes clasificados en bajo mostraron una variación sustancial entre lotes (Figura 2a). El análisis de componentes principales de las células agrupadas confirmó esta separación: los genes con grandes cargas en el componente principal 1 (PC1), que está dominado por el lote, tendieron a tener una BSS más baja y no se clasificaron como estables entre lotes, mientras que los genes con grandes cargas en el PC2, que está dominado por el tipo de célula, tendieron a tener una BSS más alta y se clasificaron como estables entre lotes (Figura S4). Para comparar la BSS con una métrica dependiente de la anotación existente, comparamos la BSS con la puntuación GTE [^[36]], que requiere anotaciones del tipo de célula. La BSS se correlacionó linealmente con la GTE en la escala logarítmica en todos los conjuntos de datos (Pearson r = 0,50–0,69; todos los p < 2,2 × 10−16; Figura S5), lo que confirma que la BSS captura una noción similar de variación entre lotes que esta métrica dependiente de la anotación, sin requerir etiquetas del tipo de célula.

IntegrateRigor clasifica los genes por BSS y utiliza un umbral de punto de inflexión adaptable a los datos [^[41]] para excluir los genes que claramente no son estables entre lotes, conservando el resto como BSG para la selección posterior de HVG. Este enfoque es deliberadamente conservador: en lugar de seleccionar agresivamente solo los genes más estables, elimina los que son claramente problemáticos y conserva un grupo amplio para la selección de HVG. En el conjunto de datos IFALD, solo 331 de 14 349 genes se excluyeron como inestables entre lotes, dejando 14 018 BSG. Tanto en la configuración predeterminada como en la restringida a BSG, se utilizaron el mismo número de 2000 HVG como características de entrada; la única diferencia fue si esos HVG se extraían del conjunto de genes completo o se restringían a los BSG. Cuando se aplicó Harmony a los 2000 HVG predeterminados, el embedding integrado mostró una clara subintegración, con las células del mismo tipo que permanecían separadas por lote. Restringir la selección de HVG a los BSG mejoró sustancialmente la recuperación de la identidad celular, aumentando el CiLISI de 1,33 a 1,47 y el ARI de 0,77 a 0,89 (Figura 2b). Se observaron mejoras consistentes en los cinco métodos de integración evaluados (FastMNN, Harmony, LIGER, scVI y Seurat-RPCA), lo que demuestra que la selección de BSG es una mejora modular y agnóstica al método que se puede conectar a las canalizaciones existentes.

Para comprender por qué los genes inestables entre lotes dificultan la recuperación de la identidad celular, examinamos su contenido biológico y técnico. Los tres donantes de IFALD difieren en edad, método de preservación de los tejidos y condición de la enfermedad (Figura 2c). Tanto en el PCA no integrado como en el UMAP de Harmony predeterminado, el donante C54 (un adulto sano, congelado) se separó claramente de los otros dos donantes, lo que sugiere que la variación dominante entre lotes refleja las características a nivel de muestra en lugar del estado de la enfermedad en sí (Figura 2b; Figura S4). En consonancia con esto, los genes inestables entre lotes se enriquecieron con términos de ontología génica relacionados con la coagulación y la fibrinolisis, las respuestas agudas e inflamatorias, y los procesos lipídicos y metabólicos (Figura 2e), procesos que probablemente estén asociados con la edad y la preservación en lugar de identidades celulares compartidas [^[42]–^[44]]. Los genes inestables entre lotes también se enriquecieron con genes mitocondriales y ribosomales, que se sabe que reflejan la calidad de la muestra, el estrés y la actividad transcripcional global en lugar de la identidad celular; por el contrario, estos genes estaban subrepresentados entre los BSG (Tabla S2). Cabe destacar que no todos los genes ribosomales eran uniformemente inestables: un subconjunto con funciones reguladoras especializadas se clasificó de forma consistente como estable entre lotes en los diferentes conjuntos de datos, lo que coincide con los análisis previos a nivel de gen [^[36]] (Material Suplementario; Figura S6 del Material Suplementario).

En conjunto, estos resultados demuestran que la selección de BSG mejora la recuperación de la identidad celular al filtrar los genes dominados por efectos específicos de la muestra o ruido técnico, al tiempo que conserva los genes que codifican la estructura de identidad celular compartida entre los lotes.

Las puntuaciones de integración guían la optimización de hiperparámetros para equilibrar la subintegración y la sobreintegración para la recuperación de la identidad celular

Los métodos de integración de última generación (por ejemplo, Seurat-RPCA, Harmony, scVI y FastMNN) producen una incrustación celular de baja dimensión como resultado de su integración, e IntegrateRigor evalúa cada dimensión de esta incrustación utilizando dos métricas complementarias independientes de la anotación: BAS, que cuantifica el grado de alineación entre lotes, y CIS, que mide la fuerza con la que la dimensión preserva la estructura de identidad celular biológicamente distinta. Para calcular BAS y CIS, IntegrateRigor utiliza un modelo de mezcla gaussiana para evaluar la distribución de los valores en cada dimensión de la incrustación. Sin embargo, las incrustaciones LIGER no se capturaron bien con este modelo (véase la sección Métodos), por lo que LIGER se excluyó de los análisis de ajuste de hiperparámetros y selección de métodos.

En el conjunto de datos de PBMC estimuladas con IFN-β [^[45]] (13.999 células, dos lotes), BAS y CIS fueron directamente interpretables a nivel de dimensión en la incrustación de Harmony: las dimensiones con un alto BAS mostraron una fuerte mezcla entre lotes, mientras que aquellas con un alto CIS resolvieron claramente los principales tipos de células inmunitarias (Figura 2f). Esto confirma que las dos puntuaciones capturan aspectos complementarios de la calidad de la integración y pueden respaldar el análisis exploratorio de las dimensiones individuales de la incrustación integrada.

Basándose en esto, IntegrateRigor calcula la media de BAS y CIS a nivel de dimensión en todas las dimensiones de la incrustación y define la puntuación de integración como su suma. En el conjunto de datos de PBMC estimuladas con IFN-β, Harmony predeterminada mostró una clara subintegración, con las células permaneciendo parcialmente separadas por lote, lo que comprometió la recuperación de la identidad celular (Figura 2g). En contraste, el valor de hiperparámetro seleccionado por la puntuación de integración, incluso sin la selección de BSG, mejoró la mezcla entre lotes (CiLISI aumentó de 1,54 a 1,65) al tiempo que mejoró sustancialmente la recuperación de la identidad celular (ARI aumentó de 0,73 a 0,94) (Figura 2g). Lo crucial es que esta optimización no requirió ninguna anotación de identidad celular y se basó enteramente en la estructura intrínseca de la incrustación integrada.

Una ventaja clave de la puntuación conjunta de la mezcla entre lotes y la preservación de la identidad celular es la capacidad de protegerse contra la sobreintegración, un modo de fallo que las métricas de mezcla entre lotes independientes de la anotación por sí solas no pueden detectar. Las métricas independientes de la anotación existentes, como iLISI, PCR, batch ASW y kBET, recompensan una mayor mezcla entre lotes, pero no evalúan si se preserva la estructura de identidad celular. Para ilustrar esta limitación, analizamos un conjunto de datos de células inmunitarias humanas con dos lotes (uno derivado de tejido de médula ósea y el otro de PBMC) [^[13]]. A continuación, evaluamos las incrustaciones de Harmony en un rango de valores de hiperparámetros. Las cuatro métricas de mezcla entre lotes mostraron tendencias similares en los hiperparámetros y, aunque kBET seleccionó un valor de hiperparámetro ligeramente diferente al de iLISI, PCR y batch ASW, todos estos ajustes favorecieron una mayor mezcla entre lotes y divergieron del hiperparámetro en el que la puntuación de integración de IntegrateRigor alcanzó su punto máximo. Específicamente, a lo largo del camino del hiperparámetro, BAS aumentó con el valor del hiperparámetro, mientras que CIS disminuyó (Figura S10), lo que demuestra el equilibrio entre la mezcla entre lotes y la preservación de la identidad celular, una consideración ausente en las métricas de mezcla entre lotes existentes (Figura S10). En el hiperparámetro seleccionado por la maximización de la puntuación de integración, las células de diferentes lotes estaban bien alineadas, mientras que las principales poblaciones inmunitarias permanecían claramente separadas (CiLISI = 1,43, cLISI = 1,16, ARI = 0,65). En contraste, el hiperparámetro seleccionado por iLISI, PCR, ASW y kBET logró una mezcla entre lotes similar (CiLISI = 1,40), pero exhibió una sobreintegración: varias identidades celulares distintas se fusionaron (cLISI aumentó a 1,26 y ARI disminuyó a 0,61) (Figura 2h). Aquí, un valor más alto de cLISI indica una peor preservación de la identidad celular.

En conjunto, estos resultados demuestran que la puntuación de integración identifica los ajustes de hiperparámetros que logran un equilibrio equilibrado entre la mezcla entre lotes y la preservación de la identidad celular para una recuperación fiable de la identidad celular, mitigando tanto la subintegración como la sobreintegración sin requerir anotaciones previas de la identidad celular.

IntegrateRigor mejora constantemente la recuperación de la identidad celular en diversos conjuntos de datos y métodos

Para evaluar la generalizabilidad de IntegrateRigor, lo evaluamos en una colección diversa de tareas de integración que abarcan múltiples fuentes de variación entre lotes: condiciones de la muestra (PBMC estimuladas con IFN-β frente a PBMC de control [^[45]]; donantes de IFALD con diferentes edades, preservación y contexto de la enfermedad [^[40]]), plataformas de secuenciación (células inmunitarias humanas [^[13]]), etapas de desarrollo (corazón de pollo [^[46]]), cortes de tejido (transcriptómica espacial del cerebro de ratón [^[47]]) y especies (páncreas humano-ratón [^[48]]). Esta diversidad nos permite probar si IntegrateRigor puede mejorar la recuperación de identidades celulares comparables entre lotes en entornos que van más allá de la corrección de lotes de una sola célula de rutina, incluidos los casos en los que las diferencias técnicas y biológicas hacen que la integración sea especialmente difícil.

En todas estas tareas, IntegrateRigor mejoró constantemente la recuperación de la identidad celular en relación con la canalización predeterminada, como lo demuestra una mejor mezcla entre lotes (CiLISI) junto con una preservación de la identidad celular mantenida o mejorada (ARI) (Figura 3a). En el conjunto de datos de PBMC estimuladas con IFN-β, la aplicación de IntegrateRigor con scVI mejoró la mezcla entre lotes (CiLISI aumentó de 1,69 a 1,78) al tiempo que se mantuvo una alta recuperación de la identidad celular (ARI de 0,82 a 0,83), lo que refleja una mejor mezcla de las células estimuladas y de control sin pérdida de la estructura de los tipos de células. En IFALD, donde Harmony predeterminada dejó las células del mismo tipo parcialmente separadas por donante, IntegrateRigor mejoró sustancialmente tanto la mezcla entre lotes como la recuperación de la identidad celular (CiLISI: 1,36 → 1,47; ARI: 0,77 → 0,88). En el conjunto de datos del corazón de pollo, FastMNN con IntegrateRigor mejoró de CiLISI = 2,49 a 2,63 y ARI de 0,62 a 0,63. En el conjunto de datos del páncreas humano-ratón, la mezcla entre lotes de Seurat-RPCA aumentó de CiLISI = 1,88 a 2,07, mientras que ARI se mantuvo alto en 0,95. Los valores métricos completos y las visualizaciones UMAP de todos los conjuntos de datos se proporcionan en las Figuras S15-S21 del Material Suplementario. En conjunto, estos resultados demuestran que IntegrateRigor mejora la recuperación de la identidad celular en una amplia gama de escenarios de conjuntos de datos, desde casos en los que la incrustación predeterminada ya es razonable hasta ejemplos con una marcada subintegración.

Las ganancias se resumen en la Figura 3b, que muestra la mejora relativa de IntegrateRigor sobre la integración predeterminada en los diferentes métodos y conjuntos de datos. La mezcla entre lotes (CiLISI) mejoró en casi todos los entornos, con las mayores ganancias en los conjuntos de datos con una fuerte variación entre lotes, como PBMC estimuladas con IFN-β, IFALD y páncreas humano-ratón. La fuerte mejora en el conjunto de datos del páncreas humano-ratón muestra que IntegrateRigor también puede optimizar la integración entre especies, donde las diferencias entre los lotes no se limitan al ruido técnico, sino que incluyen diferencias transcripómicas a nivel de especie. Este resultado apoya la aplicabilidad más amplia de IntegrateRigor: al seleccionar los genes estables entre lotes y ajustar la fuerza de la integración, el marco puede recuperar identidades celulares comparables incluso cuando los conjuntos de datos están separados por grandes diferencias biológicas y técnicas. Lo crucial es que estas mejoras en la alineación entre lotes se acompañaron de una preservación de ARI mantenida o mejorada, lo que confirma que el aumento de la mezcla no se produjo a expensas de la pérdida de la estructura de la identidad celular (Figura 3c). Cabe destacar que Harmony mostró la mayor mejora y se convirtió en el método con mejor rendimiento en general cuando fue optimizado por IntegrateRigor, a pesar de su rendimiento moderado predeterminado, lo que demuestra cómo la optimización específica del conjunto de datos puede permitir una mejor recuperación de la identidad celular.

También examinamos si la puntuación de integración, además de guiar la selección de hiperparámetros dentro de cada método, puede ayudar a distinguir entre los métodos de integración. Esto es inherentemente difícil porque la integración no tiene un estándar de oro absoluto y diferentes métricas capturan diferentes aspectos del rendimiento [^[37], ^[49]]. Por lo tanto, no afirmamos que IntegrateRigor seleccione un método óptimo único, sino que preguntamos si su clasificación coincide con las métricas de anotación dependientes establecidas de la alineación entre lotes y la preservación de la identidad celular.

En cinco de los seis conjuntos de datos, Seurat-RPCA y Harmony fueron los dos métodos con las puntuaciones de integración más altas (Figura S22). La excepción fue el conjunto de datos del corazón de pollo, para el cual Seurat-RPCA y FastMNN obtuvieron las puntuaciones más altas, con FastMNN logrando la puntuación de integración más alta (Figura S22). Este patrón es ampliamente coherente con la clasificación general en la Figura 3c, donde Seurat-RPCA y Harmony se agrupan en la región superior derecha (alto CiLISI y alto ARI), y con la fuerte separación de los tipos de células observada para FastMNN en el conjunto de datos del corazón de pollo (Figura S18). Estas observaciones respaldan la puntuación de integración como un proxy útil independiente de la anotación para la recuperación y la alineación de la identidad celular a nivel de método.

IntegrateRigor descubre nichos de interfaz tumor-inmune en el cáncer colorrectal mediante la integración de datos de transcriptómica de una sola célula y espacial

Una integración de principios para una recuperación fiable de la identidad celular puede permitir descubrimientos biológicos que permanecen ocultos bajo opciones de hiperparámetros subóptimas. Ilustramos esto en un conjunto de datos de cáncer colorrectal perfilado con transcriptómica de una sola célula y espacial [^[10]] (192.166 células; dos cortes de tejido MERSCOPE espacial y un lote de 10x Genomics de una sola célula). El estudio original [^[10]] utilizó Harmony como paso previo para obtener una incrustación integrada para Crescendo, un método de corrección de recuento posterior, y proporcionó anotaciones de tipo de célula basadas en la incrustación de Harmony (Figura 4a). Sin embargo, ni el código ni la incrustación integrada utilizada en ese estudio se publicaron. Por lo tanto, reproducimos la incrustación de Harmony con la configuración predeterminada y encontramos que era marcadamente diferente de la incrustación publicada: la incrustación original mostró una extensa mezcla entre lotes (Figura 4a), mientras que nuestra incrustación de Harmony predeterminada mostró una clara separación residual entre lotes, particularmente entre las células epiteliales cancerosas (Figura S23). Esta discrepancia destaca un desafío común: sin un criterio de principios, no está claro qué incrustación preserva mejor las identidades celulares biológicamente distintas.

La aplicación de IntegrateRigor a este conjunto de datos produjo un embedding Harmony optimizado que difería tanto del original como del resultado predeterminado, logrando un mejor equilibrio entre la mezcla de lotes y la preservación de la identidad celular (Figura 4a). Cuando examinamos si las anotaciones de tipo celular del estudio original se mantenían consistentes con el embedding optimizado, encontramos una discordancia sustancial: un subconjunto de células originalmente anotadas como células epiteliales cancerosas se agrupó entre las poblaciones inmunitarias en el embedding de IntegrateRigor. Esto nos motivó a realizar una anotación independiente de los tipos celulares en el embedding optimizado utilizando el clustering de Louvain y la expresión de genes marcadores.

En nuestra reanotación, las células que pasaron de ser células epiteliales cancerosas a células inmunitarias se volvieron a mapear a sus ubicaciones físicas en los datos de transcriptómica espacial (Figura 4b). Estas células se localizaron específicamente en la frontera entre las células epiteliales cancerosas canónicas (anotadas como epiteliales tanto en las anotaciones originales como en las revisadas) y las células inmunitarias canónicas (anotadas como inmunitarias en ambas anotaciones). Nos referimos a estas poblaciones de frontera localizadas espacialmente como nichos de interfaz cáncer-inmunidad.

Confirmamos su identidad como poblaciones de interfaz mediante un análisis de expresión diferencial utilizando la expresión génica observada de un solo lote antes de la integración, para evitar artefactos de integración (Figura 4c). En comparación con las células epiteliales cancerosas canónicas, cada nicho se enriqueció con marcadores de su tipo de célula inmunitaria correspondiente; en comparación con las células inmunitarias canónicas, cada nicho conservó una fuerte expresión de marcadores epiteliales (KRT8, MUC13, PERP). Esta firma dual (similar a la inmunitaria en relación con el epitelio, similar al epitelio en relación con la inmunidad) define la identidad del nicho de interfaz. Específicamente, el nicho cáncer-B expresó CD19, MS4A1 y CXCR5; el nicho cáncer-mast expresó CPA3, MS4A2 y KIT; el nicho cáncer-T expresó CCL5, CD3D y KLRB1; el nicho cáncer-plasma expresó MANF, CCR10 y MZB1; y el nicho cáncer-mieloide expresó CSF3R, FCGR3B y IL1A. Se espera que estos perfiles transcripcionales mixtos se encuentren en los datos de transcriptómica espacial, donde el contacto directo entre células, la segmentación celular imperfecta y la contaminación local de transcritos pueden contribuir.

El análisis de enriquecimiento funcional de los genes expresados diferencialmente entre cada nicho y su población inmunitaria canónica correspondiente (excluyendo los marcadores epiteliales compartidos) reveló programas específicos del tipo de célula inmunitaria, consistentes con funciones activas en la interfaz cáncer-inmunidad (Figura 4d; Figura S24) [^[50]–^[52]]. El nicho cáncer-mieloide se enriqueció con la "vía de señalización mediada por citocinas", la "regulación positiva de la producción de citocinas", la "adhesión célula-célula de leucocitos" y la "regulación de la activación de las células T", lo que es consistente con un programa de comunicación inmunitaria y coordinación local en el microambiente tumoral [^[53], ^[54]]. El nicho cáncer-B estuvo dominado por términos de procesamiento y presentación de antígenos, lo que es consistente con un estado especializado para el manejo de antígenos en la interfaz del cáncer epitelial [^[55]–^[57]]. El nicho cáncer-T se enriqueció con citotoxicidad mediada por leucocitos, producción de interleucina-10 y adhesión célula-célula, lo que sugiere un estado de célula T inmunorreguladora moldeado por el contacto epitelial [^[54], ^[58]]. El nicho cáncer-plasma se enriqueció con programas de presentación de antígenos y respuesta inmunitaria activada más allá de la función canónica de las células plasmáticas [^[59]]. El nicho cáncer-mast se enriqueció con términos relacionados con la formación de memoria inmunológica y la presentación de antígenos, lo que sugiere su participación en una regulación inmunitaria adaptativa más amplia [^[60]–^[62]]. En todos los nichos, los programas enriquecidos reflejan las funciones características del tipo de célula inmunitaria correspondiente, lo que es consistente con sus posiciones en la interfaz cáncer-inmunidad.

Estos nichos de interfaz no eran detectables ni en el embedding original ni en el embedding Harmony predeterminado. El embedding original los fusionó en un grupo amplio, enmascarando su identidad intermedia y provocando que se anotaran solo como células epiteliales cancerosas. El embedding Harmony predeterminado no logró alinearlos de manera consistente entre los lotes, lo que impidió que se reconocieran como tipos específicos de nichos. En contraste, el embedding Harmony optimizado de IntegrateRigor reveló estos nichos de interfaz cáncer-inmunidad y respaldó su interpretación a través de la localización espacial, las firmas duales de marcadores cáncer-inmunidad y el enriquecimiento funcional específico del tipo de célula inmunitaria. Este ejemplo demuestra que la elección de los hiperparámetros de integración no es simplemente un detalle técnico: al optimizar la integración para la recuperación de la identidad celular, IntegrateRigor puede revelar poblaciones celulares de interés biológico que, de otro modo, permanecerían ocultas.

Discusión

IntegrateRigor reformula la integración de datos, pasando de ser un paso de preprocesamiento heurístico a un problema de optimización guiado estadísticamente y basado en el conjunto de datos para la recuperación de la identidad celular. En lugar de introducir un nuevo método de integración, proporciona un marco agnóstico del método para mejorar los métodos existentes a través de la selección de genes estable entre lotes y la optimización de hiperparámetros sin necesidad de anotaciones. En diversas tareas de transcriptómica unicelular y espacial, esta estrategia mejoró la recuperación de la identidad celular al equilibrar mejor la mezcla de lotes y la preservación de la identidad celular, reduciendo tanto la subintegración como la sobreintegración, y revelando una estructura biológicamente significativa que no era aparente con la configuración predeterminada.

Nuestros resultados destacan dos aspectos poco apreciados de las prácticas de integración. Primero, la selección de genes es en sí misma un determinante importante de la recuperación de la identidad celular, un punto que recibe mucha menos atención que la elección del método de integración, a pesar de que los genes seleccionados definen la información sobre la que opera el método al alinear las células entre los lotes. Al identificar los genes estables entre lotes (BSG) y excluir los genes dominados por la variación entre lotes, IntegrateRigor mejora la capacidad de diversos métodos para recuperar identidades celulares comparables entre los lotes. En segundo lugar, el rendimiento de la integración depende no solo del método elegido, sino también de cómo se configure ese método para un conjunto de datos determinado. Harmony mostró solo un rendimiento moderado con la configuración predeterminada, pero se convirtió en uno de los métodos con mejor rendimiento después de la optimización basada en IntegrateRigor. Esto ilustra que las clasificaciones de los métodos pueden cambiar sustancialmente cuando la selección de genes y la optimización de hiperparámetros se realizan de manera más sistemática. En términos más generales, estos hallazgos abogan por pasar de las canalizaciones de integración estáticas y universales a la optimización específica del conjunto de datos, donde la pregunta central es cómo la configuración del método apoya mejor la recuperación de la identidad celular para un conjunto de datos en particular.

La calidad del embedding integrado tiene consecuencias que se extienden más allá del propio paso de integración. Dado que la recuperación de la identidad celular es un requisito previo para cualquier análisis posterior basado en la identidad celular, la fiabilidad del embedding integrado determina la calidad de todos los análisis que siguen. Por ejemplo, Crescendo [^[10]] y CellANOVA [^[26]] toman el embedding integrado directamente como entrada para realizar la corrección de lotes a nivel de recuento y el análisis de expresión diferencial específico del tipo celular en función de diseños experimentales particulares, respectivamente. Dado que sus representaciones y descomposiciones del estado celular se estiman a partir del embedding integrado, la fiabilidad de sus resultados depende directamente de la calidad de la recuperación de la identidad celular en ese embedding. Un embedding subintegrado conserva la estructura residual del lote que puede confundir la estimación del estado celular; un embedding sobreintegrado difumina los estados biológicamente distintos, lo que provoca una mezcla de tipos celulares que corrompe los análisis específicos del tipo celular. Proporcionar a estos métodos posteriores un embedding bien calibrado, que respalde una recuperación fiable de la identidad celular entre los lotes, es, por lo tanto, un requisito previo para una inferencia posterior precisa.

Vale la pena aclarar la relación entre IntegrateRigor y CellANOVA [^[26]]. CellANOVA recupera las señales biológicas asociadas con la variación entre lotes dentro de la identidad celular que es ortogonal al diseño experimental, una variación que puede eliminarse inadvertidamente durante la integración. En ese sentido, aborda una forma de sobreintegración, pero su objetivo es distinto del nuestro. IntegrateRigor aborda la tarea previa de la recuperación de la identidad celular: alinear de forma fiable las células con la misma identidad entre los lotes antes de realizar cualquier análisis posterior. Para esta tarea de recuperación de la identidad celular, la variación entre lotes dentro de la identidad celular se trata como una subintegración residual que debe corregirse, no como una señal que debe recuperarse. Además, CellANOVA aún depende de un embedding integrado existente para estimar su codificación del estado celular, por lo que su rendimiento sigue siendo sensible a la calidad de ese embedding previo. IntegrateRigor aborda esta dependencia proporcionando un embedding optimizado como entrada.

La importancia del embedding integrado para la recuperación de la identidad celular se ilustró directamente en la aplicación del cáncer colorrectal. La aplicación de IntegrateRigor para obtener un embedding Harmony optimizado permitió una recuperación más fiable de la identidad celular, lo que a su vez reveló nichos de interfaz cáncer-inmunidad previamente no identificados en los que distintas poblaciones de células inmunitarias mostraron firmas transcripcionales y espaciales consistentes con la interacción directa con las células epiteliales malignas. Sin una selección cuidadosa de genes y una optimización de parámetros, el embedding Harmony predeterminado dejó una estructura residual del lote que impidió una recuperación coherente de la identidad celular, mientras que el embedding del estudio original sobrecorrigió y oscureció estas poblaciones de interfaz. Este ejemplo ilustra cómo la optimización de la integración para la recuperación de la identidad celular da forma directamente a qué poblaciones celulares se pueden distinguir.

Se podrían realizar varias ampliaciones para fortalecer aún más este marco. En primer lugar, la formulación actual se centra en la preservación de la identidad celular como el objetivo principal de la optimización de la integración, tratando la variación entre lotes en la identidad celular como el objetivo de la corrección. Un siguiente paso lógico es incorporar covariables conocidas a nivel de muestra (por ejemplo, condición de tratamiento, estado de la enfermedad o etapa de desarrollo) para distinguir los efectos genéricos de lote de la variación biológica estructurada de interés. Con esta información del diseño experimental, los mismos principios de optimización podrían extenderse más allá de la tarea de anotación de la identidad celular para optimizar directamente la incrustación integrada para análisis posteriores, como los realizados por Crescendo y CellANOVA: por ejemplo, optimizando la integración que maximiza la potencia para la expresión diferencial específica de la condición, al tiempo que se controla el lote. Esta extensión refleja el orden natural del análisis: primero, obtener identidades celulares fiables y, luego, aprovecharlas para la inferencia basada en la identidad celular. En segundo lugar, IntegrateRigor se basa en dos supuestos de modelado: un modelo de mezcla binomial negativa (NBMM) para los recuentos de expresión por gen y un modelo de mezcla gaussiana (GMM) para los valores por dimensión en la incrustación integrada. En la práctica, estos supuestos deben verificarse antes de aplicar IntegrateRigor a un nuevo conjunto de datos, para garantizar que ambos modelos proporcionen aproximaciones adecuadas y que las puntuaciones resultantes sean interpretables. Esto es particularmente importante para el GMM subyacente a la puntuación de integración, ya que las incrustaciones producidas por algunos métodos, como LIGER, pueden desviarse de este supuesto. Por lo tanto, el trabajo futuro podría explorar modelos latentes más flexibles o estimadores de densidad no paramétricos que capturen mejor la geometría de diferentes incrustaciones. En tercer lugar, nuestros hallazgos también sugieren que los estudios comparativos futuros deberían evaluar los métodos no solo con la configuración predeterminada, sino también después de la optimización específica del conjunto de datos, ya que el rendimiento predeterminado puede subestimar sustancialmente la calidad que se puede lograr para la recuperación de la identidad celular con una configuración bien definida. En cuarto lugar, los mismos principios podrían extenderse a la integración entre modalidades, donde la preservación de la identidad celular compartida entre diferentes espacios de características presenta desafíos adicionales, y a entornos semi-supervisados en los que están disponibles anotaciones de referencia parciales, de manera similar a los enfoques guiados por referencia, como scANVI [^[24], ^[63]]. En quinto lugar, hasta ahora hemos supuesto que un único lote de referencia contiene todas las identidades celulares relevantes, un supuesto que es cada vez más razonable dada la disponibilidad de grandes atlas que pueden servir como referencia. En entornos donde cada lote contiene solo un subconjunto de identidades celulares, una extensión natural sería construir una referencia sintética mediante la fusión y la alineación de los componentes del modelo de mezcla entre lotes, en lugar de depender de un único lote de referencia; esto requerirá un desarrollo metodológico en el trabajo futuro.

Métodos

Primero, consideramos un entorno de dos lotes que consta de un lote de referencia y un lote de consulta, y luego extendemos el marco a múltiples lotes. El lote de referencia sirve como ancla para la estructura de identidad celular compartida y se supone que contiene un conjunto más diverso de identidades celulares, mientras que el lote de consulta puede contener solo un subconjunto de identidades celulares.

Sea la matriz de recuento de expresión génica del lote de referencia denotada por y la del lote de consulta por , donde y denotan el número de células en los lotes de referencia y consulta, respectivamente, y denota el número de genes.

Selección de genes estables entre lotes (BSG) para la integración de datos

Modelado estadístico para la expresión génica en los lotes de referencia y consulta

Nos centramos en un solo gen y suprimimos el índice del gen para simplificar la notación. Sea denota los recuentos de expresión en el lote , con los tamaños de biblioteca correspondientes . Aquí, se refiere al lote de referencia y se refiere al lote de consulta.

Para identificar los genes que preservan la identidad celular entre los lotes, modelamos la expresión génica utilizando un modelo de mezcla binomial negativa (NB), con componentes de mezcla que representan distintas identidades celulares. Este modelo se basa en amplias pruebas previas de que la distribución NB modela eficazmente los datos de recuento de expresión génica sobre-dispersos dentro de una identidad celular homogénea [^[64]]. La validación de este supuesto del modelo se proporciona en la subsección Métodos "Validación de los modelos estadísticos subyacentes a IntegrateRigor". Específicamente, suponemos que para cualquier gen dado

y

donde . Aquí, denota la proporción de mezcla del componente en el lote , con , y y denotan los parámetros de media y dispersión correspondientes para ese componente NB, respectivamente. Suponemos que estos parámetros específicos del componente son distintos entre los componentes de la mezcla; es decir,

Intuitivamente, si un gen refleja la identidad celular en lugar de la variación entre lotes, su patrón de expresión entre las identidades celulares debe conservarse entre los lotes, incluso cuando las proporciones de esas identidades celulares difieren. En este caso, se supone que la distribución de expresión génica en el lote de consulta comparte un subconjunto de los componentes de la mezcla en el lote de referencia, al tiempo que solo se permite que las proporciones de mezcla difieran.

Específicamente, en ausencia de variación entre lotes, los componentes de la mezcla en el lote de consulta (2) son un subconjunto de los componentes de la mezcla en el lote de referencia (1). Entonces, la distribución de expresión génica del lote de consulta se puede volver a parametrizar utilizando , al tiempo que solo se permite que las proporciones de los componentes varíen:

a lo que nos referimos como el modelo restringido por referencia. Aquí, denota la proporción de mezcla del componente de referencia en el modelo restringido por referencia para el lote de consulta. Lo distinguimos de , que se refiere a la proporción de mezcla del componente de consulta en el modelo (2), y los componentes de referencia y los componentes de consulta pueden no compartir los mismos índices.

Por el contrario, la variación entre lotes induce cambios en los parámetros de los componentes, de modo que . En consecuencia, la distribución de la consulta ya no puede explicarse adecuadamente solo con la estructura de la mezcla de referencia y, en cambio, requiere el modelo completo (2). Por lo tanto, definimos la puntuación de estabilidad entre lotes (BSS) para el gen en función de qué tan bien el lote de consulta puede explicarse mediante el modelo restringido por referencia en relación con el modelo completo.

Definición de la puntuación de estabilidad entre lotes (BSS)

Para distinguir estos dos escenarios, comparamos dos modelos, el modelo restringido por referencia (3) y el modelo completo (2), para el lote de consulta, ambos ajustados utilizando el algoritmo de Expectativa-Maximización (EM) [^[65]]. Bajo el modelo restringido por referencia (3), los parámetros de los componentes se estiman primero a partir del lote de referencia mediante EM y luego se tratan como fijos al evaluar el lote de consulta, con solo las proporciones de mezcla reestimadas por EM. Bajo el modelo completo (2), todos los parámetros se estiman directamente a partir del lote de consulta utilizando EM.

Sea denota la función de verosimilitud del modelo restringido por referencia y de la función de verosimilitud del modelo completo, ambos evaluados en función de las estimaciones de los parámetros del lote de consulta. Entonces, definimos la puntuación de estabilidad entre lotes (BSS) como

que toma un valor negativo ya que es estrictamente menor que ya que el espacio de parámetros del modelo restringido por referencia es un subespacio estricto del espacio de parámetros del modelo completo.

Una BSS alta (menos negativa) indica que el modelo restringido por referencia explica bien los datos de la consulta, lo que sugiere que el patrón de expresión del gen se transfiere entre los lotes y, por lo tanto, es estable entre lotes. Por el contrario, una BSS baja (más negativa) indica que el modelo completo proporciona un ajuste sustancialmente mejor, lo que sugiere una variación entre lotes. Los genes estables entre lotes (BSG) se seleccionan en función de la distribución empírica de los valores de BSS utilizando un umbral adaptativo de datos determinado por un criterio de codo [^[41]].

Puntuación de integración para un conjunto de datos y un método de integración dados

Modelado estadístico para la incrustación integrada

Sea y denotan las incrustaciones de baja dimensión de los lotes de referencia y consulta de una incrustación integrada producida por un método de integración dado, donde denota la j-ésima dimensión de la incrustación para el lote , donde denota el lote de referencia y denota el lote de consulta.

Para cada dimensión de incrustación , suponemos que

y

donde . Aquí, denota la proporción de mezcla del componente en la dimensión de incrustación para el lote , con , y y denotan los parámetros de media y varianza correspondientes para ese componente gaussiano, respectivamente. Dada nuestra observación de que el modelo de mezcla gaussiana refleja distintas identidades celulares en algunas dimensiones de la incrustación integrada (ver la subsección Métodos "Validación de los modelos estadísticos subyacentes a IntegrateRigor"), suponemos que los componentes de la mezcla en este modelo corresponden a identidades celulares y que los parámetros específicos del componente son distintos; es decir,

Definición de la puntuación de alineación entre lotes (BAS)

La BAS cuantifica la mezcla de lotes en el espacio de incrustación integrada. Para cada dimensión de incrustación , siguiendo la misma lógica que la definición de BSS, comparamos un modelo restringido por referencia,

con el modelo completo (4) para las incrustaciones del lote de consulta en la misma dimensión. Bajo el modelo restringido por referencia, los parámetros de los componentes se estiman a partir del lote de referencia y luego se tratan como fijos al evaluar el lote de consulta, con solo las proporciones de mezcla reestimadas. Bajo el modelo completo, todos los parámetros se estiman directamente a partir del lote de consulta.

Sea denota la función de verosimilitud del modelo restringido por referencia y denota la función de verosimilitud del modelo completo, ambas evaluadas en las incrustaciones de la consulta en la dimensión . Definimos la BAS por dimensión como

La BAS general se obtiene promediando sobre todas las dimensiones de incrustación:

Una BAS más alta (menos negativa) indica que el modelo restringido por referencia explica bien las incrustaciones de la consulta, lo que sugiere una mejor alineación entre los lotes y una eliminación más eficaz de la variación entre lotes. Por el contrario, una BAS más baja (más negativa) indica que el modelo completo proporciona un ajuste sustancialmente mejor, lo que sugiere una variación residual entre lotes y una posible subintegración.

Definición de la puntuación de identidad celular (CIS)

La CIS cuantifica la medida en que se conservan las identidades celulares biológicamente distintas en la incrustación integrada. Intuitivamente, una dimensión de incrustación que captura las identidades celulares debe exhibir una distribución multimodal, cuyos modos (es decir, componentes de mezcla) corresponden a distintas identidades celulares [^[66]]. Por el contrario, si una dimensión de incrustación exhibe una distribución unimodal, no revela identidades celulares.

Por lo tanto, para cada dimensión de incrustación , comparamos un modelo de mezcla gaussiana con un modelo gaussiano único, ambos ajustados dentro de cada lote. Sea denota la función de verosimilitud del modelo de mezcla gaussiana y denota la función de verosimilitud del modelo gaussiano único evaluado en el lote y la dimensión . Definimos la CIS por dimensión como

La CIS general se obtiene promediando sobre todas las dimensiones de la incrustación integrada:

Una CIS más alta indica que el modelo de mezcla gaussiana explica bien los valores de incrustación en todas las dimensiones, lo que sugiere una mayor preservación de la estructura de identidad celular, mientras que una CIS baja sugiere la pérdida de identidades celulares y una posible sobreintegración.

Definición de la puntuación de integración

La puntuación de integración a nivel de conjunto de datos se define como

donde BAS penaliza la subintegración y CIS penaliza la sobreintegración. Por lo tanto, la puntuación de integración equilibra el equilibrio entre la mezcla de lotes y la preservación de la identidad celular, proporcionando un criterio bien definido y sin anotaciones. Al maximizar la puntuación de integración, IntegrateRigor optimiza la configuración de hiperparámetros de un método de integración dado en un conjunto de datos dado.

Hiperparámetro para los métodos de integración

En los seis conjuntos de datos de integración que analizamos en la subsección "IntegrateRigor mejora constantemente la integración en diversos conjuntos de datos y métodos", evaluamos los siguientes hiperparámetros y valores para cada método de integración.

Para Harmony, realizamos una búsqueda exhaustiva de dos hiperparámetros: , que controla la intensidad de la penalización del agrupamiento de diversidad, y , que controla el número de grupos utilizados en el agrupamiento suave. Esta búsqueda exhaustiva resultó en ajustes de parámetros para cada conjunto de datos. Para Seurat-RPCA, exploramos el hiperparámetro , que determina el número de vecinos considerados al ponderar los anclajes y está asociado con el grado de eliminación de la variación entre lotes. Para scVI, exploramos el hiperparámetro , que especifica el número de nodos en cada capa oculta de la red neuronal. Para FastMNN, exploramos el hiperparámetro , que especifica el número de vecinos más cercanos utilizados para el emparejamiento de vecinos más cercanos entre lotes.

Generalización a múltiples lotes

Supongamos que existen lotes de consulta, además de un lote de referencia. Sea

denota la matriz de conteo de genes para el lote , y sea

denota la matriz de incrustación de baja dimensión correspondiente. Aquí, se refiere al lote de referencia, mientras que se refieren a los lotes de consulta.

Para cada lote de consulta , al compararlo con el lote de referencia, se obtiene una diversidad intra-gen (BSS) y una diversidad intra-dimensional (BAS), denominadas BSSb y BASb, respectivamente. Luego, definimos la diversidad intra-gen general para cada gen y la diversidad intra-dimensional general para cada dimensión promediando entre los lotes de consulta:

El índice de similitud cruzada (CIS) se generaliza de forma natural al entorno de múltiples lotes promediando la mejora de la verosimilitud por lote en todos los lotes, incluido el lote de referencia:

Finalmente, la puntuación de integración a nivel de conjunto de datos sigue estando definida como

Validación de los modelos estadísticos subyacentes a IntegrateRigor

Para validar los supuestos de modelado subyacentes a IntegrateRigor, examinamos la bondad de ajuste de los dos modelos probabilísticos: el modelo de mezcla de Poisson negativos (NBMM) para los conteos de expresión a nivel de gen y el modelo de mezcla gaussiana (GMM) para cada dimensión de la incrustación integrada. Del conjunto de datos de enfermedad hepática asociada con insuficiencia intestinal (IFALD) [^[40]], seleccionamos tres donantes que comprenden 32.397 genes y 15.895 células. Para la validación a nivel de incrustación, aplicamos cinco métodos de integración (Seurat-RPCA, Harmony, scVI, FastMNN y LIGER) con la configuración predeterminada para integrar los tres lotes y evaluamos el ajuste del GMM a la incrustación integrada resultante (a nivel de dimensión) para cada método. Para la validación a nivel de gen, aplicamos el NBMM a un lote antes de la integración, como ejemplo ilustrativo.

Cuantificamos la bondad de ajuste utilizando la estadística de Kolmogorov-Smirnov (KS) [^[67]], definida como la diferencia absoluta máxima entre la función de distribución empírica y la función de distribución acumulada ajustada, para medir la magnitud de la discrepancia de distribución. No utilizamos el valor p para evaluar cualitativamente la bondad de ajuste, porque es bien sabido que para tamaños de muestra grandes, como los de los conjuntos de datos de una sola célula a gran escala, el valor p de KS tiende a volverse infinitesimalmente pequeño, de modo que incluso desviaciones menores y prácticamente insignificantes de la distribución ajustada pueden conducir al rechazo del modelo [^[68]].

Para el modelado de los conteos de expresión por gen, los histogramas observados frente a los ajustados muestran que el NBMM se ajusta estrechamente a las distribuciones de conteo para los genes representativos en IFALD. Esto es consistente con la amplia evidencia previa de que la distribución NB proporciona un modelo eficaz para los datos de conteo de expresión génica sobre dispersos en poblaciones celulares homogéneas (células con la misma identidad subyacente) [^[64]]. En los ejemplos mostrados, la mayoría de los genes tienen valores de la estadística KS muy pequeños, a menudo cercanos a cero y generalmente inferiores a 0,02, con solo unos pocos genes que tienen valores alrededor de 0,03 (Figura S2). Este patrón se extiende más allá de los ejemplos mostrados. En 2000 HVG, la distribución de los valores de la estadística KS se concentra cerca de cero, con solo una cola superior modesta y un número limitado de valores atípicos (Figura S1a). En general, estas observaciones respaldan el NBMM como un modelo razonable para el modelado de la distribución a nivel de gen.

Para el modelado de los valores de incrustación de integración por dimensión, el modelado gaussiano está motivado tanto empíricamente como teóricamente: los resultados clásicos sobre la reducción de dimensión lineal implican que muchas incrustaciones de baja dimensión de datos de alta dimensión son aproximadamente gaussianas [^[69]]. Para evaluar este supuesto de GMM empíricamente, ajustamos un GMM a cada una de las primeras cinco dimensiones de la incrustación integrada producida por cada uno de los cinco métodos de integración y evaluamos la bondad de ajuste utilizando la estadística KS (Figura S3). Para scVI, Harmony, Seurat-RPCA y FastMNN, los valores de la estadística KS son consistentemente pequeños en las primeras cinco dimensiones, generalmente del orden de 0,004 a 0,016, lo que indica que sus valores de incrustación a nivel de dimensión están bien aproximados por el GMM. El resumen del diagrama de caja entre dimensiones muestra la misma tendencia, con los valores de la estadística KS concentrados cerca de cero (Figura S1b).

En contraste, la incrustación integrada de LIGER exhibe valores de KS sustancialmente mayores, alrededor de 0,223 a 0,332 en las primeras cinco dimensiones, lo que indica una bondad de ajuste más pobre para el GMM en relación con las incrustaciones integradas de los otros cuatro métodos. Este resultado es consistente con el hecho de que LIGER se basa en la factorización de matrices no negativas, que tiende a producir incrustaciones con distribuciones no negativas, sesgadas y menos gaussianas. Como resultado, excluimos LIGER de los análisis de ajuste de hiperparámetros y selección de métodos que se basan en nuestra puntuación de integración, que asume un GMM.

En conjunto, estos resultados respaldan los dos supuestos de modelado centrales en IntegrateRigor: el NBMM proporciona un buen ajuste a los conteos de expresión por gen, y el GMM proporciona una aproximación razonable para los valores por dimensión en la incrustación integrada. En la práctica, estos supuestos de distribución deben verificarse antes de aplicar IntegrateRigor, para garantizar que el NBMM y el GMM sigan siendo apropiados y que los resultados de IntegrateRigor sean válidos para un conjunto de datos determinado.

Estabilidad y convergencia

El NBMM se ajustó utilizando un algoritmo EM personalizado, mientras que el GMM se ajustó utilizando la función Mclust en el paquete R mclust. Los parámetros clave para IntegrateRigor son el número de componentes de la mezcla, es decir, , , , y . Sin embargo, para el ajuste del modelo de mezcla, siempre que el número de componentes de la mezcla especificado no sea menor que el número real de componentes de la mezcla subyacentes, se espera que nuestras puntuaciones (BSS, BAS y CIS) sigan siendo estables en diferentes elecciones de , , , y [^[70]]. De forma predeterminada, el número de componentes de la mezcla se estableció en 5 para el NBMM y el GMM. Aunque el número total de identidades celulares (por ejemplo, tipos de células) en un conjunto de datos suele ser mayor que 5, la distribución de un solo gen o una sola dimensión de incrustación normalmente no exhibe muchos modos distintos, porque cada identidad celular se define conjuntamente por múltiples genes y múltiples dimensiones de incrustación. Empíricamente, observamos que el número de modos por gen o dimensión de incrustación es modesto, y establecer el número de componentes de la mezcla en 5 es suficiente en la práctica.

Evaluamos además la estabilidad de IntegrateRigor con respecto al número especificado de componentes de la mezcla. Específicamente, examinamos la estabilidad del BSS, los BSG seleccionados y la puntuación de integración en diferentes números de componentes de la mezcla: . En relación con el ajuste del NBMM a nivel de gen utilizando nuestra implementación de algoritmo EM personalizada, variamos el número de iteraciones EM de 10, 20 a 50, con 10 iteraciones utilizadas como valor predeterminado. En estas configuraciones, los valores de BSS resultantes, cuyo cálculo se basa en el ajuste del NBMM, estuvieron consistentemente altamente correlacionados con los obtenidos en la configuración predeterminada (con la correlación de Spearman cercana a 1), y los BSG seleccionados mostraron una superposición sustancial con el conjunto de BSG predeterminado (con el índice de Jaccard para la similitud de conjuntos cercano a 1). En relación con el ajuste del GMM a nivel de incrustación con diferentes números de componentes de la mezcla, la puntuación de integración seleccionó consistentemente el mismo hiperparámetro óptimo. Estos resultados sugieren que IntegrateRigor es un procedimiento estable. Se proporcionan más detalles en las Figuras suplementarias S11 a S14.

Métricas de evaluación basadas en la anotación

Índice local inverso de Simpson consciente del tipo de célula (CiLISI)

El índice local inverso de Simpson de integración (iLISI) es una métrica basada en la vecindad local para la mezcla de lotes. Para cada célula, iLISI mide la diversidad de células de diferentes lotes entre sus vecinos más cercanos en la incrustación integrada. Un valor de iLISI más alto indica que la vecindad local de la célula contiene una mezcla más equilibrada de lotes, mientras que un valor más bajo indica que la vecindad está dominada por un subconjunto de lotes.

Para evaluar la mezcla de lotes teniendo en cuenta la preservación de la identidad celular, utilizamos una versión consciente del tipo de célula de iLISI [^[13], ^[23], ^[37]], denominada CiLISI. Para cada tipo de célula, restringimos la incrustación integrada a las células anotadas como ese tipo de célula y calculamos iLISI utilizando las etiquetas de lote dentro de este subconjunto de células. Por lo tanto, el valor de iLISI de cada célula refleja la diversidad del lote entre sus vecinos más cercanos del mismo tipo de célula, evitando la confusión de las diferencias de los distintos tipos de células.

Luego, promediamos los valores a nivel de célula dentro de cada tipo de célula para obtener un iLISI específico del tipo de célula. Estos iLISIs específicos del tipo de célula se promediaron además entre los tipos de células, ponderados por las proporciones de los tipos de células, para obtener la puntuación final de CiLISI. Un CiLISI más alto indica una mejor mezcla de lotes dentro del mismo tipo de célula.

Índice de Rand ajustado (ARI)

Para evaluar la preservación de la estructura de la identidad celular, utilizamos el ARI entre los grupos no supervisados obtenidos de la incrustación integrada y las etiquetas de identidad celular conocidas [^[13], ^[37]]. Para cada incrustación integrada, realizamos el agrupamiento de Louvain de Seurat en un rango de parámetros de resolución de 0,1 a 1,0, en incrementos de 0,05. En cada resolución, calculamos el ARI entre las asignaciones de grupos resultantes y las etiquetas de identidad celular conocidas. El ARI máximo en todas las resoluciones probadas se tomó como la puntuación final de ARI para esa incrustación. Este procedimiento evita la dependencia de una única resolución de agrupamiento y captura la medida en que la incrustación integrada refleja la estructura de la identidad celular.

Índice local inverso del tipo de célula (cLISI)

Para cada célula, cLISI mide la diversidad de las etiquetas del tipo de célula entre sus vecinos más cercanos. Un valor de cLISI más bajo indica que la vecindad local está dominada por células del mismo tipo de célula, lo que refleja una mejor preservación de la estructura de la identidad celular. Un valor de cLISI más alto indica una mezcla local más fuerte de diferentes tipos de células, lo que puede sugerir una pérdida de la separación del tipo de célula o una sobreintegración.

Primero, calculamos LISI utilizando las etiquetas del tipo de célula en todas las células, sin condicionar en el lote, y nos referimos a estas puntuaciones como valores de cLISI a nivel de célula. El valor de cLISI de cada célula refleja la diversidad local de las etiquetas del tipo de célula entre sus células vecinas, lo que permite que cLISI detecte la mezcla de diferentes tipos de células entre los lotes y, por lo tanto, capture la posible sobreintegración. Luego, promediamos los valores de cLISI a nivel de célula en todas las células para obtener el valor de cLISI a nivel de conjunto de datos.

Selección de genes estables por lote (BSG) para la integración de datos

Modelado estadístico de la expresión génica en los lotes de referencia y consulta

Nos centramos en un solo gen y suprimimos el índice de gen para simplificar la notación. Sea denota los conteos de expresión en el lote , con los tamaños de biblioteca correspondientes . Aquí, se refiere al lote de referencia, y se refiere al lote de consulta.

Para identificar los genes que conservan la identidad celular entre diferentes lotes, modelamos la expresión génica utilizando un modelo de mezcla de distribución binomial negativa (DBN), con los componentes de la mezcla que representan distintas identidades celulares. Este modelo se basa en la amplia evidencia previa de que la distribución DBN modela eficazmente los datos de recuento de expresión génica con sobredispersión dentro de una identidad celular homogénea [^[64]]. La validación de este supuesto del modelo se proporciona en la Subsección Métodos “Validación de los modelos estadísticos subyacentes a IntegrateRigor”. Específicamente, asumimos que para cualquier gen dado

y

donde . Aquí, denota la proporción de mezcla del componente en el lote , con , y y denotan los parámetros de media y dispersión correspondientes para ese componente DBN, respectivamente. Asumimos que estos parámetros específicos del componente son distintos entre los diferentes componentes de la mezcla; es decir,

Intuitivamente, si un gen refleja la identidad celular en lugar de la variación entre lotes, su patrón de expresión en las diferentes identidades celulares debe conservarse entre los lotes, incluso cuando las proporciones de esas identidades celulares difieren. En este caso, se asume que la distribución de la expresión génica en el lote de consulta comparte un subconjunto de los componentes de la mezcla en el lote de referencia, al tiempo que solo se permite que las proporciones de mezcla difieran.

Específicamente, en ausencia de variación entre lotes, los componentes de la mezcla en el modelo del lote de consulta (2) son un subconjunto de los componentes de la mezcla en el modelo del lote de referencia (1). Entonces, la distribución de la expresión génica del lote de consulta se puede volver a parametrizar utilizando , al tiempo que solo se permite que las proporciones de los componentes varíen:

a la que nos referimos como el modelo restringido por referencia. Aquí, denota la proporción de mezcla del componente de referencia en el modelo restringido por referencia para el lote de consulta. Lo distinguimos de , que se refiere a la proporción de mezcla del componente de consulta en el modelo (2), y los componentes de referencia y los componentes de consulta pueden no compartir los mismos índices.

Por el contrario, la variación entre lotes induce cambios en los parámetros de los componentes, de modo que . En consecuencia, la distribución de la consulta ya no puede explicarse adecuadamente solo con la estructura de la mezcla de referencia y, en cambio, requiere el modelo completo
(2). Por lo tanto, definimos el BSS para el gen en función de qué tan bien el lote de consulta puede explicarse mediante el modelo restringido por referencia en relación con el modelo completo.

Modelado estadístico para la expresión génica en los lotes de referencia y consulta

Nos centramos en un solo gen y suprimimos el índice del gen para simplificar la notación. Sea denota los recuentos de expresión en el lote , con los tamaños de biblioteca correspondientes . Aquí, se refiere al lote de referencia y se refiere al lote de consulta.

Para identificar los genes que conservan la identidad celular entre los lotes, modelamos la expresión génica utilizando un modelo de mezcla de distribución binomial negativa (DBN), con los componentes de la mezcla que representan distintas identidades celulares. Este modelo se basa en la amplia evidencia previa de que la distribución DBN modela eficazmente los datos de recuento de expresión génica con sobredispersión dentro de una identidad celular homogénea [^[64]]. La validación de este supuesto del modelo se proporciona en la Subsección Métodos “Validación de los modelos estadísticos subyacentes a IntegrateRigor”. Específicamente, asumimos que para cualquier gen dado

y

donde . Aquí, denota la proporción de mezcla del componente en el lote , con , y y denotan los parámetros de media y dispersión correspondientes para ese componente DBN, respectivamente. Asumimos que estos parámetros específicos del componente son distintos entre los diferentes componentes de la mezcla; es decir,

Intuitivamente, si un gen refleja la identidad celular en lugar de la variación entre lotes, su patrón de expresión en las diferentes identidades celulares debe conservarse entre los lotes, incluso cuando las proporciones de esas identidades celulares difieren. En este caso, se asume que la distribución de la expresión génica en el lote de consulta comparte un subconjunto de los componentes de la mezcla en el lote de referencia, al tiempo que solo se permite que las proporciones de mezcla difieran.

Específicamente, en ausencia de variación entre lotes, los componentes de la mezcla en el modelo del lote de consulta (2) son un subconjunto de los componentes de la mezcla en el modelo del lote de referencia (1). Entonces, la distribución de la expresión génica del lote de consulta se puede volver a parametrizar utilizando , al tiempo que solo se permite que las proporciones de los componentes varíen:

a la que nos referimos como el modelo restringido por referencia. Aquí, denota la proporción de mezcla del componente de referencia en el modelo restringido por referencia para el lote de consulta. Lo distinguimos de , que se refiere a la proporción de mezcla del componente de consulta en el modelo (2), y los componentes de referencia y los componentes de consulta pueden no compartir los mismos índices.

Por el contrario, la variación entre lotes induce cambios en los parámetros de los componentes, de modo que . En consecuencia, la distribución de la consulta ya no puede explicarse adecuadamente solo con la estructura de la mezcla de referencia y, en cambio, requiere el modelo completo
(2). Por lo tanto, definimos el BSS para el gen en función de qué tan bien el lote de consulta puede explicarse mediante el modelo restringido por referencia en relación con el modelo completo.

Definición de la puntuación de estabilidad del lote (BSS)

Para distinguir estos dos escenarios, comparamos dos modelos, el modelo restringido por referencia (3) y el modelo completo (2), para el lote de consulta, ambos ajustados utilizando el algoritmo de Expectativa-Maximización (EM) [^[65]]. Bajo el modelo restringido por referencia (3), los parámetros de los componentes se estiman primero a partir del lote de referencia mediante EM y luego se tratan como fijos al evaluar el lote de consulta, con solo las proporciones de mezcla reestimadas por EM. Bajo el modelo completo (2), todos los parámetros se estiman directamente a partir del lote de consulta utilizando EM.

Sea denota la función de verosimilitud del modelo restringido por referencia y de la función de verosimilitud del modelo completo, ambos evaluados en función de las estimaciones de los parámetros del lote de consulta. Entonces, definimos la puntuación de estabilidad del lote (BSS) como

que toma un valor negativo ya que es estrictamente menor que ya que el espacio de parámetros del modelo restringido por referencia es un subespacio estricto del espacio de parámetros del modelo completo.

Una BSS alta (menos negativa) indica que el modelo restringido por referencia explica bien los datos de la consulta, lo que sugiere que el patrón de expresión del gen se transfiere entre los lotes y, por lo tanto, es estable entre lotes. Por el contrario, una BSS baja (más negativa) indica que el modelo completo proporciona un ajuste sustancialmente mejor, lo que sugiere una variación entre lotes. Los genes estables entre lotes se seleccionan en función de la distribución empírica de los valores de BSS utilizando un umbral adaptativo de datos determinado por un criterio de codo [^[41]].

Puntuación de integración para un conjunto de datos y un método de integración dados. Modelado estadístico para la incrustación integrada

Sea y denotan las incrustaciones de baja dimensión de los lotes de referencia y consulta de una incrustación integrada producida por un método de integración dado, donde denota la j-ésima dimensión de la incrustación para el lote , donde denota el lote de referencia y denota el lote de consulta.

Para cada dimensión de incrustación , asumimos que

y

donde . Aquí, denota la proporción de mezcla del componente en la dimensión de incrustación para el lote , con , y y denotan los parámetros de media y varianza correspondientes para ese componente gaussiano, respectivamente. Dada nuestra observación de que el modelo de mezcla gaussiana refleja distintas identidades celulares en algunas dimensiones de la incrustación integrada (ver la Subsección Métodos “Validación de los modelos estadísticos subyacentes a IntegrateRigor”), asumimos que los componentes de la mezcla en este modelo corresponden a identidades celulares y que los parámetros específicos del componente son distintos; es decir,

Definición de la puntuación de alineación del lote (BAS)

La BAS cuantifica la mezcla de lotes en el espacio de incrustación integrada. Para cada dimensión de incrustación , siguiendo la misma lógica que la definición de BSS, comparamos un modelo restringido por referencia,

con el modelo completo (4) para las incrustaciones del lote de consulta en la misma dimensión. Bajo el modelo restringido por referencia, los parámetros de los componentes se estiman a partir del lote de referencia y luego se tratan como fijos al evaluar el lote de consulta, con solo las proporciones de mezcla reestimadas. Bajo el modelo completo, todos los parámetros se estiman directamente a partir del lote de consulta.

Sea denota la función de verosimilitud del modelo restringido por referencia y denota la función de verosimilitud del modelo completo, ambos evaluados en las incrustaciones de la consulta en la dimensión . Definimos la BAS por dimensión como

La BAS general se obtiene promediando sobre todas las dimensiones de incrustación:

Una BAS más alta (menos negativa) indica que el modelo restringido por referencia explica bien las incrustaciones de la consulta, lo que sugiere una mejor alineación entre los lotes y una eliminación más eficaz de la variación entre lotes. Por el contrario, una BAS más baja (más negativa) indica que el modelo completo proporciona un ajuste sustancialmente mejor, lo que sugiere una variación residual entre lotes y una posible subintegración.

Definición de la puntuación de identidad celular (CIS)

La CIS cuantifica la medida en que se conservan las identidades celulares biológicamente distintas en la incrustación integrada. Intuitivamente, una dimensión de incrustación que captura las identidades celulares debe exhibir una distribución multimodal, cuyos modos (es decir, componentes de la mezcla) corresponden a distintas identidades celulares [^[66]]. Por el contrario, si una dimensión de incrustación exhibe una distribución unimodal, no revela identidades celulares.

Por lo tanto, para cada dimensión de incrustación , comparamos un modelo de mezcla gaussiana con un modelo gaussiano único, ambos ajustados dentro de cada lote. Sea denota la función de verosimilitud del modelo de mezcla gaussiana y denota la función de verosimilitud del modelo gaussiano único evaluado en el lote y la dimensión . Definimos la CIS por dimensión como

La CIS general se obtiene promediando sobre todas las dimensiones de la incrustación integrada:

Una CIS más alta indica que el modelo de mezcla gaussiana explica bien los valores de incrustación en todas las dimensiones, lo que sugiere una mayor preservación de la estructura de identidad celular, mientras que una CIS baja sugiere la pérdida de identidades celulares y una posible sobreintegración.

Definición de la puntuación de integración

La puntuación de integración a nivel de conjunto de datos se define como

donde BAS penaliza la subintegración y CIS penaliza la sobreintegración. Por lo tanto, la puntuación de integración equilibra el equilibrio entre la mezcla de lotes y la preservación de la identidad celular, proporcionando un criterio bien definido y sin anotaciones. Al maximizar la puntuación de integración, IntegrateRigor optimiza la configuración de hiperparámetros de un método de integración dado en un conjunto de datos dado.

Hiperparámetro para los métodos de integración

En los seis conjuntos de datos de integración que analizamos en la Subsección “IntegrateRigor mejora constantemente la integración en diversos conjuntos de datos y métodos”, evaluamos los siguientes hiperparámetros y valores para cada método de integración.

Para Harmony, realizamos una búsqueda en cuadrícula sobre dos hiperparámetros: , que controla la fuerza de la penalización de agrupamiento de diversidad y , que controla el número de grupos utilizados en el agrupamiento suave. Esta búsqueda en cuadrícula dio como resultado configuraciones de parámetros para cada conjunto de datos. Para Seurat-RPCA, buscamos en el hiperparámetro, , que determina el número de vecinos considerados al ponderar los anclajes y está asociado con el grado de eliminación de la variación entre lotes. Para scVI, buscamos en el hiperparámetro, , que especifica el número de nodos en cada capa oculta de la red neuronal. Para FastMNN, buscamos en el hiperparámetro, , que especifica el número de vecinos más cercanos utilizados para el emparejamiento de vecinos más cercanos entre lotes.

Generalización a múltiples lotes

Supongamos que hay lotes de consulta además de un lote de referencia. Sea

denota la matriz de recuento de genes para el lote , y sea

denota la matriz de incrustación de baja dimensión correspondiente. Aquí, se refiere al lote de referencia, mientras que se refieren a los lotes de consulta.

Para cada lote de consultas, la comparación con el lote de referencia produce una medida BSS a nivel de gen y una medida BAS a nivel de dimensión, denominadas respectivamente BSSb y BASb. A continuación, definimos la medida BSS general para cada gen y la medida BAS general para cada dimensión promediando los valores de los diferentes lotes de consultas:

El CIS se generaliza de forma natural al entorno de múltiples lotes promediando la mejora de la función de verosimilitud por lote en todos los lotes, incluido el lote de referencia:

Finalmente, la puntuación de integración a nivel de conjunto de datos se define como:

Validación de los modelos estadísticos subyacentes a IntegrateRigor

Para validar los supuestos de modelado subyacentes a IntegrateRigor, examinamos la bondad de ajuste de los dos modelos probabilísticos: el modelo de mezcla NB (NBMM) para los recuentos de expresión a nivel de gen y el modelo de mezcla gaussiana (GMM) para cada dimensión del espacio de incrustación integrado. Del conjunto de datos de enfermedad hepática asociada con insuficiencia intestinal (IFALD) [^[40]], seleccionamos tres donantes que comprenden 32.397 genes y 15.895 células. Para la validación a nivel de incrustación, aplicamos cinco métodos de integración (Seurat-RPCA, Harmony, scVI, FastMNN y LIGER) con la configuración predeterminada para integrar los tres lotes y evaluamos el ajuste del GMM al espacio de incrustación integrado resultante (a nivel de dimensión) para cada método. Para la validación a nivel de gen, aplicamos el NBMM a un lote antes de la integración, como ejemplo ilustrativo.

Cuantificamos la bondad de ajuste utilizando la estadística de Kolmogorov-Smirnov (KS) [^[67]], definida como la diferencia absoluta máxima entre la función de distribución empírica y la función de distribución acumulada ajustada, para medir la magnitud de la discrepancia de la distribución. No utilizamos el valor p para evaluar cualitativamente la bondad de ajuste, porque es bien sabido que para tamaños de muestra grandes, como los de los conjuntos de datos de una sola célula a gran escala, el valor p de KS tiende a volverse infinitesimalmente pequeño, de modo que incluso desviaciones menores y prácticamente insignificantes de la distribución ajustada pueden conducir al rechazo del modelo [^[68]].

Para el modelado de los recuentos de expresión por gen, los histogramas observados frente a los ajustados muestran que el NBMM se ajusta estrechamente a las distribuciones de recuento para los genes representativos en IFALD. Esto es consistente con la amplia evidencia previa de que la distribución NB proporciona un modelo eficaz para los datos de recuento de expresión génica sobredispersos en poblaciones celulares homogéneas (células con la misma identidad subyacente) [^[64]]. En los ejemplos mostrados, la mayoría de los genes tienen valores de la estadística KS muy pequeños, a menudo cercanos a cero y generalmente inferiores a 0,02, con solo unos pocos genes que tienen valores alrededor de 0,03 (Figura S2). Este patrón se extiende más allá de los ejemplos mostrados. En 2000 HVG, la distribución de los valores de la estadística KS se concentra cerca de cero, con solo una cola superior modesta y un número limitado de valores atípicos (Figura S1a). En general, estas observaciones respaldan el NBMM como un modelo razonable para el modelado de la distribución a nivel de gen.

Para el modelado de los valores de incrustación de integración por dimensión, el modelado gaussiano está motivado tanto empíricamente como teóricamente: los resultados clásicos sobre la reducción de la dimensión lineal implican que muchas incrustaciones de baja dimensión de datos de alta dimensión son aproximadamente gaussianas [^[69]]. Para evaluar este supuesto del GMM empíricamente, ajustamos un GMM a cada una de las primeras cinco dimensiones del espacio de incrustación integrado producido por cada uno de los cinco métodos de integración y evaluamos la bondad de ajuste utilizando la estadística KS (Figura S3). Para scVI, Harmony, Seurat-RPCA y FastMNN, los valores de la estadística KS son consistentemente pequeños en las primeras cinco dimensiones, generalmente del orden de 0,004 a 0,016, lo que indica que sus valores de incrustación a nivel de dimensión están bien aproximados por el GMM. El diagrama de caja resumen en las dimensiones muestra la misma tendencia, con los valores de la estadística KS concentrados cerca de cero (Figura S1b).

En contraste, el espacio de incrustación integrado de LIGER exhibe valores de KS sustancialmente mayores, alrededor de 0,223 a 0,332 en las primeras cinco dimensiones, lo que indica un ajuste de GMM más deficiente en relación con los espacios de incrustación integrados de los otros cuatro métodos. Este resultado es consistente con el hecho de que LIGER se basa en la factorización de matrices no negativas, lo que tiende a producir incrustaciones con distribuciones no negativas, sesgadas y menos gaussianas. Como resultado, excluimos LIGER de los análisis de ajuste de hiperparámetros y selección de métodos que se basan en nuestra puntuación de integración, que asume un GMM.

En conjunto, estos resultados respaldan los dos supuestos de modelado centrales en IntegrateRigor: el NBMM proporciona un buen ajuste a los recuentos de expresión por gen, y el GMM proporciona una aproximación razonable para los valores por dimensión en el espacio de incrustación integrado. En la práctica, estos supuestos de distribución deben verificarse antes de aplicar IntegrateRigor, para garantizar que el NBMM y el GMM sigan siendo apropiados y que los resultados de IntegrateRigor sean válidos para un conjunto de datos determinado.

Estabilidad y convergencia

El NBMM se ajustó utilizando un algoritmo EM personalizado, mientras que el GMM se ajustó utilizando la función Mclust en el paquete R mclust. Los parámetros clave para IntegrateRigor son el número de componentes de la mezcla, es decir, , , , y . Sin embargo, para el ajuste del modelo de mezcla, siempre que el número de componentes de la mezcla especificado no sea menor que el número real de componentes de la mezcla subyacentes, se espera que nuestras puntuaciones (BSS, BAS y CIS) sigan siendo estables en diferentes elecciones de , , , y [^[70]]. De forma predeterminada, el número de componentes de la mezcla se estableció en 5 para el NBMM y el GMM. Aunque el número total de identidades celulares (por ejemplo, tipos de células) en un conjunto de datos suele ser mayor que 5, la distribución de un solo gen o una sola dimensión de incrustación normalmente no exhibe muchos modos distintos, porque cada identidad celular se define conjuntamente por múltiples genes y múltiples dimensiones de incrustación. Empíricamente, observamos que el número de modos por gen o dimensión de incrustación es modesto, y establecer el número de componentes de la mezcla en 5 es suficiente en la práctica.

Evaluamos además la estabilidad de IntegrateRigor con respecto al número especificado de componentes de la mezcla. Específicamente, examinamos la estabilidad de la BSS, los BSG seleccionados y la puntuación de integración en diferentes números de componentes de la mezcla: . En relación con el ajuste del NBMM a nivel de gen utilizando nuestra implementación personalizada del algoritmo EM, variamos el número de iteraciones EM de 10, 20 a 50, con 10 iteraciones utilizadas como valor predeterminado. En estas configuraciones, los valores de BSS resultantes, cuyo cálculo se basa en el ajuste del NBMM, estuvieron consistentemente altamente correlacionados con los obtenidos con la configuración predeterminada (con la correlación de Spearman cercana a 1), y los BSG seleccionados mostraron una superposición sustancial con el conjunto de BSG predeterminado (con el índice de Jaccard para la similitud de conjuntos cercano a 1). En relación con el ajuste del GMM a nivel de incrustación con diferentes números de componentes de la mezcla, la puntuación de integración seleccionó consistentemente el mismo hiperparámetro óptimo. Estos resultados sugieren que IntegrateRigor es un procedimiento estable. Se proporcionan más detalles en las Figuras Suplementarias S11 a S14.

Métricas de evaluación basadas en la anotación

Índice local inverso de Simpson consciente del tipo de célula (CiLISI)

El índice local inverso de Simpson de integración (iLISI) es una métrica basada en la vecindad local para la mezcla de lotes. Para cada célula, iLISI mide la diversidad de células de diferentes lotes entre sus vecinos más cercanos en el espacio de incrustación integrado. Un valor de iLISI más alto indica que la vecindad local de la célula contiene una mezcla más equilibrada de lotes, mientras que un valor más bajo indica que la vecindad está dominada por un subconjunto de lotes.

Para evaluar la mezcla de lotes teniendo en cuenta la preservación de la identidad celular, utilizamos una versión consciente del tipo de célula de iLISI [^[13], ^[23], ^[37]], denominada CiLISI. Para cada tipo de célula, restringimos el espacio de incrustación integrado a las células anotadas como ese tipo de célula y calculamos iLISI utilizando las etiquetas de lote dentro de este subconjunto de células. Por lo tanto, el valor de iLISI de cada célula refleja la diversidad del lote entre sus vecinos más cercanos del mismo tipo de célula, evitando la confusión de las diferencias de los distintos tipos de células.

Luego, promediamos los valores a nivel de célula dentro de cada tipo de célula para obtener un iLISI específico del tipo de célula. Estos iLISIs específicos del tipo de célula se promediaron además en los tipos de células, ponderados por las proporciones de los tipos de células, para obtener la puntuación CiLISI final. Un CiLISI más alto indica una mejor mezcla de lotes dentro del mismo tipo de célula.

Índice de Rand ajustado (ARI)

Para evaluar la preservación de la estructura de la identidad celular, utilizamos el ARI entre los grupos no supervisados obtenidos del espacio de incrustación integrado y las etiquetas de identidad celular conocidas [^[13], ^[37]]. Para cada espacio de incrustación integrado, realizamos el agrupamiento de Louvain de Seurat en un rango de parámetros de resolución de 0,1 a 1,0, en incrementos de 0,05. En cada resolución, calculamos el ARI entre las asignaciones de grupos resultantes y las etiquetas de identidad celular conocidas. El ARI máximo en todas las resoluciones probadas se tomó como la puntuación ARI final para ese espacio de incrustación. Este procedimiento evita la dependencia de una única resolución de agrupamiento y captura la medida en que el espacio de incrustación integrado refleja la estructura de la identidad celular.

Índice local inverso de tipo de célula (cLISI)

Para cada célula, cLISI mide la diversidad de las etiquetas de tipo de célula entre sus vecinos más cercanos. Un valor de cLISI más bajo indica que la vecindad local está dominada por células del mismo tipo de célula, lo que refleja una mejor preservación de la estructura de la identidad celular. Un valor de cLISI más alto indica una mezcla local más fuerte de diferentes tipos de células, lo que puede sugerir una pérdida de la separación del tipo de célula o una sobreintegración.

Primero, calculamos LISI utilizando las etiquetas de tipo de célula en todas las células, sin condicionar en el lote, y nos referimos a estas puntuaciones como valores de cLISI a nivel de célula. El valor de cLISI de cada célula refleja la diversidad local de las etiquetas de tipo de célula entre sus células vecinas, lo que permite que cLISI detecte la mezcla de diferentes tipos de células en los lotes y, por lo tanto, capture la posible sobreintegración. Luego, promediamos los valores de cLISI a nivel de célula en todas las células para obtener el valor de cLISI a nivel de conjunto de datos.

Índice local inverso de Simpson de integración consciente del tipo de célula (CiLISI)

El índice local inverso de Simpson de integración (iLISI) es una métrica basada en la vecindad local para la mezcla de lotes. Para cada célula, iLISI mide la diversidad de células de diferentes lotes entre sus vecinos más cercanos en el espacio de incrustación integrado. Un valor de iLISI más alto indica que la vecindad local de la célula contiene una mezcla más equilibrada de lotes, mientras que un valor más bajo indica que la vecindad está dominada por un subconjunto de lotes.

Para evaluar la mezcla de lotes teniendo en cuenta la preservación de la identidad celular, utilizamos una versión consciente del tipo de célula de iLISI [^[13], ^[23], ^[37]], denominada CiLISI. Para cada tipo de célula, restringimos el espacio de incrustación integrado a las células anotadas como ese tipo de célula y calculamos iLISI utilizando las etiquetas de lote dentro de este subconjunto de células. Por lo tanto, el valor de iLISI de cada célula refleja la diversidad del lote entre sus vecinos más cercanos del mismo tipo de célula, evitando la confusión de las diferencias de los distintos tipos de células.

Luego, promediamos los valores a nivel de célula dentro de cada tipo de célula para obtener un iLISI específico del tipo de célula. Estos iLISIs específicos del tipo de célula se promediaron además en los tipos de células, ponderados por las proporciones de los tipos de células, para obtener la puntuación CiLISI final. Un CiLISI más alto indica una mejor mezcla de lotes dentro del mismo tipo de célula.

Índice de Rand ajustado (ARI)

Para evaluar la preservación de la estructura de identidad celular, utilizamos el índice de similitud ajustado (ARI) entre los grupos obtenidos mediante el análisis de agrupamiento no supervisado a partir de la incrustación integrada y las etiquetas conocidas de identidad celular [^[13], ^[37]]. Para cada incrustación integrada, realizamos un análisis de agrupamiento de Louvain con Seurat en un rango de parámetros de resolución de 0,1 a 1,0, en incrementos de 0,05. En cada resolución, calculamos el ARI entre las asignaciones de grupos resultantes y las etiquetas conocidas de identidad celular. El ARI máximo de todas las resoluciones probadas se tomó como la puntuación final del ARI para esa incrustación. Este procedimiento evita la dependencia de una única resolución de agrupamiento y captura el grado en que la incrustación integrada refleja la estructura de identidad celular.

Índice local inverso de Simpson para tipos de células (cLISI)

Para cada célula, el cLISI mide la diversidad de las etiquetas de tipo de célula entre sus vecinos más cercanos. Un valor de cLISI más bajo indica que la vecindad local está dominada por células del mismo tipo de célula, lo que refleja una mejor preservación de la estructura de identidad celular. Un valor de cLISI más alto indica una mayor mezcla local de diferentes tipos de células, lo que puede sugerir una pérdida de separación de tipos de células o una sobreintegración.

Primero, calculamos el LISI utilizando las etiquetas de tipo de célula en todas las células, sin condicionar por lote, y nos referimos a estas puntuaciones como valores de cLISI a nivel de célula. El valor de cLISI de cada célula refleja la diversidad local de las etiquetas de tipo de célula entre sus células vecinas, lo que permite que el cLISI detecte la mezcla de diferentes tipos de células entre lotes y, por lo tanto, capture una posible sobreintegración. Luego, promediamos los valores de cLISI a nivel de célula en todas las células para obtener el valor de cLISI a nivel de conjunto de datos.