Massively parallel digital transcriptional profiling of single cells

High-speed imaging of gel beads and cells in GEMs

A microscope (Nikon Ti-E, × 10 objective) and a high-speed video camera (Photron SA5, frame rate=4,000 s-1) were used to image every GEM as they were generated in the microfluidic chip. Un software di analisi personalizzato è stato utilizzato per contare il numero di GEMs generato e il numero di perline presenti in ogni GEM, sulla base di rilevamento del bordo e il contrasto tra i bordi di perline e bordi GEM e il liquido adiacente. I risultati dell’analisi sono riassunti in Fig. 1c. Per stimare la distribuzione delle cellule in GEMs, conteggio manuale è stato utilizzato per ∼28k fotogrammi di un video su un sottoinsieme di GEMs. I risultati indicano un’approssimativa aderenza a una distribuzione di Poisson. Tuttavia, la percentuale di incapsulamenti cellulari multipli è stata del 16% superiore al valore atteso, probabilmente a causa di un errore di sottocampionamento o di interazioni cellula-cellula (alcuni grumi di due cellule sono stati osservati durante il conteggio manuale) (Fig. 1b).

Line cellulari e campioni di pazienti trapiantati

Le cellule Jurkat (ATCC TIB-152), 293T (ATCC CRL-11268) e 3T3 (ATCC CRL-1658) sono state acquisite da ATCC e coltivate secondo le linee guida ATCC. Le PBMC fresche, le PBMC congelate e le BMMC sono state acquistate dalla ALLCELLS. Le PBMC congelate del donatore A sono state ottenute da PBMC fresche del donatore A mescolando 1e6 cellule in mezzo di congelamento (15% dimetilsolfossido (DMSO) in mezzi Dulbecco modificati di Iscove contenenti 20% FBS) delicatamente, e raffreddate in CoolCell FTS30 (BioCision) a -80 °C per almeno 4 h prima del trasferimento in azoto liquido per la conservazione per 3 settimane.

L’Institutional Review Board del Fred Hutchinson Cancer Research Center ha approvato lo studio sui campioni da trapianto. Le procedure seguite erano conformi alla Dichiarazione di Helsinki del 1975 e alla Common Rule. I campioni sono stati ottenuti dopo che i pazienti hanno fornito il consenso informato scritto sulle analisi molecolari. Abbiamo identificato i pazienti con AML sottoposti a trapianto allogenico di cellule staminali ematopoietiche presso il Fred Hutchinson Cancer Research Center. La diagnosi di AML è stata stabilita secondo i criteri rivisti dell’Organizzazione Mondiale della Sanità33.

Aspirati di midollo osseo sono stati ottenuti per i test clinici standard 20-30 giorni prima del trapianto e in serie post-trapianto secondo il protocollo di trattamento. Le aliquote di aspirato di midollo osseo sono state elaborate entro 2 ore dal prelievo. Le BMMC sono state isolate mediante centrifugazione attraverso un gradiente di Ficoll (Histopaque-1077; Sigma Life Science, St Louis, MO, USA). Le BMMC sono state raccolte dall’interfaccia siero-Ficoll con una pipetta Pasteur monouso e trasferite in una provetta conica da 50 ml con il 2% di siero del paziente in 1 × PBS. Le BMMC sono state contate con un emocitometro e la vitalità è stata valutata con il Trypan blu. Le BMMC sono state risospese in 90% FBS, 10% DMSO media di congelamento e congelate utilizzando un Thermo Scientific Nalgene Mr Frosty (Thermo Scientific) in un congelatore a -80 °C per 24 ore prima di essere trasferite in azoto liquido per la conservazione a lungo termine.

Stima del contenuto di RNA per cellula

La quantità di RNA per tipo di cellula è stata determinata quantificando (Qubit; Invitrogen) l’RNA estratto (Maxwell RSC simplyRNA Cells Kit) da diversi numeri noti di cellule.

Preparazione delle cellule

Le cellule fresche sono state raccolte, lavate con 1 × PBS e risospese a 1 × 106 cellule per ml in 1 × PBS e 0,04% sieroalbumina bovina. Le PBMC fresche sono state congelate a 10 × risospendendo le PBMC in DMEM+40% FBS+10% DMSO, congelando a -by °C in un CoolCell® FTS30 (BioCision) e poi poste in azoto liquido per la conservazione.

Le fiale di cellule congelate da studi ALLCELLS e trapianto sono state scongelate rapidamente in un bagno d’acqua a 37 °C per ∼2 min. Le fiale sono state rimosse quando è rimasto un piccolo cristallo di ghiaccio. Le PBMC scongelate sono state lavate due volte nel mezzo e poi risospese in 1 × PBS e 0,04% sieroalbumina bovina a temperatura ambiente. Le cellule sono state centrifugate a 300 r.c.f. per 5 minuti ogni volta. Le BMMC scongelate sono state lavate e risospese in 1 × PBS e 20% FBS. La concentrazione finale di cellule scongelate era 1 × 106 cellule per ml.

Costruzione di librerie di sequenziamento utilizzando la piattaforma GemCode

Le sospensioni cellulari sono state caricate su un GemCode Single-Cell Instrument (10x Genomics, Pleasanton, CA, USA) per generare GEMs singola cella. Singola cella RNA-Seq librerie sono state preparate utilizzando GemCode Single-Cell 3′ Gel Bead e Library Kit (ora venduto come P / N 120230, 120231, 120232, 10x Genomics). GEM-RT è stata eseguita in un C1000 Touch Thermal cycler con 96-Deep Well Reaction Module (Bio-Rad; P / N 1851197): 55 ° C per 2 ore, 85 ° C per 5 min; tenuto a 4 ° C. Dopo RT, GEMs sono stati rotti e il cDNA a singolo filamento è stato pulito con DynaBeads MyOne Silane Beads (Thermo Fisher Scientific; P/N 37002D) e SPRIselect Reagent Kit (0,6 × SPRI; Beckman Coulter; P/N B23318). cDNA è stato amplificato utilizzando il C1000 Touch termociclatore con 96-Deep Well Reaction Module: 98 °C per 3 min; 14 cicli: 98 °C per 15 s, 67 °C per 20 s e 72 °C per 1 min; 72 °C per 1 min; mantenuto a 4 °C. Il prodotto cDNA amplificato è stato pulito con il kit di reagenti SPRIselect (0,6 × SPRI). Il cDNA è stato successivamente tosato a ∼200 bp utilizzando un sistema Covaris M220 (Covaris; P/N 500295). Le librerie di sequenziamento indicizzate sono state costruite utilizzando i reagenti del kit GemCode Single-Cell 3′ Library, seguendo questi passaggi: (1) riparazione dell’estremità e A-tailing; (2) legatura adattatore; (3) pulizia postligazione con SPRIselect; (4) PCR indice campione e pulizia. Le librerie di sequenziamento del codice a barre sono state quantificate mediante PCR quantitativa (KAPA Biosystems Library Quantification Kit per Illumina piattaforme P/N KK4824). Librerie di sequenziamento sono stati caricati a 2,1 pM su un Illumina NextSeq500 con 2 × 75 kit paired-end utilizzando la seguente lunghezza di lettura: 98 bp Read1, 14 bp I7 Index, 8 bp I5 Index e 10 bp Read2. Alcune librerie precedenti sono stati fatti con 5 nt UMI, e 5 bp Read2 è stato ottenuto invece. Queste librerie sono stati documentati nella tabella supplementare 1.

ERCC test

ERCC sintetico spike-in RNA (Thermo Fisher Scientific; P / N 4456740) sono stati diluiti (1:10 o 1:50) e caricati in un GemCode strumento a cella singola, sostituendo le cellule normalmente utilizzate per generare GEMs. Spike-in Mix1 e Mix2 sono stati entrambi testati. Un protocollo leggermente modificato è stato utilizzato come solo una piccola frazione di GEMs sono stati raccolti per RT e amplificazione cDNA. Dopo il completamento di GEM-RT, 1,25 microlitri di emulsione è stato rimosso e aggiunto a una miscela bifasica di agente di recupero (125 microlitri) (P / N 220016) e 25 mM additivo 1 (30 microlitri) (P / N 220074, 10x Genomics). L’agente di recupero è stato poi rimosso e la soluzione acquosa rimanente è stata pulita con il kit reagente SPRISelect (0,8 × SPRI). cDNA è stato amplificato utilizzando il C1000 Touch Thermal cycler con 96-Deep Well Reaction Module: 98 °C per 3 min; 14 cicli: 98 °C per 15 s, 67 °C per 20 s e 72 °C per 1 min; 72 °C per 1 min; mantenuto a 4 °C. Prodotto cDNA amplificato è stato pulito con il kit reagente SPRIselect (0,8 × ) cDNA è stato successivamente tosato a ∼200 bp utilizzando un sistema Covaris M220 per costruire campioni indicizzati librerie con adattatori 10x Genomics. I conteggi delle molecole ERCC attesi sono stati calcolati in base alla quantità di molecole ERCC utilizzate e ai fattori di diluizione del campione. I conteggi sono stati confrontati con i conteggi delle molecole rilevate (conteggi UMI) per calcolare l’efficienza di conversione.

saggio ddPCR

Le cellule Jurkat sono state utilizzate nei saggi ddPCR per stimare l’efficienza di conversione come segue: (1) la quantità di RNA per cellula Jurkat è stata determinata quantificando (Qubit, Invitrogen) l’RNA estratto (Maxwell RNA Purification Kits) da un diverso numero noto di cellule Jurkat. (2) La RT-ddPCR di massa (Bio-Rad One-Step RT-ddPCR Advanced Kit for Probes 1864021) è stata eseguita sull’RNA estratto per determinare il numero di copie per cella di otto geni selezionati. (3) Circa 5.000 cellule Jurkat sono state elaborate utilizzando la piattaforma GemCode Single-Cell 3′, e il cDNA a singolo filamento è stato raccolto dopo RT in GEMs seguendo i protocolli elencati nella sezione ‘Sequencing library construction using the GemCode platform’. Le copie di cDNA degli otto geni sono state determinate utilizzando ddPCR (Bio-Rad ddPCR Supermix for Probes (no dUTP) P/N 1863024). Il conteggio effettivo delle cellule Jurkat è stato trovato sequenziando un sottoinsieme delle reazioni GEM-RT su un MiSeq. L’efficienza di conversione è il rapporto tra le copie di cDNA per cella (passo 3) e le copie di RNA per cella da massa RT-ddPCR (passo 2), assumendo un 50% di efficienza in RT-ddPCR34.

Le sequenze di sonde per il test ddPCR sono i seguenti: SERAC1_f, 5′-CACGAGCCGCCAGC-3′ e SERAC1_r, 5′-TCTGCAACAGATGACGCAATAAG-3′; SERAC1_p: /56-FAM/CGCCTGCCG/ZEN/GCAGAATGTC/3IABkFQ/. AP1S3_f, 5′-GAAGCAGCCATGGTCTAAGC-3′ e AP1S3_r, 5′-CCTTGTCGACTGAAGCAATATG-3′; AP1S3_p: /56-FAM/CGGCCCAGC/ZEN/CACGATGATACAT/3IABkFQ/OR. AOV1_f, 5′-CCGGAAGTGGGTCGTOR-3′ e AOV1_r, 5′-TTCTTCATAGCCTTCCCGATACCOR-3′; AOV1_p: /56-FAM/TCGTGATGG/ZEN/CGGATGAGAGGTTTCA/3IABkFQ/. DOLPP1_f, 5′-ATGGCAGCGGACGGA-3′ e DOLPP1_r, 5′-GGCTCAGGTAGGCAAGGA-3′; DOLPP1_p: /56-FAM/CCACGTCGA/ZEN/ATATCCTGCAGGTGATCT/3IABkFQ/. KPNA6_f, 5′-TGAAAGCTGCCGCTGAAG-3′ e KPNA6_r, 5′-CCCTGGGCTCGCCAT-3′; KPNA6_p: /56-FAM/CGGACCCGC/ZEN/GATGGAGACC/3IABkFQ/. ITSN2_f, 5′-GTGACAGGCTACGCAACAG-3′ e ITSN2_r, 5′-TCCTGAGTTTTCCTTGCTAGCT-3′; ITSN2_p: /56-FAM/AGGGCGGCCA/ZEN/GATGGCTGA/3IABkFQ/. LCMT1_f, 5′-GTCGACCCCGCTTCCA-3′ e LCMT1_r, 5′-GGTCATGCCAGTAGCCAATG-3′; LCMT1_p: /56-FAM/ATGCTTCCC/ZEN/TGTGCAAGAGGTTTGC/3IABkFQ/. AP2M1_f, 5′-GCAGCGGGCAGACG-3′ e AP2M1_r, 5′-ATGGCGGCAGATCAGTCT-3′; AP2M1_p: /56-FAM/CATCGCTCT/ZEN/GAGAACAGACCTGGTG/3IABkFQ/.

Calcolo dell’efficienza di cattura delle cellule

L’efficienza è calcolata prendendo il rapporto tra il numero di cellule rilevate dal sequenziamento rispetto al numero di cellule caricate nel chip. Quest’ultimo è determinato da (volume aggiunto × concentrazione di ingresso delle cellule). La concentrazione in ingresso delle cellule è stata determinata utilizzando un contatore cellulare automatizzato Countess II (Thermo Fisher Scientific). Vale la pena notare che c’è un errore del 15-20% nel conteggio delle cellule, che potrebbe spiegare almeno una parte della variabilità nelle efficienze calcolate.

Saggio di chimerismo

PowerPlex 16 System (Promega) è stato utilizzato in combinazione con un Applied Biosystems (Life Technologies) 3130xl analizzatore genetico. Donatore BMMCs sono stati utilizzati come la linea di base di riferimento.

Allineamento, assegnazione codice a barre e conteggio UMI

Il Cell Ranger Single-Cell Software Suite è stato utilizzato per eseguire demultiplexing campione, elaborazione codice a barre e singola cella 3′ gene conteggio (http://software.10xgenomics.com/single-cell/overview/welcome). In primo luogo, demultiplexing campione è stato eseguito sulla base del campione 8 bp indice letto per generare FASTQs per il Read1 e Read2 paired-end legge, così come il codice a barre GemCode 14 bp. Dieci tag UMI basepair sono stati estratti da Read2 (14 librerie sono state fatte con 5 bp tag UMI, come indicato nella tabella 1 supplementare, a causa di una precedente iterazione dei metodi. Per questi campioni, i tag UMI da 5 bp sono stati estratti da Read2). Poi, Read1, che contiene l’inserto cDNA, è stato allineato a un genoma di riferimento appropriato utilizzando STAR35. Per le cellule di topo, è stato usato mm10 e per le cellule umane, hg19. Per i campioni con miscele di cellule murine e umane, è stata usata l’unione di hg19 e mm10. Per i campioni ERCC, è stato utilizzato il riferimento ERCC (https://tools.thermofisher.com/content/sfs/manuals/cms_095047.txt).

In seguito, sono stati filtrati i codici a barre GemCode e gli UMI. Sono stati considerati tutti gli elenchi noti di codici a barre che si trovano a 1 distanza di Hamming da un codice a barre osservato. Poi, viene calcolata la probabilità posteriore che il codice a barre osservato sia stato prodotto da un errore di sequenziamento, date le qualità di base del codice a barre osservato e la probabilità anteriore di osservare il codice a barre candidato (presa dalla distribuzione generale del conteggio dei codici a barre). Se la probabilità posteriore per qualsiasi codice a barre candidato è almeno 0,975, allora il codice a barre viene corretto al codice a barre candidato con la più alta probabilità posteriore. Se tutte le sequenze candidate sono ugualmente probabili, allora viene scelta quella che appare per prima in ordine lessicale.

UMI con punteggio di qualità di sequenziamento >10 sono stati considerati validi se non erano omopolimeri. Qual=10 implica un’accuratezza di chiamata della base del 90%. Un UMI che è 1-Hamming-distante da un altro UMI (con più letture) per lo stesso codice a barre e lo stesso gene è corretto all’UMI con più letture. Questo approccio è quasi identico a quello di Jaitin et al.4, ed è simile a quello di Klein et al.8 (anche se Klein et al.8 ha anche usato gli UMI per risolvere le letture multimappate, che non è stato implementato qui).

Infine, i duplicati PCR sono stati contrassegnati se due serie di coppie di letture condividevano una sequenza di codice a barre, un tag UMI e un ID del gene (sono stati utilizzati gli Ensembl GTF GRCh37.82, ftp://ftp.ensembl.org/pub/grch37/release-84/gtf/homo_sapiens/Homo_sapiens.GRCh37.82.gtf.gz e GRCm38.84, ftp://ftp.ensembl.org/pub/release-84/gtf/mus_musculus/Mus_musculus.GRCm38.84.gtf.gz). Solo i duplicati mappati con sicurezza (MAPQ=255), non-PCR con codici a barre e UMI validi sono stati utilizzati per generare la matrice gene-barcode.

I codici a barre delle cellule sono stati determinati in base alla distribuzione dei conteggi UMI. Tutti i codici a barre superiori entro lo stesso ordine di grandezza (>10% del codice a barre superiore n, dove n è l’1% del conteggio delle cellule recuperate previsto) sono stati considerati codici a barre delle cellule. Il numero di letture che forniscono informazioni significative è calcolato come il prodotto di quattro metriche: (1) codici a barre validi; (2) UMI valido; (3) associato a un codice a barre delle cellule; e (4) mappato con sicurezza agli esoni.

Negli esperimenti di miscelazione del topo e dell’uomo, il tasso di multiplet è stato definito come il doppio del tasso di codici a barre delle cellule con conteggi UMI significativi sia dal topo che dall’uomo, dove il top 1% dei conteggi UMI era considerato significativo. L’entità del crosstalk codice a barre è stato valutato dalla frazione di mouse legge in codici a barre umani, o viceversa.

I campioni elaborati da più canali possono essere combinati concatenando matrici gene-cella-barcode. Questa funzionalità è fornita nel kit Cell Ranger R (http://support.10xgenomics.com/single-cell/software/pipelines/latest/rkit). I dati di sequenziamento da più corse di sequenziamento di una libreria possono essere combinati contando le letture non duplicate. Questa funzionalità è fornita nella pipeline Cell Ranger. Inoltre, i dati di sequenziamento possono essere sottocampionati per ottenere un determinato numero di conteggi UMI per cella. Questa funzionalità è anche fornita nel Cell Ranger R Kit, ed è utile quando si combinano i dati da più campioni per il confronto.

Analisi PCA di miscelazione di cellule Jurkat e 293T

Matrice gene-cellula-barcode da ciascuno dei quattro campioni è stato concatenato. Sono stati utilizzati solo i geni con almeno un conteggio UMI rilevato in almeno una cella. La normalizzazione UMI è stata eseguita dividendo prima i conteggi UMI per i conteggi totali UMI in ogni cella, seguita dalla moltiplicazione con la mediana dei conteggi totali UMI in tutte le cellule. Poi, abbiamo preso il log naturale dei conteggi UMI. Infine, ogni gene è stato normalizzato in modo che il segnale medio per ogni gene sia 0 e la deviazione standard sia 1. La PCA è stata eseguita sulla matrice gene-barcode normalizzata. I conteggi UMI normalizzati di ogni gene sono utilizzati per mostrare l’espressione di un marcatore in una trama tSNE.

Analisi SNV di Jurkat e 293T scRNA-seq dati

SNVs sono stati chiamati eseguendo Freebayes 1.0.2 (rif. 36) sul genoma BAM prodotto da Cell Ranger. SNV di alta qualità (SNV chiamando Qual>=100 con almeno 10 conteggi UMI da almeno due cellule; gli indel sono ignorati) che sono stati osservati solo in cellule Jurkat o 293T (ma non entrambi) sono stati selezionati. Le cellule sono state etichettate come Jurkat o 293T sulla base di Jurkat- e 293T-specifici conteggi SNV, dove la frazione di conteggi da altre specie è <0,2. Le cellule con una frazione di SNV da entrambe le specie tra 0,2 e 0,8 sono considerate multiplette. Il tasso di multipletto dedotto è 2* tasso di multipletto osservato (per tenere conto dei multipletti Jurkat:Jurkat e 293T:293T).

PCA e analisi tSNE di PBMCs

Sono usati geni con almeno un conteggio UMI rilevato in almeno una cellula. I primi 1.000 geni più variabili sono stati identificati in base alla loro media e dispersione (varianza/media), che è simile all’approccio utilizzato da Macoscko et al.7 I geni sono stati inseriti in 20 bins in base alla loro espressione media. Dispersione normalizzata è calcolato come la differenza assoluta tra dispersione e dispersione mediana della media di espressione, normalizzata dalla deviazione assoluta mediana all’interno di ogni bin.

PCA è stato eseguito sulla matrice normalizzata gene-barcode dei primi 1.000 geni più variabili per ridurre il numero di caratteristiche (gene) dimensioni. La normalizzazione degli UMI è stata eseguita dividendo prima i conteggi degli UMI per i conteggi totali degli UMI in ogni cella, seguita dalla moltiplicazione con la mediana dei conteggi totali degli UMI nelle celle. Poi, abbiamo preso il log naturale dei conteggi UMI. Infine, ogni gene è stato normalizzato in modo che il segnale medio per ogni gene sia 0 e la deviazione standard sia 1. La PCA è stata eseguita sulla matrice gene-barcode normalizzata. Dopo aver eseguito la PCA, Barnes-hut37 approssimazione a t-SNE16 è stata eseguita sui primi 50 PC per visualizzare le cellule in uno spazio bidimensionale. Cinquanta PC sono stati utilizzati perché: (1) usare tutte le PC richiederebbe molto tempo con l’analisi tSNE; (2) hanno spiegato ∼25% della varianza totale. K-means15 clustering è stato eseguito per raggruppare le cellule per l’analisi di clustering. k = 10 è stato selezionato sulla base della somma degli errori quadratici scree plot (Supplementary Fig. 5d).

Identificazione dei geni specifici del cluster e classificazione basata su marcatori

Per identificare i geni che sono arricchiti in un cluster specifico, l’espressione media di ogni gene è stata calcolata in tutte le cellule del cluster. Poi ogni gene del cluster è stato confrontato con l’espressione mediana dello stesso gene nelle cellule di tutti gli altri cluster. I geni sono stati classificati in base alla loro differenza di espressione, e sono stati selezionati i primi 10 geni arricchiti da ogni cluster. Per il clustering gerarchico, è stata calcolata la correlazione a coppie tra ogni cluster, e l’espressione centrata di ogni gene è stata utilizzata per la visualizzazione tramite heatmap.

La classificazione delle PBMC è stata dedotta dall’annotazione dei geni specifici del cluster. Nel caso del cluster 10, è stata rilevata l’espressione del marcatore di più tipi di cellule (per esempio, B, dendritiche e T). Poiché la dimensione relativa dei cluster di B, dendritici e T è 5,7%, 6,6% e 81%, rispettivamente, ci aspetteremmo che il cluster 10 (che è solo 0,5%) contenga multipli composti principalmente da B:dendritici (0.36%) e B:dendritico:T (0,3%).

Selezione di sottopopolazioni purificate di PBMC

Ogni popolazione di PBMC purificata è stata sottocampionata a ∼16k letture per cella. PCA, tSNE e k-means clustering sono stati eseguiti per ogni matrice downsampled, seguendo gli stessi passaggi delineati in PCA e t-SNE analisi di PBMCs. Solo un cluster è stato rilevato nella maggior parte dei campioni, coerente con le analisi FACS (Fig. 6 supplementare). Per i campioni con più di un cluster, solo i cluster che hanno visualizzato l’espressione del gene marcatore previsto sono stati selezionati per l’analisi a valle. Per i monociti CD14+, due cluster sono stati osservati e identificati come monociti CD14+ e cellule dendritiche in base all’espressione dei geni marcatori FTL e CLEC9A, rispettivamente.

Analisi di classificazione delle cellule utilizzando PBMC purificate

Ogni popolazione di PBMC purificate è stata sottocampionata a ∼16k letture mappate con sicurezza per cella. Poi, è stato calcolato un profilo di espressione genica medio per tutte le cellule. Successivamente, l’espressione genica di ogni cellula della popolazione complessa è stata confrontata con i profili di espressione genica di popolazioni purificate di PBMC mediante la correlazione di Spearman. Alla cellula è stato assegnato l’ID della popolazione purificata se aveva la più alta correlazione con quella popolazione. Si noti che la differenza tra la più alta e la seconda più alta correlazione era piccola per alcune cellule (per esempio, la differenza tra le cellule T citotossiche e NK), suggerendo che l’assegnazione delle cellule non era così sicura per queste cellule. Alcune delle popolazioni PBMC purificate si sono sovrapposte tra loro. Per esempio, le cellule T-helper CD4+ includono tutte le cellule CD4+. Ciò significa che le cellule di questo campione si sovrappongono alle cellule dei campioni che contengono cellule CD4+, tra cui CD4+/CD25+ T reg, CD4+/CD45RO+ T memory, CD4+/CD45RA+/CD25- naive T. Pertanto, quando a una cellula è stato assegnato l’ID di cellula T-helper CD4+ in base al punteggio di correlazione, è stata controllata la correlazione successiva più alta per vedere se era uno dei campioni CD4+. In caso affermativo, l’ID della cellula è stato aggiornato al tipo di cellula con la successiva correlazione più alta. La stessa procedura è stata eseguita per i T citotossici CD8+ e i T citotossici naive CD8+/CD45RA+ (che è un sottoinsieme dei T citotossici CD8+).

Il codice R utilizzato per analizzare 68k PBMC e PBMC purificati può essere trovato qui: https://github.com/10XGenomics/single-cell-3prime-paper.

Classificazione e clustering delle cellule con Seurat

La matrice gene-cellula-barcode di 68k PBMCs è stata log-trasformata come input per Seurat. I primi 469 geni più variabili selezionati da Seurat sono stati utilizzati per calcolare le PC. I primi 22 PC sono risultati significativi (P<0,01) in base all’analisi jackstraw integrata e sono stati utilizzati per la visualizzazione tSNE. La classificazione delle cellule è stata presa dall’analisi di classificazione delle cellule utilizzando PBMC purificate.

Confronto tra PBMC fresche e congelate

I dati di sequenziamento di 68k PBMC fresche e 3k PBMC congelate sono state sottocampionate in modo che ogni campione abbia ∼14k letture mappate con sicurezza per cella. Solo i geni rilevati in almeno una cellula sono stati inclusi per il confronto, che utilizza la media di ogni gene in tutte le cellule.

Per il confronto della classificazione cellulare tra PBMC purificati e congelati, abbiamo raggruppato tutte le cellule etichettate come cellule T o natural killer. Questo perché le sottopopolazioni all’interno di T e tra T e cellule natural killer sono a volte difficili da raggruppare separatamente. Non volevamo che il confronto tra cellule fresche e congelate fosse influenzato dai metodi di clustering utilizzati.

Assegnazione dei genotipi basata sugli SNV

Gli SNV sono stati chiamati eseguendo Freebayes 1.0.2 (rif. 36) sul genoma BAM prodotto da Cell Ranger. Sono stati inclusi solo gli SNV supportati da almeno due codici a barre di cellule, con un punteggio minimo SNV Qual >=30, minimo SNV base Qual>=1. I conteggi degli alleli di riferimento (R) e alternativi (A) sono stati calcolati per ogni SNV, producendo una matrice di conteggi UMI di riferimento delle cellule e conteggi UMI degli alleli alternativi delle cellule. Queste matrici sono state modellate come una miscela di due genomi dove la probabilità di uno qualsiasi dei tre genotipi (R/R, R/A o A/A) in un sito è stata presa per essere distribuita binomialmente con un tasso di errore fisso dello 0,1%. Per ogni campione, due modelli sono stati dedotti in parallelo, uno in cui è presente un solo genoma (K=1) e un altro in cui sono presenti due genomi (K=2). L’inferenza dei parametri del modello (assegnazioni cellula-genoma e i set K di genotipi) è stata eseguita utilizzando un campionatore di Gibbs per approssimare le loro distribuzioni posteriori. Per migliorare il problema del label-switching nell’inferenza Monte Carlo dei modelli di miscela, la rietichettatura delle assegnazioni campionate cellula-genoma è stata eseguita secondo Stephens et al.38

Negli esperimenti di miscelazione delle cellule in silico, quando il modello K = 2 non è riuscito a separare adeguatamente i due genomi, ha riportato una distribuzione di probabilità posteriore vicino a 0,5 per le chiamate cellula-genoma, indicando una mancanza di fiducia in quelle chiamate. Abbiamo applicato il requisito che il 90% delle cellule abbia una probabilità posteriore >75% per selezionare il modello K=2 rispetto al modello K=1. Selezionare K=1 indica che la frazione di miscela è al di sotto del livello di rilevazione del metodo, che negli esperimenti di miscelazione in silico è stato determinato essere il 4% di 6.000 cellule.

Confronto dei genotipi con il campione puro

Per accertare l’assegnazione dei genotipi agli individui, sono stati considerati solo gli SNV condivisi tra il gruppo di genotipi e il campione puro. Poi, il genotipo medio di tutte le cellule è stato confrontato con quello del campione puro. Per ottenere una linea di base per la % di sovrapposizione dei genotipi tra individui diversi, abbiamo eseguito il confronto a coppie di genotipi chiamati dagli stessi individui (11 confronti a coppie) o da individui diversi (15 confronti a coppie). La sovrapposizione percentuale dei genotipi tra gli stessi individui ha una media di ∼98±0,3%, mentre la sovrapposizione percentuale dei genotipi tra individui diversi ha una media di ∼73±2%.

PCA e analisi tSNE di BMMCs

Sono stati utilizzati i dati di sei campioni: due controlli sani, AML027 pre e post-trapianto e AML035 pre e post-trapianto. Ogni campione è stato ricampionato a ∼10k letture mappate con sicurezza per cella. Poi la matrice del codice a barre gene-cellula di ogni campione è stata concatenata. PCA, tSNE e k-means clustering sono stati eseguiti sulla matrice raggruppata, seguendo gli stessi passaggi delineati in PCA e tSNE analisi di PBMCs. Per il clustering k-means, K=10 è stato utilizzato sulla base della curva nella somma degli errori quadratici scree plot.

I geni cluster-specifici sono stati identificati seguendo i passaggi delineati in ‘Identificazione dei geni cluster-specifici e classificazione basata su marcatori’. La classificazione è stata assegnata in base ai geni specifici del cluster e in base all’espressione di alcuni noti marcatori dei tipi di cellule immunitarie. Blasti e Ery 1 immaturi” si riferisce al cluster 4, che esprime CD34, un marcatore dei progenitori ematopoietici39 , e Gata2, un marcatore degli eritroidi precoci40. Ery 2 immaturo” si riferisce ai cluster 5 e 8, che mostrano l’espressione di Gata1, un fattore di trascrizione essenziale per l’eritropoiesi41, ma non CD71, che si trovano spesso in cellule eritroidi più impegnate39. Ery 3 immaturo” si riferisce al cluster 1, che mostra l’espressione di CD71. Il termine “Ery maturi” si riferisce al cluster 2. HBA1, un marcatore di cellule eritroidi mature, è rilevato preferenzialmente nel cluster 2. Il cluster 3 è stato assegnato come ‘Granulociti immaturi’ a causa dell’espressione di marcatori di granulociti precoci come AZU1 e IL8 (rif. 42), e la mancanza di espressione di CD16. Il cluster 7 è stato assegnato come ‘Monociti’ a causa dell’espressione di CD14 e FCN1, per esempio. ‘B’ si riferisce ai cluster 6 e 9 a causa di marcatori come CD19 e CD79A. ‘T’ si riferisce al cluster 10, a causa di marcatori come CD3D e CD8A.

Disponibilità dei dati

Tutti i dati rilevanti sono disponibili presso gli autori. I dati RNA-seq delle singole cellule sono stati depositati nello Short Read Archive con il numero di adesione SRP073767. I dati sono disponibili anche a http://support.10xgenomics.com/single-cell/datasets. Il codice di analisi per l’analisi 68k PBMC è disponibile a https://github.com/10XGenomics/single-cell-3prime-paper.

.