Perfiles transcripcionales digitales masivamente paralelos de células individuales

Imágenes de alta velocidad de microesferas de gel y células en GEMs

Se utilizó un microscopio (Nikon Ti-E, objetivo × 10) y una cámara de vídeo de alta velocidad (Photron SA5, velocidad de fotogramas=4.000 s-1) para obtener imágenes de cada GEM a medida que se generaban en el chip microfluídico. Se utilizó un software de análisis personalizado para contar el número de GEMs generados y el número de cuentas presentes en cada GEM, basándose en la detección de bordes y el contraste entre los bordes de las cuentas y los GEMs y el líquido adyacente. Los resultados del análisis se resumen en la Fig. 1c. Para estimar la distribución de las células en los GEM, se utilizó el recuento manual para ∼28k fotogramas de un vídeo en un subconjunto de GEM. Los resultados indican una adhesión aproximada a una distribución de Poisson. Sin embargo, el porcentaje de encapsulaciones celulares múltiples fue un 16% superior al valor esperado, posiblemente debido a un error de submuestreo o a las interacciones célula-célula (se observaron algunos grupos de dos células durante el recuento manual) (Supplementary Fig. 1b).

Líneas celulares y muestras de pacientes trasplantados

Las células Kurkat (ATCC TIB-152), 293T (ATCC CRL-11268) y 3T3 (ATCC CRL-1658) se adquirieron de ATCC y se cultivaron según las directrices de ATCC. Las PBMC frescas, las PBMC congeladas y las BMMC se adquirieron en ALLCELLS. Las PBMC congeladas del donante A se hicieron a partir de PBMC frescas del donante A mezclando 1e6 células en medio de congelación (15% de dimetilsulfóxido (DMSO) en medio de Dulbecco modificado de Iscove que contenía 20% de FBS) suavemente, y se enfriaron en CoolCell FTS30 (BioCision) en -80 °C durante al menos 4 h antes de transferirlas a nitrógeno líquido para su almacenamiento durante 3 semanas.

La Junta de Revisión Institucional del Centro de Investigación del Cáncer Fred Hutchinson aprobó el estudio de las muestras de trasplante. Los procedimientos seguidos se ajustaron a la Declaración de Helsinki de 1975 y a la Regla Común. Las muestras se obtuvieron después de que los pacientes dieran su consentimiento informado por escrito sobre los análisis moleculares. Se identificaron pacientes con LMA sometidos a trasplante alogénico de células madre hematopoyéticas en el Fred Hutchinson Cancer Research Center. El diagnóstico de LMA se estableció de acuerdo con los criterios revisados de la Organización Mundial de la Salud33.

Los aspirados de médula ósea se obtuvieron para las pruebas clínicas estándar 20-30 días antes del trasplante y en serie después del trasplante de acuerdo con el protocolo de tratamiento. Las alícuotas de los aspirados de médula ósea se procesaron en las 2 horas siguientes a la extracción. Las BMMC se aislaron mediante centrifugación a través de un gradiente de Ficoll (Histopaque-1077; Sigma Life Science, St Louis, MO, USA). Las BMMC se recogieron de la interfaz suero-Ficoll con una pipeta Pasteur desechable y se transfirieron al tubo cónico de 50 ml con suero del paciente al 2% en 1 × PBS. Las BMMC se contaron con un hemocitómetro y la viabilidad se evaluó con azul tripán. Las BMMC se resuspendieron en medio de congelación 90% FBS, 10% DMSO y se congelaron utilizando un Thermo Scientific Nalgene Mr Frosty (Thermo Scientific) en un congelador a -80 °C durante 24 h antes de ser transferidas a nitrógeno líquido para su almacenamiento a largo plazo.

Estimación del contenido de ARN por célula

La cantidad de ARN por tipo de célula se determinó cuantificando (Qubit; Invitrogen) el ARN extraído (Maxwell RSC simplyRNA Cells Kit) de varios números diferentes conocidos de células.

Preparación de las células

Se cosecharon células frescas, se lavaron con 1 × PBS y se resuspendieron a 1 × 106 células por ml en 1 × PBS y 0,04% de albúmina de suero bovino. Las PBMC frescas se congelaron a 10 × resuspendiendo las PBMC en DMEM+40% FBS+10% DMSO, congelando a -by °C en un CoolCell® FTS30 (BioCision) y colocando después en nitrógeno líquido para su almacenamiento.

Los viales de células congeladas de los estudios ALLCELLS y de trasplante se descongelaron rápidamente en un baño de agua a 37 °C durante ∼2 min. Los viales se retiraron cuando quedó un pequeño cristal de hielo. Las PBMC descongeladas se lavaron dos veces en el medio y luego se resuspendieron en 1 × PBS y 0,04% de albúmina de suero bovino a temperatura ambiente. Las células se centrifugaron a 300 r.c.f. durante 5 minutos cada vez. Las BMMC descongeladas se lavaron y resuspendieron en 1 × PBS y 20% de FBS. La concentración final de las células descongeladas fue de 1 × 106 células por ml.

Construcción de bibliotecas de secuenciación utilizando la plataforma GemCode

Las suspensiones celulares se cargaron en un GemCode Single-Cell Instrument (10x Genomics, Pleasanton, CA, USA) para generar GEMs unicelulares. Las bibliotecas de ARN-Seq unicelulares se prepararon utilizando el kit GemCode Single-Cell 3′ Gel Bead and Library Kit (ahora vendido como P/N 120230, 120231, 120232, 10x Genomics). La GEM-RT se realizó en un termociclador C1000 Touch con 96-Deep Well Reaction Module (Bio-Rad; P/N 1851197): 55 °C durante 2 h, 85 °C durante 5 min; se mantuvo a 4 °C. Después de la RT, se rompieron los GEM y el ADNc monocatenario se limpió con DynaBeads MyOne Silane Beads (Thermo Fisher Scientific; P/N 37002D) y el kit de reactivos SPRIselect (0,6 × SPRI; Beckman Coulter; P/N B23318). El ADNc se amplificó utilizando el termociclador C1000 Touch con el módulo de reacción de 96 pozos profundos: 98 °C durante 3 minutos; ciclos de 14 × : 98 °C durante 15 s, 67 °C durante 20 s y 72 °C durante 1 min; 72 °C durante 1 min; se mantuvo a 4 °C. El producto de ADNc amplificado se limpió con el kit de reactivos SPRIselect (0,6 × SPRI). Posteriormente, el ADNc se cortó a ∼200 pb utilizando un sistema Covaris M220 (Covaris; P/N 500295). Las bibliotecas de secuenciación indexadas se construyeron utilizando los reactivos del kit de bibliotecas GemCode Single-Cell 3′, siguiendo estos pasos: (1) reparación de extremos y cola A; (2) ligadura de adaptadores; (3) limpieza posterior a la ligadura con SPRIselect; (4) PCR de índice de muestra y limpieza. Las bibliotecas de secuenciación con código de barras se cuantificaron mediante PCR cuantitativa (Kit de cuantificación de bibliotecas de KAPA Biosystems para plataformas Illumina P/N KK4824). Las bibliotecas de secuenciación se cargaron a 2,1 pM en un Illumina NextSeq500 con kits de 2 × 75 paired-end utilizando la siguiente longitud de lectura: 98 bp Read1, 14 bp I7 Index, 8 bp I5 Index y 10 bp Read2. Algunas bibliotecas anteriores se hicieron con UMI de 5 nt, y en su lugar se obtuvo Read2 de 5 pb. Estas bibliotecas se han documentado en la Tabla Suplementaria 1.

Ensayo de ERCC

Los ARN sintéticos de spike-in (Thermo Fisher Scientific; P/N 4456740) se diluyeron (1:10 o 1:50) y se cargaron en un GemCode Single-Cell Instrument, sustituyendo a las células normalmente utilizadas para generar GEMs. Se probaron las mezclas Spike-in Mix1 y Mix2. Se utilizó un protocolo ligeramente modificado, ya que sólo se recogió una pequeña fracción de GEM para la RT y la amplificación del ADNc. Tras la finalización de la GEM-RT, se retiraron 1,25 μl de la emulsión y se añadieron a una mezcla bifásica de agente de recuperación (125 μl) (P/N 220016) y aditivo 1 de 25 mM (30 μl) (P/N 220074, 10x Genomics). A continuación, se eliminó el agente de recuperación y se limpió la solución acuosa restante con el kit de reactivos SPRISelect (0,8 × SPRI). El ADNc se amplificó utilizando el termociclador C1000 Touch con módulo de reacción de 96 pozos profundos: 98 °C durante 3 minutos; ciclos de 14 × : 98 °C durante 15 s, 67 °C durante 20 s y 72 °C durante 1 min; 72 °C durante 1 min; se mantuvo a 4 °C. El producto de ADNc amplificado se limpió con el kit de reactivos SPRIselect (0,8 × ) El ADNc se cortó posteriormente a ∼200 pb utilizando un sistema Covaris M220 para construir bibliotecas indexadas por muestra con adaptadores 10x Genomics. Los recuentos de moléculas de ERCC esperados se calcularon en función de la cantidad de moléculas de ERCC utilizadas y de los factores de dilución de la muestra. Los recuentos se compararon con los recuentos de moléculas detectadas (recuentos UMI) para calcular la eficiencia de conversión.

ensayo ddPCR

Las células Jurkat se utilizaron en ensayos ddPCR para estimar la eficiencia de conversión de la siguiente manera: (1) se determinó la cantidad de ARN por célula Jurkat cuantificando (Qubit, Invitrogen) el ARN extraído (Kits de purificación de ARN Maxwell) de varios números diferentes conocidos de células Jurkat. (2) Se realizó una RT-ddPCR masiva (Bio-Rad One-Step RT-ddPCR Advanced Kit for Probes 1864021) sobre el ARN extraído para determinar el número de copias por célula de ocho genes seleccionados. (3) Se procesaron aproximadamente 5.000 células Jurkat utilizando la plataforma GemCode Single-Cell 3′, y se recogió ADNc monocatenario después de la RT en GEMs siguiendo los protocolos indicados en la sección «Construcción de bibliotecas de secuenciación utilizando la plataforma GemCode». Se determinaron las copias de ADNc de los ocho genes utilizando ddPCR (Bio-Rad ddPCR Supermix for Probes (no dUTP) P/N 1863024). El recuento real de células Jurkat se determinó secuenciando un subconjunto de las reacciones GEM-RT en un MiSeq. La eficiencia de conversión es la relación entre las copias de ADNc por célula (paso 3) y las copias de ARN por célula de la RT-ddPCR masiva (paso 2), asumiendo una eficiencia del 50% en la RT-ddPCR34.

Las secuencias de la sonda para el ensayo de ddPCR son las siguientes: SERAC1_f, 5′-CACGAGGCCAGC-3′ y SERAC1_r, 5′-TCTGCAACAGATGACGCAATAAG-3′; SERAC1_p: /56-FAM/CGCCTGCCG/ZEN/GCAGAATGTC/3IABkFQ/. AP1S3_f, 5′-GAAGCAGCCATGTCTAAGC-3′ y AP1S3_r, 5′-CCTTGTCGACTGAAGCAATATG-3′; AP1S3_p: /56-FAM/CGGCCCAGC/ZEN/CACGATGATACAT/3IABkFQ/OR. AOV1_f, 5′-CCGGAAGTGGTCTCGTOR-3′ y AOV1_r, 5′-TTCTTCATAGCCTTCCCGATACCOR-3′; AOV1_p: /56-FAM/TCGTGATGG/ZEN/CGGATGAGGTTTCA/3IABkFQ/. DOLPP1_f, 5′-ATGGCAGCGACGGA-3′ y DOLPP1_r, 5′-GGCTCAGGTAGCAAGGA-3′; DOLPP1_p: /56-FAM/CCACGTCGA/ZEN/ATATCCTGCAGGTGATCT/3IABkFQ/. KPNA6_f, 5′-TGAAAGCTGCCGCTGAAG-3′ y KPNA6_r, 5′-CCCTGGCTCGCCAT-3′; KPNA6_p: /56-FAM/CGGACCCGC/ZEN/GATGGAGACC/3IABkFQ/. ITSN2_f, 5′-GTGACAGGCTACGCAACAG-3′ e ITSN2_r, 5′-TCCTGAGTTTTCCTTGCTAGCT-3′; ITSN2_p: /56-FAM/AGGGCGCCA/ZEN/GATGGCTGA/3IABkFQ/. LCMT1_f, 5′-GTCGACCGCTTCCA-3′ y LCMT1_r, 5′-GGTCATGCCAGTAGCCAATG-3′; LCMT1_p: /56-FAM/ATGCTTCCC/ZEN/TGTGCAAGGTTTGC/3IABkFQ/. AP2M1_f, 5′-GCAGCGGCAGACG-3′ y AP2M1_r, 5′-ATGGCGCAGATCAGTCT-3′; AP2M1_p: /56-FAM/CATCGCTCT/ZEN/GAGAACAGACCTGGTG/3IABkFQ/.

Cálculo de la eficiencia de captura de células

La eficiencia se calcula tomando la relación entre el número de células detectadas por secuenciación frente al número de células cargadas en el chip. Este último se determina a partir de (volumen añadido × concentración de entrada de células). La concentración de células de entrada se determinó utilizando un contador celular automatizado Countess II (Thermo Fisher Scientific). Cabe señalar que existe un error del 15-20% en el recuento de células, lo que podría explicar al menos parte de la variabilidad en las eficiencias calculadas.

Asunto de quimerismo

Sistema PowerPlex 16 (Promega) se utilizó junto con un analizador genético 3130xl de Applied Biosystems (Life Technologies). Se utilizaron BMMC de donantes como línea de base de referencia.

Alineación, asignación de códigos de barras y recuento de UMI

Se utilizó el Cell Ranger Single-Cell Software Suite para realizar el demultiplexado de las muestras, el procesamiento de los códigos de barras y el recuento de genes 3′ de una sola célula (http://software.10xgenomics.com/single-cell/overview/welcome). En primer lugar, se llevó a cabo la demultiplexación de la muestra basándose en la lectura de índice de muestra de 8 pb para generar FASTQs para las lecturas de fin de pares de Read1 y Read2, así como el código de barras GemCode de 14 pb. Se extrajeron etiquetas UMI de diez pares de bases de Read2 (14 bibliotecas se hicieron con etiquetas UMI de 5 pb, como se indica en la Tabla Suplementaria 1, debido a una iteración anterior de los métodos. Para estas muestras, las etiquetas UMI de 5 pb se extrajeron de Read2). A continuación, Read1, que contiene el inserto de ADNc, se alineó con un genoma de referencia apropiado utilizando STAR35. Para las células de ratón, se utilizó mm10 y para las células humanas, hg19. Para las muestras con mezclas de células de ratón y humanas, se utilizó la unión de hg19 y mm10. Para las muestras ERCC, se utilizó la referencia ERCC (https://tools.thermofisher.com/content/sfs/manuals/cms_095047.txt).

A continuación, se filtraron los códigos de barras GemCode y los UMI. Se consideraron todas las listas conocidas de códigos de barras que están a una distancia de 1-Hamming de un código de barras observado. A continuación, se calcula la probabilidad posterior de que el código de barras observado se haya producido por un error de secuenciación, dadas las calidades de base del código de barras observado y la probabilidad previa de observar el código de barras candidato (tomada de la distribución global de recuento de códigos de barras). Si la probabilidad posterior de cualquier código de barras candidato es de al menos 0,975, el código de barras se corrige al código de barras candidato con la mayor probabilidad posterior. Si todas las secuencias candidatas son igualmente probables, entonces se elige la que aparece primero por orden léxico.

Los UMIs con puntuación de calidad de secuenciación >10 se consideraron válidos si no eran homopolímeros. Qual=10 implica un 90% de precisión en la llamada de bases. Una UMI que está a una distancia de 1-Hamming de otra UMI (con más lecturas) para el mismo código de barras de la célula y el gen se corrige a la UMI con más lecturas. Este enfoque es casi idéntico al de Jaitin et al.4, y es similar al de Klein et al.8 (aunque Klein et al.8 también utilizó UMIs para resolver lecturas multimapping, lo que no se implementó aquí).

Por último, se marcaron los duplicados de PCR si dos conjuntos de pares de lecturas compartían una secuencia de código de barras, una etiqueta UMI y un ID de gen (se utilizaron los GTFs de Ensembl GRCh37.82, ftp://ftp.ensembl.org/pub/grch37/release-84/gtf/homo_sapiens/Homo_sapiens.GRCh37.82.gtf.gz y GRCm38.84, ftp://ftp.ensembl.org/pub/release-84/gtf/mus_musculus/Mus_musculus.GRCm38.84.gtf.gz). Sólo se utilizaron duplicados mapeados con confianza (MAPQ=255), no-PCR con códigos de barras y UMIs válidos para generar la matriz gen-código de barras.

Los códigos de barras de las células se determinaron basándose en la distribución de los recuentos de UMI. Todos los códigos de barras superiores dentro del mismo orden de magnitud (>10% del código de barras superior n, donde n es el 1% del recuento de células recuperadas previsto) se consideraron códigos de barras de células. El número de lecturas que proporcionan información significativa se calcula como el producto de cuatro métricas: (1) códigos de barras válidos; (2) UMI válidos; (3) asociados con un código de barras de células; y (4) mapeados con confianza a los exones.

En los experimentos de mezcla de ratones y humanos, la tasa de multipletes se definió como el doble de la tasa de códigos de barras de células con recuentos de UMI significativos tanto de ratones como de humanos, donde el 1% superior de los recuentos de UMI se consideró significativo. El grado de interferencia de los códigos de barras se evaluó por la fracción de lecturas de ratón en los códigos de barras humanos, o viceversa.

Las muestras procesadas de múltiples canales pueden combinarse concatenando matrices de genes-células-códigos de barras. Esta funcionalidad se proporciona en el Cell Ranger R Kit (http://support.10xgenomics.com/single-cell/software/pipelines/latest/rkit). Los datos de secuenciación de múltiples ejecuciones de secuenciación de una biblioteca pueden combinarse contando las lecturas no duplicadas. Esta funcionalidad se proporciona en el pipeline de Cell Ranger. Además, los datos de secuenciación pueden ser submuestreados para obtener un número determinado de recuentos UMI por célula. Esta funcionalidad también se proporciona en el Cell Ranger R Kit, y es útil cuando se combinan datos de múltiples muestras para su comparación.

Análisis PCA de la mezcla de células Jurkat y 293T

Se concatenó la matriz de códigos de barras de células de genes de cada una de las cuatro muestras. Sólo se utilizaron los genes con al menos un recuento de UMI detectado en al menos una célula. La normalización de UMI se realizó dividiendo primero los recuentos de UMI por los recuentos totales de UMI en cada célula, seguido de una multiplicación por la mediana de los recuentos totales de UMI en todas las células. A continuación, se tomó el logaritmo natural de los recuentos de UMI. Por último, cada gen fue normalizado de tal manera que la señal media para cada gen es 0, y la desviación estándar es 1. PCA se ejecutó en la matriz normalizada de genes y códigos de barras. Los recuentos UMI normalizados de cada gen se utilizan para mostrar la expresión de un marcador en un gráfico tSNE.

Análisis de SNV de los datos scRNA-seq de Jurkat y 293T

Los SNVs se llamaron ejecutando Freebayes 1.0.2 (ref. 36) en el BAM del genoma producido por Cell Ranger. Se seleccionaron los SNV de alta calidad (SNV calling Qual>=100 con al menos 10 recuentos de UMI de al menos dos células; se ignoran los indels) que sólo se observaron en células Jurkat o 293T (pero no en ambas). Las células se etiquetaron como Jurkat o 293T basándose en los recuentos de SNV específicos de Jurkat y 293T, donde la fracción de recuentos de la otra especie es <0,2. Las células con una fracción de SNV de cualquiera de las dos especies entre 0,2 y 0,8 se consideran multipletes. La tasa de multipletes inferida es 2* la tasa de multipletes observada (para tener en cuenta los multipletes Jurkat:Jurkat y 293T:293T).

Análisis PCA y tSNE de PBMCs

Se utilizan genes con al menos un recuento UMI detectado en al menos una célula. Se identificaron los 1.000 genes más variables en función de su media y dispersión (varianza/media), lo que es similar al enfoque utilizado por Macoscko et al.7 Los genes se colocaron en 20 intervalos en función de su expresión media. La dispersión normalizada se calcula como la diferencia absoluta entre la dispersión y la dispersión mediana de la media de expresión, normalizada por la desviación absoluta mediana dentro de cada bin.

PCA se ejecutó en la matriz normalizada gen-código de barras de los 1.000 genes más variables para reducir el número de dimensiones de características (genes). La normalización de UMI se realizó dividiendo primero los recuentos de UMI por los recuentos totales de UMI en cada celda, seguido de una multiplicación con la mediana de los recuentos totales de UMI en todas las celdas. A continuación, se tomó el logaritmo natural de los recuentos de UMI. Por último, cada gen se normalizó de forma que la señal media de cada gen fuera 0 y la desviación estándar fuera 1. Se ejecutó el PCA en la matriz normalizada de genes y códigos de barras. Después de ejecutar el PCA, se realizó la aproximación de Barnes-hut37 a t-SNE16 en los primeros 50 PCs para visualizar las células en un espacio bidimensional. Se utilizaron 50 PCs porque (1) utilizar todos los PCs llevaría mucho tiempo con el análisis tSNE; (2) explicaban ∼25% de la varianza total. Se ejecutó la agrupación K-means15 para agrupar las células para el análisis de agrupación. Se seleccionó k=10 basándose en el gráfico de la suma del error cuadrado (Fig. 5d suplementaria).

Identificación de los genes específicos del clúster y clasificación basada en marcadores

Para identificar los genes que están enriquecidos en un clúster específico, se calculó la expresión media de cada gen en todas las células del clúster. A continuación, se comparó cada gen del clúster con la expresión media del mismo gen de las células de todos los demás clústeres. Los genes se clasificaron en función de su diferencia de expresión y se seleccionaron los 10 genes más enriquecidos de cada clúster. Para la agrupación jerárquica, se calculó la correlación por pares entre cada clúster y se utilizó la expresión centrada de cada gen para su visualización mediante un mapa de calor.

La clasificación de las PBMC se dedujo de la anotación de los genes específicos de los clústeres. En el caso del clúster 10, se detectó la expresión de marcadores de múltiples tipos de células (por ejemplo, B, dendríticas y T). Dado que el tamaño relativo de los clústeres B, dendríticos y T es del 5,7%, 6,6% y 81%, respectivamente, cabría esperar que el clúster 10 (que es sólo el 0,5%) contuviera multipletes formados principalmente por B:dendríticos (0.36%) y B:dendrítico:T (0,3%).

Selección de subpoblaciones purificadas de PBMCs

Cada población de PBMCs purificadas fue muestreada a ∼16k lecturas por célula. Se realizaron PCA, tSNE y k-means clustering para cada matriz downsampled, siguiendo los mismos pasos descritos en el análisis PCA y t-SNE de PBMCs. Sólo se detectó un clúster en la mayoría de las muestras, en consonancia con los análisis FACS (Fig. 6 suplementaria). En el caso de las muestras con más de un clúster, sólo se seleccionaron para el análisis posterior los clústeres que mostraban la expresión del gen marcador esperado. En el caso de los monocitos CD14+, se observaron dos agrupaciones que se identificaron como monocitos CD14+ y células dendríticas en función de la expresión de los genes marcadores FTL y CLEC9A, respectivamente.

Análisis de clasificación de células utilizando PBMCs purificadas

Cada población de PBMCs purificadas se sometió a un muestreo descendente de ∼16k lecturas mapeadas con confianza por célula. A continuación, se calculó un perfil de expresión génica promedio (media) en todas las células. A continuación, se comparó la expresión génica de cada célula de la población compleja con los perfiles de expresión génica de las poblaciones purificadas de PBMC mediante la correlación de Spearman. Se asignó a la célula el ID de la población purificada si tenía la mayor correlación con esa población. Obsérvese que la diferencia entre la correlación más alta y la segunda más alta era pequeña para algunas células (por ejemplo, la diferencia entre las células T citotóxicas y las células NK), lo que sugiere que la asignación celular no era tan segura para estas células. Algunas de las poblaciones de PBMC purificadas se solapaban entre sí. Por ejemplo, las células T auxiliares CD4+ incluyen todas las células CD4+. Esto significa que las células de esta muestra se solaparán con las células de las muestras que contienen células CD4+, incluyendo CD4+/CD25+ T reg, CD4+/CD45RO+ T memory, CD4+/CD45RA+/CD25- naive T. Por lo tanto, cuando a una célula se le asignó el ID de célula T-helper CD4+ basándose en la puntuación de correlación, se comprobó la siguiente correlación más alta para ver si era una de las muestras CD4+. Si lo era, el ID de la célula se actualizaba al tipo de célula con la siguiente correlación más alta. El mismo procedimiento se llevó a cabo para los T citotóxicos CD8+ y los T citotóxicos ingenuos CD8+/CD45RA+ (que es un subconjunto de los T citotóxicos CD8+).

El código R utilizado para analizar las PBMC de 68k y las PBMC purificadas puede encontrarse aquí: https://github.com/10XGenomics/single-cell-3prime-paper.

Clasificación y agrupación de células con Seurat

La matriz gen-celda-código de barras de 68k PBMCs fue transformada log como entrada a Seurat. Los 469 genes más variables seleccionados por Seurat se utilizaron para calcular las PC. Los primeros 22 PCs fueron significativos (P<0.01) basados en el análisis jackstraw incorporado, y se utilizaron para la visualización tSNE. La clasificación de las células se tomó del análisis de clasificación de las células utilizando PBMCs purificadas.

Comparación entre PBMCs frescas y congeladas

Los datos de secuenciación de 68k PBMCs frescas y 3k PBMCs congeladas fueron muestreados de tal manera que cada muestra tiene ∼14k lecturas mapeadas con confianza por célula. Sólo se incluyeron en la comparación los genes detectados en al menos una célula, que utiliza la media de cada gen en todas las células.

Para la comparación de la clasificación celular entre PBMC purificadas y congeladas, agrupamos todas las células etiquetadas como células T o natural killer. Esto se debe a que las subpoblaciones dentro de las células T y entre las células T y natural killer son a veces difíciles de agrupar por separado. No queríamos que la comparación entre las células frescas y las congeladas se viera afectada por los métodos de agrupación utilizados.

Asignación de genotipos basada en SNV

Los SNVs fueron llamados ejecutando Freebayes 1.0.2 (ref. 36) en el BAM del genoma producido por Cell Ranger. Sólo se incluyeron los SNV con apoyo de al menos dos códigos de barras de células, con una puntuación mínima de SNV Qual >=30, base mínima de SNV Qual>=1. Se calcularon los recuentos de alelos de referencia (R) y alternativos (A) en cada SNV, produciendo una matriz de recuentos de UMI de referencia de células y de UMI de alelos alternativos de células. Estas matrices se modelaron como una mezcla de dos genomas en la que la probabilidad de cualquiera de los tres genotipos (R/R, R/A o A/A) en un sitio se tomó como una distribución binomial con una tasa de error fija del 0,1%. Para cada muestra, se infirieron dos modelos en paralelo, uno en el que sólo está presente un genoma (K=1) y otro en el que están presentes dos genomas (K=2). La inferencia de los parámetros del modelo (asignaciones de célula a genoma y los K conjuntos de genotipos) se realizó utilizando un muestreador de Gibbs para aproximar sus distribuciones posteriores. Para mejorar el problema del cambio de etiquetas en la inferencia de Monte Carlo de los modelos de mezcla, el reetiquetado de las asignaciones de célula a genoma muestreadas se realizó según Stephens et al.38

En los experimentos de mezcla de células in silico, cuando el modelo K=2 no pudo separar adecuadamente los dos genomas, informó de una distribución de probabilidades posteriores cercana a 0,5 para las llamadas de célula a genoma, indicando una falta de confianza en esas llamadas. Aplicamos el requisito de que el 90% de las células tengan una probabilidad posterior >75% para seleccionar el modelo K=2 sobre el modelo K=1. La selección de K=1 indica que la fracción de mezcla está por debajo del nivel de detección del método, que en los experimentos de mezcla in silico se determinó que era el 4% de 6.000 células.

Comparación de genotipos con la muestra pura

Para determinar la asignación de genotipos a los individuos, sólo se consideraron los SNV compartidos entre el grupo de genotipos y la muestra pura. A continuación, se comparó el genotipo medio de todas las células con el de la muestra pura. Para obtener una línea de base para el % de solapamiento de genotipos entre diferentes individuos, realizamos una comparación por pares de genotipos llamados de los mismos individuos (11 comparaciones por pares) o de diferentes individuos (15 comparaciones por pares). El porcentaje de solapamiento de genotipos entre los mismos individuos tiene un promedio de ∼98±0,3%, mientras que el porcentaje de solapamiento de genotipos entre los diferentes individuos tiene un promedio de ∼73±2%.

Análisis PCA y tSNE de BMMCs

Se utilizaron datos de seis muestras: dos controles sanos, AML027 antes y después del trasplante, y AML035 antes y después del trasplante. Cada muestra fue muestreada hasta ∼10k lecturas mapeadas con confianza por célula. A continuación, se concatenó la matriz de códigos de barras de genes de cada muestra. Se realizó un PCA, un tSNE y una agrupación de k-means en la matriz combinada, siguiendo los mismos pasos descritos en el análisis PCA y tSNE de PBMC. Para la agrupación de k-means, se utilizó K=10 basándose en la curvatura de la suma del error cuadrático del diagrama de dispersión.

Los genes específicos de la agrupación se identificaron siguiendo los pasos descritos en «Identificación de genes específicos de la agrupación y clasificación basada en marcadores». La clasificación se asignó basándose en los genes específicos de los clusters y en la expresión de algunos marcadores bien conocidos de los tipos de células inmunitarias. Blasto y Ery 1 Inmaduro’ se refiere al cluster 4, que expresa CD34, un marcador de progenitores hematopoyéticos39, y Gata2, un marcador de eritroides tempranos40. El «Ery 2 inmaduro» se refiere a los grupos 5 y 8, que muestran la expresión de Gata1, un factor de transcripción esencial para la eritropoyesis41, pero no de CD71, que suele encontrarse en las células eritroides más comprometidas39. El término «Ery 3 inmaduro» se refiere al grupo 1, que muestra la expresión de CD71. El término «Ery maduro» se refiere al grupo 2. HBA1, un marcador de células eritroides maduras, se detecta preferentemente en el grupo 2. El cluster 3 fue asignado como ‘Granulocitos Inmaduros’ debido a la expresión de marcadores de granulocitos tempranos como AZU1 e IL8 (ref. 42), y la falta de expresión de CD16. El cluster 7 fue asignado como ‘Monocitos’ debido a la expresión de CD14 y FCN1, por ejemplo. ‘B’ se refiere a los clusters 6 y 9 debido a marcadores como CD19 y CD79A. ‘T’ se refiere al clúster 10, debido a marcadores como CD3D y CD8A.

Disponibilidad de datos

Todos los datos relevantes están disponibles en los autores. Los datos de RNA-seq de una sola célula se han depositado en el Short Read Archive con el número de acceso SRP073767. Los datos también están disponibles en http://support.10xgenomics.com/single-cell/datasets. El código de análisis para el análisis de 68k PBMC está disponible en https://github.com/10XGenomics/single-cell-3prime-paper.