Proceso completo de análisis unicelular: normalización de datos
En este punto, debes prestar atención a la diferencia entre normalización y corrección por lotes. La normalización es independiente de la estructura del lote y solo tiene en cuenta las desviaciones técnicas, mientras que la corrección del lote solo ocurre entre lotes y se deben considerar tanto las desviaciones técnicas como las diferencias biológicas. Los sesgos técnicos tienden a afectar a los genes de manera similar o al menos en relación con sus características biofísicas (por ejemplo, longitud, contenido de GC), y las diferencias biológicas entre lotes pueden ser muy impredecibles. Por lo tanto, 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 de cada 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 el valor medio esperado para la célula. El factor de tamaño para cada celda representa una estimación del sesgo relativo en esa celda, por lo que dividir su recuento por el factor de tamaño debería eliminar el 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 el 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 de escala. Definimos el tamaño de la biblioteca como la suma de todos los recuentos de genes en cada célula, asumiendo que su valor esperado es proporcional a cualquier sesgo específico de la célula. Luego, sujeto a una constante de escala definida, el "factor de tamaño del grupo" de cada celda es proporcional al tamaño de su grupo, de modo que el factor de tamaño promedio en todas las celdas es igual a 1. Esta definición garantiza que los valores de la expresión normalizada estén en la misma escala que los recuentos originales, lo cual es muy útil para la interpretación, especialmente cuando se trabaja con datos transformados.
En los datos del cerebro de Zeisel, la mayor diferencia en el factor de tamaño de la biblioteca entre células fue de 10 veces. 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 "desequilibrio" en los (de)genes expresados diferencialmente entre cualquier par de células. Es decir, cualquier regulación positiva de un subconjunto de genes puede verse compensada por la misma regulación negativa de un subconjunto diferente de genes. Al evitar los efectos de la síntesis, esto garantiza que el tamaño de la biblioteca sea una estimación imparcial en relación con la especificidad celular. Sin embargo, la DE equilibrada generalmente no existe 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 normalizada 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, solo el rango de 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 propósito es identificar poblaciones de células y definir marcadores superiores para cada población de células.
Como se mencionó anteriormente, el sesgo compositivo ocurre cuando existe alguna expresión diferencial desequilibrada entre muestras. Tomemos el ejemplo de dos células, donde un único gen X está regulado positivamente en la célula A en relación con la célula B. Esta regulación positiva significa (I) que se dedican más recursos de secuenciación cuando el tamaño se determina experimentalmente (por ejemplo, debido a la cuantificación de la biblioteca); de otros genes no diferenciales disminuye, o (ii) cuando se asignan más lecturas o UMI a X, el tamaño de la biblioteca de A aumenta, aumentando así el factor de tamaño de la biblioteca, y se producen valores de expresión normalizados más pequeños para todos los genes no DE . En ambos casos, el resultado final es que los genes distintos de DE en A están erróneamente regulados a la baja en comparación con B.
Para el análisis de grandes cantidades de datos de secuenciación de ARN, la eliminación del sesgo de componentes es un problema bien estudiado.
La normalización se puede realizar utilizando la función estimaSizeFactorsFromMatrix() en el paquete DESeq2 o la función calcNormFactors() en el paquete edgeR. Estas hipótesis sugieren que la mayoría de los genes no están segregados entre células. Se supone que cualquier diferencia sistemática en el tamaño del recuento entre las dos células para la mayoría de los genes que no son DE 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 por lotes a datos de una sola celda puede resultar problemático debido a la gran cantidad de recuentos bajos y cero. Para superar este problema, sumamos los recuentos de muchas celdas para estimar con precisión el 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.
Utilizamos quickCluster() para un paso previo al agrupamiento en el que las celdas de cada grupo se normalizan por separado y el factor de tamaño se reescala a Comparable. 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 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.
Podemos ver que el factor de tamaño de deconvolución y el factor de tamaño de biblioteca en la Figura 7.2 muestran sesgos específicos del tipo de celda. Esto es consistente con la presencia de un sesgo compositivo introducido por una fuerte expresión diferencial entre tipos de células. El factor de tamaño de deconvolución se utiliza para ajustar estos sesgos y mejorar la precisión de la normalización para aplicaciones posteriores.
La normalización precisa es muy importante para el proceso involucrado en la estimación e interpretación de las estadísticas de cada gen. Por ejemplo, el sesgo compositivo puede socavar el análisis DE al desplazar sistemáticamente los cambios logarítmicos 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 de picos se basa en la suposición de que se añade la misma cantidad de ARN de picos a cada célula. Las diferencias sistemáticas en la cobertura de la transcripción de picos se deben simplemente a sesgos específicos de las células, como la eficiencia de captura o la profundidad de secuenciación. Para eliminar estos sesgos, equilibramos la cobertura de picos entre celdas escalando el "factor de tamaño de pico". En comparación con 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 de pico 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 necesita prestar atención a las diferencias en el contenido de ARN total en células individuales y debe preservarlo en análisis posteriores, debe utilizar la normalización estándar. Para una célula determinada, el factor de tamaño de la espiga no aumenta al aumentar el ARN endógeno total. 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 interpretan cualquier cambio en el contenido total de ARN como parte del sesgo y lo eliminan.
Por ejemplo, la normalización de picos se utiliza en diferentes conjuntos de datos que involucran la activación de células T después de la estimulación por ligandos de receptores de células T con diferentes afinidades.
Utilizamos el método CompuSpikeFactors() para estimar los factores de tamaño de pico para todas las celdas. El recuento máximo total por celda se convierte en un factor de tamaño utilizando el mismo razonamiento que en LibrarySizeFactors(). El escalado eliminará cualquier diferencia en la cobertura de picos entre celdas.
Observamos una correlación positiva entre el factor de tamaño aumentado y el factor de tamaño desconvolucionado bajo cada condición de procesamiento (Figura 7.3), lo que indica que capturan valores similares en términos de profundidad de secuenciación y eficiencia de captura. . Sin embargo, también observamos que la estimulación de los receptores de células T aumenta al aumentar la afinidad o el tiempo, lo que resulta en una disminución en el factor de pico en comparació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 de incorporación relativa en cada biblioteca (reduciendo así el factor de tamaño de incorporación) pero aumenta la cobertura de genes endógenos (aumentando así el factor de tamaño de 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 pico a la matriz de recuento, los valores de expresión en las células no estimuladas aumentarán proporcionalmente, mientras que los valores de expresión en las células estimuladas disminuirán proporcionalmente. 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 de la DE entre condiciones, como se muestra a continuación en Malat1 (Figura 7.4).
Una vez calculado el factor de tamaño, la función logNormCounts() en dispersión se puede utilizar para calcular el valor de expresión normalizado para cada celda. Esto se hace dividiendo el recuento por transcripción de gen/pico por el factor de tamaño apropiado para la célula. Esta función también transforma logarítmicamente los valores estandarizados, creando un nuevo análisis llamado "recuentos de registros". Estos valores logarítmicos se utilizarán como base para nuestro análisis posterior en secciones posteriores.
La transformación logarítmica es útil porque la diferencia en los valores logarítmicos indica el cambio logarítmico 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 procesar datos transformados logarítmicamente, nos aseguramos de que estos procedimientos midieran la distancia entre las células en función del cambio logarítmico en la expresión genética. Por ejemplo, un gen tiene un nivel de expresión promedio de 50 en células de tipo A y 100 en células de tipo B, o un gen tiene un nivel de expresión de 1100 en células de tipo A y 1000 en células de tipo B, y 1000 en células de tipo B relativamente fuerte. Se pueden mostrar diferencias en la transformación logarítmica, por lo que nos centraremos en la primera.
Cuando hacemos una transformación logarítmica, generalmente agregamos un conteo ficticio 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 en la expresión de genes de alta abundancia. Por el contrario, recuentos de errores más pequeños aumentan la contribución relativa de genes de baja abundancia. La 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 los patológicos.
Por cierto, el aumento de recuentos espurios se debe al factor de tamaño de la concentración. Esto garantiza que los valores del pseudorecuento y de la expresión normalizada estén en el mismo rango. Un recuento espurio de 1 puede interpretarse como un número adicional de lecturas por gen o UMI. De hecho, 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 múltiples cambios en la expresión (por ejemplo, basadas en diferencias en los valores logarítmicos entre grupos de celdas) se vuelvan cada vez más precisas a medida que aumenta la cobertura. Por el contrario, si se aplica un pseudorecuento constante a una métrica similar a partes por millón, entonces no importa cuánta secuenciación adicional hagamos, los cambios posteriores a una precisión multiplicada nunca mejorarán.
En casos raros, no es adecuado ajustar directamente el recuento 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 la media de los recuentos logarítmicos tiene una tendencia sistemática con respecto al tamaño del recuento. Esto generalmente muestra que incluso después de la normalización del tamaño de la biblioteca, las trayectorias están estrechamente relacionadas con el 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.
Debido a que el problema se debe a diferencias en los tamaños de conteo, la solución más sencilla es reducir la frecuencia de muestreo de las celdas de alta cobertura para que coincida 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 escalan. No intentamos reducir el muestreo al factor de tamaño más pequeño, ya que esto resultaría en una pérdida excesiva de información para celdas atípicas con factores de tamaño muy bajos. Podemos ver que esto elimina Las dos primeras PC están correlacionadas con el factor de tamaño de la biblioteca, mejorando así la resolución de diferencias conocidas basadas en proporciones de mezcla (Figura 7.6). La transformación logarítmica todavía es necesaria, pero cuando los recuentos entre celdas son similares, ya no resulta. un cambio en la media.
Si bien la reducción de resolución es una solución conveniente, estadísticamente es ineficaz porque requiere aumentar el ruido de las celdas de alta cobertura para evitar diferencias con respecto a 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 el 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.