Proceso completo de análisis unicelular: normalización de datos

A menudo se observan diferencias sistemáticas en la cobertura de secuenciación entre bibliotecas en los datos de secuenciación de ARN unicelular. A menudo son causadas por diferencias técnicas en la captura de ADNc o la eficiencia de la amplificación por PCR entre células, lo que se atribuye a la dificultad para lograr una preparación consistente de la biblioteca con un mínimo de material de partida. La normalización tiene como objetivo eliminar estas diferencias para que no interfieran con la comparación de perfiles de expresión entre células. Esto garantiza que cualquier heterogeneidad o expresión diferencial observada entre las poblaciones de células se deba a la biología y no a un sesgo técnico.

En este punto, es necesario tener en cuenta la diferencia entre normalización y corrección por lotes. La normalización se produce independientemente de la estructura del lote y sólo tiene en cuenta las desviaciones técnicas, mientras que la corrección de lote sólo se produce entre lotes y debe tener en cuenta tanto las desviaciones técnicas como las diferencias biológicas. El sesgo técnico tiende a afectar a los genes de manera similar o al menos en relación con sus propiedades biofísicas (por ejemplo, longitud, contenido de GC), mientras que las diferencias biológicas entre lotes pueden ser muy impredecibles. Como tal, las dos tareas implican suposiciones diferentes y, a menudo, métodos computacionales diferentes (aunque algunos paquetes de software están diseñados para realizar ambos pasos a la vez, como zinbwave). Por lo tanto, es importante evitar confundir datos "normalizados" y "corregidos por lotes", ya que estos datos a menudo representan cosas diferentes.

Nos centraremos principalmente en escalar la normalización, que es la estrategia de normalización más simple y más utilizada. Esto implica dividir todos los recuentos por celda por un factor de escala específico de la celda, a menudo llamado "factor de tamaño". La suposición aquí es que cualquier sesgo específico de una célula (por ejemplo, eficiencia de captura o amplificación) afectará a todos los genes por igual al escalar la media esperada para esa célula. El factor de tamaño de cada celda representa una estimación del sesgo relativo en esa celda, por lo que dividir su recuento por su factor de tamaño debería eliminar ese sesgo. Los "datos normalizados" resultantes se pueden utilizar para análisis posteriores, como agrupación y reducción de dimensionalidad. Para la demostración, utilizaremos un conjunto de datos del paquete de software scRNAseq.

La normalización del tamaño de la biblioteca es la estrategia más sencilla para realizar la normalización del escalado. Definimos el tamaño de la biblioteca como la suma de los recuentos de todos los genes en cada célula, suponiendo que su valor esperado aumenta con cualquier sesgo específico de la célula. Luego, con una constante de proporcionalidad definida, el "factor de tamaño del grupo" de cada celda es directamente proporcional al tamaño de su grupo, de modo que el factor de tamaño promedio para todas las celdas es igual a uno. Esta definición garantiza que los valores de expresión normalizados estén en la misma escala que los recuentos sin procesar, lo cual es útil para la interpretación, especialmente cuando se trabaja con datos transformados.

En los datos del cerebro de Zeisel, el factor de tamaño de la biblioteca varió hasta 10 veces entre las células. Esto es típico de la variación de la cobertura de datos de scRNA-seq.

Estrictamente hablando, el uso de factores de tamaño de biblioteca supone que no hay "desequilibrios" en los genes expresados ​​diferencialmente (DE) entre cualquier par de células. Es decir, cualquier regulación positiva en un subconjunto de genes puede compensarse con la misma magnitud de regulación negativa en un subconjunto diferente de genes. Esto garantiza que el tamaño de la biblioteca sea una estimación imparcial del sesgo relativo específico de la célula al evitar efectos de síntesis. Sin embargo, la DE equilibrada a menudo está ausente en las aplicaciones scRNA-seq, lo que significa que la normalización del tamaño de la biblioteca puede no producir valores de expresión normalizados precisos para el análisis posterior.

En la práctica, la precisión de la normalización no es una consideración importante en el análisis exploratorio de datos de scRNA-seq. El sesgo compositivo generalmente no afecta la separación de las poblaciones celulares, sino solo la magnitud de los cambios logarítmicos entre poblaciones celulares o tipos de células, en menor medida. Por lo tanto, la normalización del tamaño de la biblioteca suele ser suficiente en muchas aplicaciones donde el objetivo es identificar poblaciones de células y definir marcadores superiores para cada población de células.

Como se mencionó anteriormente, el sesgo compositivo se produce cuando existe alguna expresión diferencial desequilibrada entre muestras. Considere dos células, donde un solo gen X está regulado positivamente en la célula A en comparación con la célula B.

Esta regulación positiva significa que (i) se dedican más recursos de secuenciación a X en A y, por lo tanto, menos cobertura de otros genes no diferenciales cuando el tamaño total de la biblioteca por célula se determina experimentalmente (por ejemplo, debido a la cuantificación de la biblioteca), o (ii) cuando a X se le asignan más lecturas o UMI, el tamaño de la biblioteca de A aumenta, aumentando así el factor de tamaño de la biblioteca y produciendo valores de expresión normalizados más pequeños para todos los genes que no son DE. En ambos casos, el resultado final es que los genes distintos de DE en A estarán erróneamente regulados a la baja en comparación con B.

La eliminación del sesgo de composición es un problema bien estudiado para el análisis de grandes cantidades de datos de secuenciación de ARN. La normalización se puede realizar utilizando la función estimaSizeFactorsFromMatrix() en el paquete DESeq2 o la función calcNormFactors() en el paquete edgeR. Estos suponen que la mayoría de los genes no son DE entre células. Se supone que cualquier diferencia sistemática en el tamaño del recuento entre la mayoría de los genes que no son DE entre dos células representa un sesgo, que se utiliza para calcular un factor de tamaño apropiado para eliminarlo.

Sin embargo, aplicar estos métodos de normalización masiva a datos de una sola celda puede resultar problemático debido a la gran cantidad de recuentos bajos y cero. Para superar este problema, agregamos recuentos de muchas celdas para obtener estimaciones precisas del factor de tamaño. Luego, los factores de tamaño basados ​​en la biblioteca se "descompusieron" en factores de tamaño basados ​​en células para normalizar el perfil de expresión de cada célula. Esto se realiza utilizando la función ComputeSumFactors() de scran como se muestra a continuación.

Usamos un paso previo a la agrupación con quickCluster() donde las celdas de cada grupo se normalizan por separado y los factores de tamaño se reescalan para que sean comparables entre los grupos. Esto evita la suposición de que la mayoría de los genes de toda la población no son DE; solo se requiere una mayoría no DE entre pares de grupos, lo cual es una suposición más débil para poblaciones altamente heterogéneas. De forma predeterminada, quickCluster() utilizará un algoritmo de aproximación para PCA basado en los métodos del paquete irlba. La aproximación se basa en una inicialización aleatoria, por lo que debemos establecer una semilla aleatoria (a través de set.seed()) para lograr reproducibilidad.

Vemos que el factor de tamaño de deconvolución exhibe desviaciones específicas del tipo de celda del factor de tamaño de biblioteca en la Figura 7.2. Esto es consistente con la presencia de un sesgo compositivo introducido por una fuerte expresión diferencial entre tipos de células. El uso de factores de tamaño de deconvolución ajusta estos sesgos para mejorar la precisión de la normalización para aplicaciones posteriores.

La normalización precisa es más importante para los procesos que implican la estimación e interpretación de estadísticas por gen. Por ejemplo, el sesgo compositivo puede socavar el análisis DE al desplazar sistemáticamente el cambio logarítmico en una dirección u otra. Sin embargo, para los análisis basados ​​en celdas, como el análisis de conglomerados, a menudo proporciona menos beneficios que la simple normalización del tamaño de la biblioteca. La presencia de sesgo compositivo ya implica grandes diferencias en los perfiles de expresión, por lo que es poco probable que cambiar la estrategia de normalización afecte el resultado del proceso de agrupación.

La normalización del pico de entrada se basa en la suposición de que se añade la misma cantidad de pico de ARN a cada célula. Las diferencias sistemáticas en la cobertura de la transcripción de aumento se deben únicamente a sesgos específicos de las células, como la eficiencia de captura o la profundidad de secuenciación. Para eliminar estos sesgos, igualamos la cobertura de picos entre celdas escalando el "factor de tamaño de picos". A diferencia de los métodos anteriores, la normalización de picos no requiere suposiciones biológicas sistemáticas (es decir, no muchos genes DE). En cambio, se supone que la transcripción incorporada (i) se agrega a cada célula a un nivel constante y (ii) responde al sesgo de la misma manera relativa que el gen endógeno.

De hecho, si las diferencias en el contenido total de ARN de células individuales son motivo de preocupación y deben conservarse en análisis posteriores, se debe utilizar la normalización de picos. Para una célula determinada, un aumento en la cantidad total de ARN endógeno no aumenta su factor de tamaño de pico.

Esto garantiza que las diferencias de expresión en el contenido total de ARN entre poblaciones no se eliminen al realizar el escalado. Por el contrario, los otros métodos de normalización mencionados anteriormente simplemente interpretarán cualquier cambio en el contenido total de ARN como parte del sesgo y lo eliminarán.

Como ejemplo, utilizamos la normalización de picos en diferentes conjuntos de datos que involucran la activación de células T después de la estimulación con ligandos de receptores de células T de diferentes afinidades.

Aplicamos el método ComputeSpikeFactors() para estimar el pico. -en factores de tamaño para todas las celdas. Se define convirtiendo el recuento total de picos por celda en un factor de tamaño utilizando el mismo razonamiento que en LibrarySizeFactors(). Luego, el escalado eliminará cualquier diferencia en la cobertura de picos entre celdas.

Observamos una correlación positiva entre los factores de tamaño de pico y los factores de tamaño de deconvolución en cada condición de procesamiento (Figura 7.3), lo que indica que capturan tecnologías similares en profundidad de secuenciación y capturan el sesgo de eficiencia. Sin embargo, también observamos que el aumento de la estimulación del receptor de células T, en términos de aumento de la afinidad o del tiempo, dio como resultado una disminución en el factor de aumento en relación con el factor de tamaño de la biblioteca. Esto es consistente con un aumento en la actividad biosintética y el contenido total de ARN durante la estimulación, lo que reduce la cobertura relativa del pico en cada biblioteca (reduciendo así el factor de tamaño del pico) pero aumenta la cobertura de genes endógenos (por lo tanto, aumenta el factor de tamaño de la biblioteca) .

La diferencia entre los dos conjuntos de factores de tamaño tiene implicaciones prácticas para la interpretación posterior. Si se aplica un factor de tamaño de aumento a la matriz de recuento, los valores de expresión en las células no estimuladas se ampliarán, mientras que la expresión en las células estimuladas se reducirá. Sin embargo, si utiliza el factor de tamaño de deconvolución, sucede lo contrario. Cuando cambiamos entre estrategias de normalización, esto puede manifestarse como cambios en el tamaño y la dirección del DE entre condiciones, como se muestra a continuación para Malat1 (Figura 7.4).

Una vez calculado el factor de tamaño, el valor de la expresión normalizada se puede calcular para cada celda usando la función logNormCounts() en dispersión. Esto se hace dividiendo el recuento de cada gen/transcripción de adición por el factor de tamaño apropiado para esa célula. Esta función también transforma logarítmicamente los valores normalizados, creando un nuevo ensayo llamado "logcounts". Estos valores logarítmicos servirán como base para nuestro análisis posterior en los siguientes capítulos.

La transformación logarítmica es útil porque las diferencias en los valores logarítmicos representan cambios logarítmicos en la expresión genética. Esto es importante en los procesos posteriores basados ​​en la distancia euclidiana, que incluyen muchas formas de agrupación y reducción de dimensionalidad. Al operar con datos transformados logarítmicamente, nos aseguramos de que estos procedimientos midan las distancias entre las células en función de los cambios logarítmicos en la expresión genética. Por ejemplo, un gen con un nivel de expresión promedio de 50 en el tipo de célula A y 10 en el tipo de célula B, o un gen con un nivel de expresión promedio de 1100 en A y 1000 en B, puede mostrar una fuerte diferencia relativa mediante transformación logarítmica. , por lo que nos centraremos en el primero.

Al hacer una transformación logarítmica solemos añadir un pseudo recuento para evitar valores cero. Para genes de baja abundancia, pseudocuentas más grandes reducirán efectivamente el cambio logarítmico entre células a cero, lo que significa que los análisis posteriores de alta dimensión estarán impulsados ​​más por las diferencias de expresión de genes de alta abundancia. Por el contrario, pseudocuentas más pequeñas aumentarán la contribución relativa de genes de baja abundancia. Una práctica común es utilizar un pseudorecuento de 1 por la sencilla razón de que es práctico porque preserva la escasez en la matriz original (es decir, los ceros en la matriz original siguen siendo cero después de la transformación). Este método es eficaz en todos los casos excepto en la mayoría de las condiciones patológicas.

Por cierto, el aumento de pseudocuentas está motivado por centrar los factores de tamaño a la unidad. Esto garantiza que tanto el valor del pseudorecuento como el de la expresión normalizada estén en el mismo rango. Un pseudorecuento de 1 se puede interpretar como lecturas adicionales o UMI por gen.

En la práctica, centrar significa que el efecto de contracción del pseudoconteo disminuye a medida que aumenta la profundidad del conteo. Esto garantiza correctamente que las estimaciones de los cambios logarítmicos en la expresión (por ejemplo, basadas en diferencias en los valores logarítmicos entre grupos de células) sean cada vez más precisas a medida que aumenta la cobertura. Por el contrario, si se aplica un pseudorecuento constante a una métrica como partes por millón, la precisión de los cambios logarítmicos posteriores nunca mejorará, sin importar cuánta secuenciación adicional realicemos.

En casos raros, el escalado directo de los recuentos no es adecuado debido a los efectos descritos por A.Lun. En pocas palabras, esto se debe a la diferencia entre la media de los recuentos normalizados logarítmicamente y la media de los recuentos normalizados transformados logarítmicamente. La diferencia entre ellos depende de la media y la varianza de los recuentos brutos, por lo que existe una tendencia sistemática en la media de los recuentos logarítmicos en relación con el tamaño del recuento. Esto a menudo se manifiesta por el hecho de que las trayectorias están estrechamente relacionadas con el tamaño de la biblioteca incluso después de la normalización del tamaño de la biblioteca, como se muestra en la Figura 7.5. Los datos sintéticos de scRNA-seq generados por el método de fusión y división se muestran en la Figura 5.

Dado que el problema se debe a diferencias en el tamaño del recuento, la solución más sencilla es reducir la muestra de las celdas de alta cobertura para que coincidan con las celdas de baja cobertura. Esto utiliza el factor de tamaño para determinar la reducción de resolución de cada celda necesaria para alcanzar el primer percentil del factor de tamaño. (Solo unas pocas celdas con factores de tamaño más pequeños simplemente se amplían. No intentamos reducir la muestra al factor de tamaño más pequeño, ya que esto resultaría en una pérdida excesiva de información para una celda atípica con un factor de tamaño muy bajo). Podemos Observe que esto elimina los rastros en las dos primeras PC que están relacionados con el factor de tamaño de la biblioteca, mejorando así la resolución de las diferencias conocidas basadas en proporciones de mezcla (Figura 7.6). La transformación de registros sigue siendo necesaria, pero ya no produce cambios medios cuando los tamaños de recuento entre celdas son similares.

Si bien el submuestreo es una solución conveniente, estadísticamente es ineficaz debido a la necesidad de aumentar el ruido en las celdas de alta cobertura para evitar diferencias con las celdas de baja cobertura. También es más lento que el simple escalado. Por lo tanto, solo recomendamos utilizar este método después de que un análisis inicial de escala muestre trayectorias sospechosas que estén altamente correlacionadas con el factor de tamaño. En este caso, es sencillo volver a determinar si la trayectoria es un artefacto de la transformación logarítmica mediante reducción de resolución.