El objetivo del método propuesto es identificar los tipos de células presentes en una mezcla de células individuales. La entrada del método es la matriz de expresión génica de una sola célula(Mgene × cell) en la que las filas representan los genes y las columnas representan las células. A continuación, proporcionamos más detalles sobre los datos de entrada y los diferentes pasos del marco propuesto. El enfoque general se muestra en la Fig. 1.

El flujo de trabajo general del método propuesto. Dada la matriz de expresión génica de una sola célula, el módulo (A) elimina los genes que no se expresan en ninguna célula. Usando la matriz resultante, el módulo (B) calcula la distancia euclidiana entre las células. La salida de este módulo es una matriz de distancia en la que las filas y columnas son las celdas (Dcell×cell). El módulo (C) reduce la dimensionalidad de la matriz de distancias utilizando la técnica de incrustación de vecinos estocásticos distribuidos en t (t-SNE). En este módulo, se emplea un método de silueta media para elegir el número óptimo de clústeres k. Finalmente, en el módulo (D), la matriz de distancia de dimensión inferior y el número óptimo de clústeres k obtenidos del módulo (C) se utilizan como datos de entrada para identificar la agrupación de celdas más estable. La Figura 2 muestra los detalles del módulo D.

Fuente de datos

Los ocho conjuntos de datos scRNA-seq disponibles públicamente, así como los cinco conjuntos de datos de simulación que utilizamos en nuestro análisis, se incluyen en los Materiales Complementarios. Entre los ocho conjuntos de datos reales, todos menos tres (Klein51, Patel52, Treutlein53) se consideran «estándar dorado», ya que las etiquetas de las celdas se conocen de manera definitiva. Patel52 y Treutlein53 se denominan «estándar de plata» por Kiselev et al.28 ya que sus etiquetas celulares se determinan en base a los métodos computacionales y el conocimiento de los autores de la biología subyacente.

Obtuvimos los datos procesados del sitio web de Hemberg lab (https://hemberg-lab.github.io/scRNA.seq.datasets). Hemberg et al.54 utilice el bioconductor de un solo elemento S4 clase 55 para almacenar los datos, y el scater package56 para el control de calidad y el trazado. Los datos normalizados se depositan como un objeto de experimento de célula única (.Archivo RData) y se accede a la información del tipo de celda en la columna cell_type1 de la ranura «colData» de este objeto. Los valores de expresión génica de las células se organizan como una matriz en la que las filas son células y las columnas son los genes. En nuestro análisis, se eliminan los genes (características) que no se expresan en ninguna célula. No filtramos ninguna célula en este análisis.

Filtrado de genes

Como se muestra en la Fig. 1A, eliminamos los genes / transcripciones que no se expresan en ninguna célula (el valor de expresión es cero en todas las células). Tales genes no pueden proporcionar información útil que pueda diferenciar entre tipos de células 57. El resultado de realizar el método de filtrado en la matriz de expresión génica de una sola célula (Mgene×cell) se utiliza como entrada para el segundo módulo del marco propuesto.

Medición de la diferencia entre las celdas

La distancia entre las celdas se calcula utilizando la métrica euclidiana (Fig. 1B). La salida de este paso es la distancia (diferencia) de la matriz Dcell×celular. Reducimos la dimensión de D mediante la incrustación de vecinos estocásticos distribuidos en t (t-SNE)34,58, la técnica de reducción/visualización de dimensionalidad no lineal (Fig. 1C). Nos referiremos a la salida como D’cell×l, donde 2 ≤ l ≤ cell. En este estudio, el número de dimensiones es 2.

Agrupación en clústeres

Identificación del número óptimo de clústeres

En esta sección se describe el tercer módulo del método propuesto (Fig. 1C). En este análisis, el t-SNE se aplica repetidamente (n = 50) en la matriz de distancia Dcell×cell para obtener la matriz de distancia reducida de dimensionalidad D’cell×l. Cada vez, el número óptimo de clústeres se calcula en función del método de silueta promedio utilizando la matriz de distancia reducida de dimensionalidad D’. Para encontrar el número óptimo de clústeres k, la agrupación de k-medias se aplica a la matriz D’ utilizando un valor de rango (predeterminado = 2:20), y se selecciona la k que maximiza la medida de silueta promedio. Finalmente, el promedio de los números seleccionados k a través de diferentes repeticiones (n = 50) (redondeado al entero más cercano) se considera como el número óptimo final de grupos.

La silueta evalúa la calidad de ese agrupamiento en función de lo bien que se agrupan sus puntos de datos. Se asigna una medida de silueta a cada punto de datos que representa cuán cerca está un punto de datos de su propio clúster en comparación con otros clústeres. Para cada punto de datos i, esta medida se calcula de la siguiente manera:

donde a(i) es la distancia media entre el punto i y el resto de puntos de datos dentro del mismo grupo. b (i) es la distancia media más pequeña de i a todos los puntos de cualquier otro grupo del que i no sea miembro. s (i) toma valores de -1 a 1, donde una puntuación positiva alta muestra que el punto de datos dado está bien agrupado (cerca de otros puntos en su propio clúster y lejos de los puntos en los otros clústeres). Por el contrario, una puntuación negativa alta muestra que el punto de datos está mal agrupado.

k-significa agrupación basada en el método de remuestreo

Esta sección describe los detalles del último módulo del método propuesto. Como se muestra en la Fig. 2, utilizando la matriz de distancia reducida de dimensionalidad D’ y el número elegido de clústeres k del paso anterior, identificamos la agrupación más estable generando diferentes soluciones de agrupación (clusteringi (i ∈ )) y medimos la estabilidad de cada solución de agrupación basada en un método de remuestreo. La medida de estabilidad asignada a cada agrupación en particular (denotada como clusteringi) representa la frecuencia con la que se conservan los clústeres k pertenecientes a esa agrupación cuando los datos de entrada (D’) se vuelven a muestrear varias veces. Los conjuntos de datos remuestreados se generan a partir de D’ reemplazando aleatoriamente el 5% de los puntos de datos (celdas) con ruido. Estos conjuntos de datos ruidosos se utilizan como entrada para el algoritmo k-means. Por lo tanto,se generan varios clusterings (clusteringi, j, j ∈ ) a partir de los datos remuestreados (versiones remuestreadas de clusteringi).

la Identificación de los más estables de la agrupación. En este análisis, dada la matriz de distancias de dimensión inferior D’cell×l y el número óptimo de clusters k, calculamos n diferentes clusterings (clustering1, …, clusteringn) usando el algoritmo de clustering de k-means. A continuación, se evalúa la estabilidad de cada agrupación sobre la base de un enfoque de remuestreo (cuadro gris). Se asigna una puntuación de estabilidad a cada agrupación en función de la frecuencia con la que se recuperan sus agrupaciones cuando se perturban los datos de entrada (se vuelve a muestrear). Como solución final se selecciona un agrupamiento con la puntuación de estabilidad máxima.

Para evaluar la estabilidad de cada cluster c en el clusteringi (clusteringi original), el cluster c se compara con todos los clusters en el clustering que se obtiene a partir de los datos de remuestreo (clusteringi,j) en función de la distancia Jaccard. El coeficiente de Jaccard59, una medida de similitud entre conjuntos, se utiliza para calcular la similitud entre dos grupos de la siguiente manera:

donde el término a y B son dos grupos, que consta de unos puntos de datos X = {x1, …, xN}.

Si la similitud de Jaccard entre el clúster c (del clúster original clusteringi) y el clúster más similar del clúster remuestreado es igual o superior a 0,75, ese clúster se considera estable (preservado). Por lo tanto, la estabilidad de cada clúster en clusteringi se calcula como el porcentaje de las veces que se conserva el clúster (coeficiente de Jaccard ≥ 0.75) a través de los m diferentes remuestreos.

Luego promediamos las medidas de estabilidad de los clústeres k pertenecientes a clusteringi, y lo consideramos como la medida de estabilidad general de clusteringi. Entre n soluciones de clustering diferentes (clusteringi (i ∈ )), seleccionamos la solución de clustering con la medida de estabilidad máxima como solución de clustering final.

La Figura 3 muestra el detalle del método de remuestreo que realizamos para calcular la medida de estabilidad para cada agrupación. Los clústeres que se obtienen aplicando k-media en el conjunto de datos remuestreado se comparan con los clústeres de los datos de entrada originales solo en función de los puntos sin ruido (los puntos de datos de ruido se excluyen cuando se comparan dos clústeres en función de la métrica de similitud Jaccard.

El remuestreo marco para el cálculo de la estabilidad de la medida para cada agrupación. La entrada incluye N puntos de datos X = {x1, …, xN}, el número de clústeres k, el número de remuestreos m y el agrupamiento C que se obtiene aplicando k-means en X. Este análisis genera datos de remuestreo de m reemplazando aleatoriamente el 5% de los puntos de datos con el ruido, y calcula los agrupamientos de remuestreo de m basados en el agrupamiento de k-means. Cada clúster c en C se compara con el clúster más similar en la agrupación de remuestreo, y se calcula el coeficiente de Jaccard entre los dos clústeres, mientras que se excluyen los puntos de ruido. El porcentaje de veces que los coeficientes de Jaccard son mayores que 0.75 se considera la medida de estabilidad para el clúster c. El promedio de medidas de estabilidad para todos los clústeres que pertenecen al clúster C se calcula y se considera la medida de estabilidad general para el clúster C.

Métodos de validación

Utilizamos 13 conjuntos de datos diferentes en los que se conocen los tipos de celdas (etiquetas). Para medir el nivel de similitud entre las etiquetas de referencia y las etiquetas inferidas que se obtienen por cada método de agrupamiento, utilizamos tres métricas diferentes: índice rand ajustado (ARI), información mutua ajustada (AMI) y medida en V como se explica a continuación.

Índice de rand ajustado

Dadas las etiquetas de las celdas, el Índice de Rand Ajustado (ARI)47 se utiliza para evaluar la similitud entre el agrupamiento inferido y el agrupamiento verdadero. ARI varía de 0, para una coincidencia deficiente (una agrupación aleatoria), a 1 para un acuerdo perfecto con la agrupación verdadera. Para un conjunto de n puntos de datos, la tabla de contingencia se construye en función del número compartido de puntos de datos entre dos clústeres. Supongamos que X = {X1, X2,…, XR} e Y = {Y1, Y2, …, YC} representan dos agrupaciones diferentes con agrupaciones R y C, respectivamente. La superposición entre X e Y se puede resumir como una tabla de contingencia MR×C = , donde i = 1…R, j = 1…C. Xi e Yj denotan un clúster en agrupamientos X e Y, e i y j se refieren al número de fila y al número de columna de la tabla de contingencia, respectivamente. El ARI se define de la siguiente manera:

Información mutua ajustada