Massively parallel digital transcriptional profiling of single cells

High-speed imaging of gel beads and cells in GEMs

Ein Mikroskop (Nikon Ti-E, × 10 Objektiv) und eine Hochgeschwindigkeits-Videokamera (Photron SA5, frame rate=4.000 s-1) wurden verwendet, um jedes GEM abzubilden, während es im Mikrofluidik-Chip erzeugt wurde. Eine kundenspezifische Analysesoftware wurde verwendet, um die Anzahl der erzeugten GEMs und die Anzahl der in jedem GEM vorhandenen Beads zu zählen, und zwar auf der Grundlage der Kantenerkennung und des Kontrasts zwischen den Bead-Kanten und den GEM-Kanten und der angrenzenden Flüssigkeit. Die Ergebnisse der Analyse sind in Abb. 1c zusammengefasst. Um die Verteilung der Zellen in den GEMs abzuschätzen, wurden ∼28k Bilder eines Videos auf einer Teilmenge von GEMs manuell gezählt. Die Ergebnisse deuten auf eine annähernde Übereinstimmung mit einer Poisson-Verteilung hin. Der Prozentsatz der mehrfachen Zellverkapselungen war jedoch um 16 % höher als der erwartete Wert, was möglicherweise auf Fehler bei der Unterauswahl oder auf Zell-Zell-Interaktionen zurückzuführen ist (bei der manuellen Zählung wurden einige Zweizell-Klumpen beobachtet) (Ergänzende Abb. 1b).

Zelllinien und Proben von Transplantationspatienten

Jurkat- (ATCC TIB-152), 293T- (ATCC CRL-11268) und 3T3- (ATCC CRL-1658) Zellen wurden von ATCC erworben und gemäß den ATCC-Richtlinien kultiviert. Frische PBMCs, gefrorene PBMCs und BMMCs wurden von ALLCELLS erworben. Gefrorene PBMCs von Spender A wurden aus frischen PBMCs von Spender A durch vorsichtiges Mischen von 1e6 Zellen in Gefriermedium (15 % Dimethylsulfoxid (DMSO) in Iscove’s modifiziertem Dulbecco’s-Medium mit 20 % FBS) hergestellt und in CoolCell FTS30 (BioCision) bei -80 °C mindestens 4 Stunden lang gekühlt, bevor sie zur Lagerung für 3 Wochen in flüssigen Stickstoff überführt wurden.

Das Institutional Review Board des Fred Hutchinson Cancer Research Center genehmigte die Studie mit Transplantatproben. Die angewandten Verfahren standen im Einklang mit der Deklaration von Helsinki von 1975 und der Common Rule. Die Proben wurden entnommen, nachdem die Patienten ihr schriftliches Einverständnis zu molekularen Analysen gegeben hatten. Wir identifizierten Patienten mit AML, die sich einer allogenen hämatopoetischen Stammzelltransplantation am Fred Hutchinson Cancer Research Center unterzogen. Die AML-Diagnose wurde gemäß den überarbeiteten Kriterien der Weltgesundheitsorganisation33 gestellt.

Knochenmarkaspirate wurden 20-30 Tage vor der Transplantation für klinische Standardtests entnommen und nach der Transplantation gemäß dem Behandlungsprotokoll seriell untersucht. Die Aliquote des Knochenmarkaspirats wurden innerhalb von 2 Stunden nach der Entnahme verarbeitet. Die BMMCs wurden durch Zentrifugation über einen Ficoll-Gradienten isoliert (Histopaque-1077; Sigma Life Science, St Louis, MO, USA). Die BMMCs wurden mit einer Einwegpasteurpipette von der Serum-Ficoll-Grenzfläche gesammelt und in ein konisches 50-ml-Röhrchen mit 2 % Patientenserum in 1 × PBS überführt. Die BMMCs wurden mit einem Hämacytometer gezählt und die Lebensfähigkeit wurde mit Trypanblau beurteilt. Die BMMCs wurden in 90% FBS, 10% DMSO Gefriermedium resuspendiert und mit einem Thermo Scientific Nalgene Mr Frosty (Thermo Scientific) in einem -80 °C Gefrierschrank für 24 h eingefroren, bevor sie zur Langzeitlagerung in flüssigen Stickstoff überführt wurden.

Schätzung des RNA-Gehalts pro Zelle

Die RNA-Menge pro Zelltyp wurde durch Quantifizierung (Qubit; Invitrogen) der RNA bestimmt, die mit dem Maxwell RSC simplyRNA Cells Kit aus verschiedenen bekannten Zellzahlen extrahiert wurde.

Zellpräparation

Frische Zellen wurden geerntet, mit 1 × PBS gewaschen und mit 1 × 106 Zellen pro ml in 1 × PBS und 0,04 % Rinderserumalbumin resuspendiert. Frische PBMCs wurden bei 10 × eingefroren, indem PBMCs in DMEM+40% FBS+10% DMSO resuspendiert, in einem CoolCell® FTS30 (BioCision) auf -by °C eingefroren und dann zur Lagerung in flüssigen Stickstoff gelegt wurden.

Gefrorene Zellfläschchen aus ALLCELLS- und Transplantationsstudien wurden in einem Wasserbad bei 37 °C für ∼2 min schnell aufgetaut. Die Fläschchen wurden entfernt, sobald ein winziger Eiskristall übrig war. Aufgetaute PBMCs wurden zweimal mit dem Medium gewaschen und dann in 1 × PBS und 0,04 % Rinderserumalbumin bei Raumtemperatur resuspendiert. Die Zellen wurden bei 300 Umdrehungen pro Minute für jeweils 5 Minuten zentrifugiert. Die aufgetauten BMMCs wurden gewaschen und in 1 × PBS und 20 % FBS resuspendiert. Die Endkonzentration der aufgetauten Zellen betrug 1 × 106 Zellen pro ml.

Aufbau von Sequenzierungsbibliotheken mit der GemCode-Plattform

Die Zellsuspensionen wurden auf ein GemCode Single-Cell Instrument (10x Genomics, Pleasanton, CA, USA) geladen, um Einzelzell-GEMs zu erzeugen. Einzelzell-RNA-Seq-Bibliotheken wurden mit dem GemCode Single-Cell 3′ Gel Bead and Library Kit (jetzt erhältlich als P/N 120230, 120231, 120232, 10x Genomics) hergestellt. Die GEM-RT wurde in einem C1000 Touch Thermocycler mit 96-Deep Well Reaction Module (Bio-Rad; P/N 1851197) durchgeführt: 55 °C für 2 h, 85 °C für 5 min; gehalten bei 4 °C. Nach RT wurden die GEMs gebrochen und die einzelsträngige cDNA mit DynaBeads MyOne Silane Beads (Thermo Fisher Scientific; P/N 37002D) und SPRIselect Reagent Kit (0,6 × SPRI; Beckman Coulter; P/N B23318) aufgereinigt. cDNA wurde mit dem C1000 Touch Thermocycler mit 96-Deep Well Reaction Module amplifiziert: 98 °C für 3 min; 14 × Zyklen: 98 °C für 15 s, 67 °C für 20 s und 72 °C für 1 min; 72 °C für 1 min; gehalten bei 4 °C. Das amplifizierte cDNA-Produkt wurde mit dem SPRIselect Reagent Kit (0,6 × SPRI) aufgereinigt. Die cDNA wurde anschließend mit einem Covaris M220-System (Covaris; P/N 500295) auf ∼200 bp geschert. Indizierte Sequenzierungsbibliotheken wurden mit den Reagenzien des GemCode Single-Cell 3′ Library Kits in den folgenden Schritten erstellt: (1) Endreparatur und A-tailing; (2) Adapterligation; (3) Postligation-Cleanup mit SPRIselect; (4) Probenindex-PCR und Cleanup. Die Barcode-Sequenzierungsbibliotheken wurden mittels quantitativer PCR quantifiziert (KAPA Biosystems Library Quantification Kit for Illumina platforms P/N KK4824). Die Sequenzierbibliotheken wurden mit 2,1 pM auf ein Illumina NextSeq500 mit 2 × 75 Paired-End-Kits geladen, wobei folgende Leselängen verwendet wurden: 98 bp Read1, 14 bp I7 Index, 8 bp I5 Index und 10 bp Read2. Einige frühere Bibliotheken wurden mit 5 nt UMI hergestellt, und stattdessen wurde 5 bp Read2 erhalten. Diese Bibliotheken sind in der ergänzenden Tabelle 1 dokumentiert.

ERCC-Assay

ERCC synthetische Spike-in RNAs (Thermo Fisher Scientific; P/N 4456740) wurden verdünnt (1:10 oder 1:50) und in ein GemCode Single-Cell Instrument geladen, das die normalerweise zur Erzeugung von GEMs verwendeten Zellen ersetzt. Spike-in Mix1 und Mix2 wurden beide getestet. Es wurde ein leicht verändertes Protokoll verwendet, da nur ein kleiner Teil der GEMs für die RT- und cDNA-Amplifikation gesammelt wurde. Nach Abschluss der GEM-RT wurden 1,25 μl der Emulsion entnommen und zu einer zweiphasigen Mischung aus Recovery Agent (125 μl) (P/N 220016) und 25 mM Additiv 1 (30 μl) (P/N 220074, 10x Genomics) hinzugefügt. Anschließend wurde das Recovery Agent entfernt und die verbleibende wässrige Lösung mit dem SPRISelect Reagent Kit (0,8 × SPRI) aufgereinigt. cDNA wurde mit dem C1000 Touch Thermocycler mit 96-Deep Well Reaction Module amplifiziert: 98 °C für 3 min; 14 × Zyklen: 98 °C für 15 s, 67 °C für 20 s und 72 °C für 1 min; 72 °C für 1 min; gehalten bei 4 °C. Das amplifizierte cDNA-Produkt wurde mit dem SPRIselect Reagent Kit (0,8 × ) aufgereinigt. Anschließend wurde die cDNA mit einem Covaris M220-System auf ∼200 bp geschert, um probenindizierte Bibliotheken mit 10x Genomics-Adaptern zu erstellen. Die erwartete Anzahl der ERCC-Moleküle wurde auf der Grundlage der Menge der verwendeten ERCC-Moleküle und der Probenverdünnungsfaktoren berechnet. Die Zählungen wurden mit den Zählungen der detektierten Moleküle (UMI-Zählungen) verglichen, um die Konversionseffizienz zu berechnen.

ddPCR-Assay

Jurkat-Zellen wurden in ddPCR-Assays verwendet, um die Konversionseffizienz wie folgt zu schätzen: (1) Die RNA-Menge pro Jurkat-Zelle wurde durch Quantifizierung (Qubit, Invitrogen) von RNA bestimmt, die aus verschiedenen, bekannten Anzahlen von Jurkat-Zellen extrahiert wurde (Maxwell RNA Purification Kits). (2) Mit der extrahierten RNA wurde eine Bulk-RT-ddPCR (Bio-Rad One-Step RT-ddPCR Advanced Kit for Probes 1864021) durchgeführt, um die Kopienzahl pro Zelle von acht ausgewählten Genen zu bestimmen. (3) Etwa 5.000 Jurkat-Zellen wurden mit der GemCode Single-Cell 3′-Plattform prozessiert, und einzelsträngige cDNA wurde nach RT in GEMs gemäß den im Abschnitt „Aufbau von Sequenzierungsbibliotheken mit der GemCode-Plattform“ aufgeführten Protokollen gesammelt. cDNA-Kopien der acht Gene wurden mit ddPCR (Bio-Rad ddPCR Supermix for Probes (no dUTP) P/N 1863024) bestimmt. Die tatsächliche Jurkat-Zellzahl wurde durch Sequenzierung einer Teilmenge der GEM-RT-Reaktionen auf einem MiSeq ermittelt. Die Umwandlungseffizienz ist das Verhältnis zwischen den cDNA-Kopien pro Zelle (Schritt 3) und den RNA-Kopien pro Zelle aus der RT-ddPCR (Schritt 2), wobei eine Effizienz von 50 % bei der RT-ddPCR34 angenommen wird.

Die Sondensequenzen für den ddPCR-Assay lauten wie folgt: SERAC1_f, 5′-CACGAGCCGCCAGC-3′ und SERAC1_r, 5′-TCTGCAACAGATGACGCAATAAG-3′; SERAC1_p: /56-FAM/CGCCTGCCG/ZEN/GCAGAATGTC/3IABkFQ/. AP1S3_f, 5′-GAAGCAGCCATGGTCTAAGC-3′ und AP1S3_r, 5′-CCTTGTCGACTGAAGAGCAATATG-3′; AP1S3_p: /56-FAM/CGGCCCAGC/ZEN/CACGATGATACAT/3IABkFQ/OR. AOV1_f, 5′-CCGGAAGTGGGTCTCGTOR-3′ und AOV1_r, 5′-TTCTTCATAGCCTTCCCGATACCOR-3′; AOV1_p: /56-FAM/TCGTGATGG/ZEN/CGGATGAGAGGTTTCA/3IABkFQ/. DOLPP1_f, 5′-ATGGCAGCGGACGGA-3′ und DOLPP1_r, 5′-GGCTCAGGTAGGCAAGGA-3′; DOLPP1_p: /56-FAM/CCACGTCGA/ZEN/ATATCCTGCAGGTGATCT/3IABkFQ/. KPNA6_f, 5′-TGAAAGCTGCCGCTGAAG-3′ und KPNA6_r, 5′-CCCTGGGCTCGCCAT-3′; KPNA6_p: /56-FAM/CGGACCCGC/ZEN/GATGGAGACC/3IABkFQ/. ITSN2_f, 5′-GTGACAGGCTACGCAACAG-3′ und ITSN2_r, 5′-TCCTGAGTTTTCCTTGCTAGCT-3′; ITSN2_p: /56-FAM/AGGGCGCCA/ZEN/GATGGCTGA/3IABkFQ/. LCMT1_f, 5′-GTCGACCCCGCTTCCA-3′ und LCMT1_r, 5′-GGTCATGCCAGTAGCCAATG-3′; LCMT1_p: /56-FAM/ATGCTTCCC/ZEN/TGTGCAAGAGGTTTGC/3IABkFQ/. AP2M1_f, 5′-GCAGCGGGCAGACG-3′ und AP2M1_r, 5′-ATGGCGGCAGATCAGTCT-3′; AP2M1_p: /56-FAM/CATCGCTCT/ZEN/GAGAACAGACCTGGTG/3IABkFQ/.

Berechnung der Zellfangeffizienz

Die Effizienz wird berechnet, indem das Verhältnis zwischen der Anzahl der durch Sequenzierung nachgewiesenen Zellen und der Anzahl der in den Chip geladenen Zellen ermittelt wird. Letzteres ergibt sich aus (zugegebenes Volumen × Eingangskonzentration der Zellen). Die Eingangskonzentration der Zellen wurde mit einem Countess II Automated Cell Counter (Thermo Fisher Scientific) bestimmt. Es sei darauf hingewiesen, dass bei der Zellzählung ein Fehler von 15-20 % auftritt, der zumindest einen Teil der Schwankungen bei den berechneten Wirkungsgraden erklären könnte.

Chimärismus-Assay

Das PowerPlex 16 System (Promega) wurde in Verbindung mit einem Applied Biosystems (Life Technologies) 3130xl Genetic Analyzer verwendet. Spender-BMMCs wurden als Referenz-Basislinie verwendet.

Ausrichtung, Barcode-Zuordnung und UMI-Zählung

Die Cell Ranger Single-Cell Software Suite wurde verwendet, um Proben-Demultiplexing, Barcode-Verarbeitung und Einzelzell-3′-Gen-Zählung (http://software.10xgenomics.com/single-cell/overview/welcome) durchzuführen. Zunächst wurde das Demultiplexen der Proben auf der Grundlage des 8-Bp-Probenindex-Reads durchgeführt, um FASTQs für die Read1- und Read2-Paired-End-Reads sowie den 14-Bp-GemCode-Barcode zu erzeugen. Zehn Basenpaare UMI-Tags wurden aus Read2 extrahiert (14 Bibliotheken wurden aufgrund einer früheren Iteration der Methoden mit 5 bp UMI-Tags erstellt, wie in der ergänzenden Tabelle 1 angegeben). Für diese Proben wurden 5 bp UMI-Tags aus Read2 extrahiert.). Anschließend wurde Read1, das das cDNA-Insert enthält, mit STAR35 an ein geeignetes Referenzgenom angeglichen. Für Mauszellen wurde mm10 und für menschliche Zellen hg19 verwendet. Für Proben mit Mischungen aus Maus- und Humanzellen wurde die Vereinigung von hg19 und mm10 verwendet. Für ERCC-Proben wurde die ERCC-Referenz (https://tools.thermofisher.com/content/sfs/manuals/cms_095047.txt) verwendet.

Als nächstes wurden GemCode-Barcodes und UMIs gefiltert. Alle bekannten Barcodes, die eine 1-Hamming-Distanz von einem beobachteten Barcode entfernt sind, werden berücksichtigt. Dann wird die Posterior-Wahrscheinlichkeit berechnet, dass der beobachtete Barcode durch einen Sequenzierungsfehler erzeugt wurde, wenn man die Basisqualitäten des beobachteten Barcodes und die Prior-Wahrscheinlichkeit für die Beobachtung des Kandidatenbarcodes (aus der Gesamtverteilung der Barcodeanzahl) berücksichtigt. Wenn die Posterior-Wahrscheinlichkeit für einen beliebigen Kandidatenbarcode mindestens 0,975 beträgt, wird der Barcode auf den Kandidatenbarcode mit der höchsten Posterior-Wahrscheinlichkeit korrigiert. Wenn alle Kandidatensequenzen gleich wahrscheinlich sind, wird diejenige ausgewählt, die in lexikalischer Reihenfolge zuerst auftaucht.

UMIs mit einem Sequenzierungsqualitätsscore >10 wurden als gültig betrachtet, wenn sie keine Homopolymere waren. Qual=10 bedeutet 90 % Basenaufrufgenauigkeit. Eine UMI, die 1-Hamming-Distanz von einer anderen UMI (mit mehr Reads) für denselben Zell-Barcode und dasselbe Gen entfernt ist, wird auf die UMI mit mehr Reads korrigiert. Dieser Ansatz ist nahezu identisch mit dem von Jaitin et al.4 und ähnelt dem von Klein et al.8 (obwohl Klein et al.8 auch UMIs zur Auflösung von Multimapped Reads verwendeten, was hier nicht implementiert wurde).

Schließlich wurden PCR-Duplikate markiert, wenn zwei Sätze von Lesepaaren eine Barcode-Sequenz, ein UMI-Tag und eine Gen-ID gemeinsam hatten (die Ensembl GTFs GRCh37.82, ftp://ftp.ensembl.org/pub/grch37/release-84/gtf/homo_sapiens/Homo_sapiens.GRCh37.82.gtf.gz und GRCm38.84, ftp://ftp.ensembl.org/pub/release-84/gtf/mus_musculus/Mus_musculus.GRCm38.84.gtf.gz, wurden verwendet). Nur sicher zugeordnete (MAPQ=255), nicht-PCR-Duplikate mit gültigen Barcodes und UMIs wurden zur Erstellung einer Gen-Barcode-Matrix verwendet.

Die Zell-Barcodes wurden anhand der Verteilung der UMI-Zahlen bestimmt. Alle Top-Barcodes innerhalb der gleichen Größenordnung (>10 % des obersten n-ten Barcodes, wobei n 1 % der erwarteten wiederhergestellten Zellzahl ist) wurden als Zell-Barcodes betrachtet. Die Anzahl der Reads, die aussagekräftige Informationen liefern, wird als Produkt von vier Metriken berechnet: (1) gültige Barcodes; (2) gültige UMI; (3) assoziiert mit einem Zell-Barcode; und (4) sicher auf Exons gemappt.

In den Mausexperimenten und den menschlichen Mischexperimenten wurde die Multiplett-Rate als das Doppelte der Rate der Zell-Barcodes mit signifikanten UMI-Zählungen sowohl von der Maus als auch vom Menschen definiert, wobei das oberste 1% der UMI-Zählungen als signifikant angesehen wurde. Das Ausmaß des Barcode-Crosstalk wurde anhand des Anteils der Maus-Reads in den Human-Barcodes oder umgekehrt bewertet.

Proben aus mehreren Kanälen können durch Verkettung von Gen-Zell-Barcode-Matrizen kombiniert werden. Diese Funktion ist im Cell Ranger R Kit (http://support.10xgenomics.com/single-cell/software/pipelines/latest/rkit) enthalten. Sequenzierungsdaten aus mehreren Sequenzierungsläufen einer Bibliothek können durch Zählen nicht-duplizierter Reads kombiniert werden. Diese Funktion ist in der Cell Ranger-Pipeline enthalten. Darüber hinaus können Sequenzierungsdaten unterabgetastet werden, um eine bestimmte Anzahl von UMI-Zählungen pro Zelle zu erhalten. Diese Funktion ist ebenfalls im Cell Ranger R-Kit enthalten und ist nützlich, wenn Daten aus mehreren Proben zum Vergleich kombiniert werden sollen.

PCA-Analyse der Mischung von Jurkat- und 293T-Zellen

Die Gen-Zell-Barcode-Matrix aus jeder der vier Proben wurde verkettet. Es wurden nur Gene verwendet, bei denen mindestens ein UMI-Wert in mindestens einer Zelle nachgewiesen wurde. Die UMI-Normalisierung wurde durchgeführt, indem zunächst die UMI-Zählungen durch die gesamten UMI-Zählungen in jeder Zelle geteilt wurden, gefolgt von einer Multiplikation mit dem Median der gesamten UMI-Zählungen in allen Zellen. Dann wurde der natürliche Logarithmus der UMI-Zahlen genommen. Schließlich wurde jedes Gen so normalisiert, dass das mittlere Signal für jedes Gen 0 und die Standardabweichung 1 ist. Die PCA wurde mit der normalisierten Gen-Barcode-Matrix durchgeführt. Die normalisierten UMI-Zahlen jedes Gens werden verwendet, um die Expression eines Markers in einem tSNE-Plot darzustellen.

SNV-Analyse von Jurkat- und 293T-scRNA-seq-Daten

SNVs wurden durch Ausführen von Freebayes 1.0.2 (Ref. 36) auf dem von Cell Ranger erstellten Genom-BAM aufgerufen. Hochwertige SNVs (SNV calling Qual>=100 mit mindestens 10 UMI-Zählungen aus mindestens zwei Zellen; Indels werden ignoriert), die nur in Jurkat- oder 293T-Zellen (aber nicht in beiden) beobachtet wurden, wurden ausgewählt. Die Zellen wurden auf der Grundlage der Jurkat- und 293T-spezifischen SNV-Zählungen als Jurkat oder 293T gekennzeichnet, wobei der Anteil der Zählungen der anderen Spezies <0,2 beträgt. Zellen mit einem Anteil an SNV von einer der beiden Spezies zwischen 0,2 und 0,8 werden als Multipletts betrachtet. Die abgeleitete Multiplett-Rate ist 2* beobachtete Multiplett-Rate (um Jurkat:Jurkat und 293T:293T-Multipletts zu berücksichtigen).

PCA- und tSNE-Analyse von PBMCs

Gene mit mindestens einer UMI-Zählung, die in mindestens einer Zelle nachgewiesen wurde, werden verwendet. Die 1.000 variabelsten Gene wurden auf der Grundlage ihres Mittelwerts und ihrer Streuung (Varianz/Mittelwert) identifiziert, was dem von Macoscko et al.7 verwendeten Ansatz ähnelt.7 Die Gene wurden auf der Grundlage ihrer mittleren Expression in 20 Bins eingeteilt. Die normalisierte Streuung wird als absolute Differenz zwischen der Streuung und der medianen Streuung des Mittelwerts der Expression berechnet, normalisiert durch die mediane absolute Abweichung innerhalb jedes Bins.

PCA wurde auf der normalisierten Gen-Barcode-Matrix der 1.000 variabelsten Gene durchgeführt, um die Anzahl der Merkmalsdimensionen (Gene) zu reduzieren. Die UMI-Normalisierung wurde durchgeführt, indem zunächst die UMI-Zahlen durch die gesamten UMI-Zahlen in jeder Zelle geteilt wurden, gefolgt von einer Multiplikation mit dem Median der gesamten UMI-Zahlen in den Zellen. Dann wurde der natürliche Logarithmus der UMI-Zahlen genommen. Schließlich wurde jedes Gen so normalisiert, dass das mittlere Signal für jedes Gen 0 und die Standardabweichung 1 ist. Die PCA wurde auf der normalisierten Gen-Barcode-Matrix durchgeführt. Nach der Durchführung der PCA wurde die Barnes-Hut37-Näherung an t-SNE16 für die ersten 50 PCs durchgeführt, um die Zellen in einem zweidimensionalen Raum zu visualisieren. Es wurden fünfzig PCs verwendet, weil: (1) die Verwendung aller PCs würde bei der tSNE-Analyse sehr viel Zeit in Anspruch nehmen; (2) sie erklärten ∼25 % der Gesamtvarianz. K-means15 Clustering wurde durchgeführt, um die Zellen für die Clusteranalyse zu gruppieren. k=10 wurde auf der Grundlage der Summe der quadrierten Fehler im Scree Plot ausgewählt (ergänzende Abb. 5d).

Identifizierung von clusterspezifischen Genen und markerbasierte Klassifizierung

Um Gene zu identifizieren, die in einem bestimmten Cluster angereichert sind, wurde die mittlere Expression jedes Gens über alle Zellen im Cluster berechnet. Dann wurde jedes Gen aus dem Cluster mit der mittleren Expression desselben Gens aus Zellen in allen anderen Clustern verglichen. Die Gene wurden auf der Grundlage ihrer Expressionsdifferenz in eine Rangfolge gebracht, und die 10 am stärksten angereicherten Gene aus jedem Cluster wurden ausgewählt. Für das hierarchische Clustering wurde die paarweise Korrelation zwischen den einzelnen Clustern berechnet, und die zentrierte Expression jedes Gens wurde zur Visualisierung mittels Heatmap verwendet.

Die Klassifizierung der PBMCs wurde aus der Annotation der clusterspezifischen Gene abgeleitet. Im Fall von Cluster 10 wurde die Expression von Markern mehrerer Zelltypen (z. B. B, dendritisch und T) festgestellt. Da die relative Clustergröße von B, dendritisch und T 5,7 %, 6,6 % bzw. 81 % beträgt, würden wir erwarten, dass Cluster 10 (das nur 0,5 % ausmacht) Multiplets enthält, die hauptsächlich aus B:dendritisch (0.36%) und B:dendritisch:T (0,3%) besteht.

Auswahl gereinigter Subpopulationen von PBMCs

Jede Population gereinigter PBMCs wurde auf ∼16k Reads pro Zelle heruntergestuft. PCA, tSNE und k-means Clustering wurden für jede downsampled Matrix durchgeführt, wobei die gleichen Schritte wie bei der PCA- und t-SNE-Analyse von PBMCs befolgt wurden. In den meisten Proben wurde nur ein Cluster entdeckt, was mit den FACS-Analysen übereinstimmt (ergänzende Abb. 6). Bei Proben mit mehr als einem Cluster wurden nur die Cluster, die die erwartete Expression der Markergene aufwiesen, für die nachgeschaltete Analyse ausgewählt. Bei CD14+ Monozyten wurden zwei Cluster beobachtet und anhand der Expression der Markergene FTL und CLEC9A als CD14+ Monozyten bzw. dendritische Zellen identifiziert.

Zellklassifizierungsanalyse mit gereinigten PBMCs

Jede Population gereinigter PBMCs wurde auf ∼16k sicher zugeordnete Reads pro Zelle heruntergerechnet. Dann wurde ein durchschnittliches (mittleres) Genexpressionsprofil für alle Zellen berechnet. Anschließend wurde die Genexpression jeder Zelle der komplexen Population mit den Genexpressionsprofilen gereinigter Populationen von PBMCs durch die Spearman-Korrelation verglichen. Der Zelle wurde die ID der gereinigten Population zugewiesen, wenn sie die höchste Korrelation mit dieser Population aufwies. Beachten Sie, dass der Unterschied zwischen der höchsten und der zweithöchsten Korrelation bei einigen Zellen gering war (z. B. der Unterschied zwischen zytotoxischen T- und NK-Zellen), was darauf hindeutet, dass die Zellzuordnung bei diesen Zellen nicht so sicher war. Einige der gereinigten PBMC-Populationen überlappten sich. Die CD4+ T-Helferzellen umfassen beispielsweise alle CD4+ Zellen. Das bedeutet, dass sich Zellen aus dieser Probe mit Zellen aus Proben überschneiden, die CD4+ Zellen enthalten, einschließlich CD4+/CD25+ T reg, CD4+/CD45RO+ T memory, CD4+/CD45RA+/CD25- naive T. Wenn also einer Zelle auf der Grundlage des Korrelationswertes die ID CD4+ T-Helferzelle zugewiesen wurde, wurde die nächsthöhere Korrelation überprüft, um festzustellen, ob es sich um eine der CD4+ Proben handelte. War dies der Fall, wurde die ID der Zelle auf den Zelltyp mit der nächsthöheren Korrelation aktualisiert. Das gleiche Verfahren wurde für CD8+ zytotoxische T-Zellen und CD8+/CD45RA+ naive zytotoxische T-Zellen (eine Untergruppe der CD8+ zytotoxischen T-Zellen) durchgeführt.

Den R-Code, der zur Analyse von 68k PBMCs und gereinigten PBMCs verwendet wurde, finden Sie hier: https://github.com/10XGenomics/single-cell-3prime-paper.

Zellclusterung und Klassifizierung mit Seurat

Die Gen-Zell-Barcode-Matrix von 68k PBMCs wurde als Eingabe für Seurat log-transformiert. Die 469 variabelsten Gene, die von Seurat ausgewählt wurden, wurden für die Berechnung der PCs verwendet. Die ersten 22 PCs waren signifikant (P<0,01), basierend auf der integrierten Jackstraw-Analyse, und wurden für die tSNE-Visualisierung verwendet. Die Zellklassifizierung wurde von der Zellklassifizierungsanalyse mit gereinigten PBMCs übernommen.

Vergleich zwischen frischen und gefrorenen PBMCs

Die Sequenzierungsdaten von 68k frischen PBMCs und 3k gefrorenen PBMCs wurden so heruntergestuft, dass jede Probe ∼14k sicher zugeordnete Reads pro Zelle aufweist. Für den Vergleich wurden nur Gene berücksichtigt, die in mindestens einer Zelle nachgewiesen wurden, wobei der Mittelwert jedes Gens über alle Zellen hinweg verwendet wurde.

Für den Vergleich der Zellklassifizierung zwischen gereinigten und gefrorenen PBMCs wurden alle Zellen, die als T- oder natürliche Killerzellen gekennzeichnet waren, zusammengefasst. Der Grund dafür ist, dass die Subpopulationen innerhalb der T-Zellen und zwischen T- und natürlichen Killerzellen manchmal schwer getrennt zu klassifizieren sind. Wir wollten nicht, dass der Vergleich zwischen frischen und gefrorenen Zellen durch die verwendeten Clustermethoden beeinträchtigt wird.

SNV-basierte Genotypzuweisung

SNVs wurden aufgerufen, indem Freebayes 1.0.2 (Ref. 36) auf der von Cell Ranger erstellten Genom-BAM ausgeführt wurde. Nur SNVs mit Unterstützung von mindestens zwei Zell-Barcodes, mit einem minimalen SNV-Qual-Score >=30, minimaler SNV-Base-Qual>=1 wurden berücksichtigt. Die Anzahl der Referenzallele (R) und der alternativen Allele (A) wurde für jede SNV berechnet, so dass eine Matrix von Zell-Referenz-UMI-Zahlen und Zell-Alternativallel-UMI-Zahlen entstand. Diese Matrizen wurden als eine Mischung aus zwei Genomen modelliert, wobei die Wahrscheinlichkeit eines der drei Genotypen (R/R, R/A oder A/A) an einer Stelle als binomialverteilt mit einer festen Fehlerrate von 0,1 % angenommen wurde. Für jede Probe wurden zwei Modelle parallel abgeleitet, eines, bei dem nur ein Genom vorhanden ist (K=1), und ein anderes, bei dem zwei Genome vorhanden sind (K=2). Die Ableitung der Modellparameter (Zell-zu-Genom-Zuordnungen und K-Sätze von Genotypen) wurde mit Hilfe eines Gibbs-Samplers durchgeführt, um ihre Posterior-Verteilungen zu approximieren. Um das „Label-Switching“-Problem bei der Monte-Carlo-Inferenz von Mischungsmodellen zu verbessern, wurde eine Neuetikettierung der gesampelten Zell-Genom-Zuordnungen gemäß Stephens et al.38

In in silico-Zellmischungsexperimenten wurde, wenn das K=2-Modell die beiden Genome nicht angemessen trennen konnte, eine Verteilung der Posterior-Wahrscheinlichkeiten in der Nähe von 0,5 für die Zell-Genom-Aufrufe gemeldet, was auf einen Mangel an Vertrauen in diese Aufrufe hinweist. Wir wendeten die Bedingung an, dass 90 % der Zellen eine Posteriorwahrscheinlichkeit von >75 % haben müssen, um das K=2-Modell dem K=1-Modell vorzuziehen. Die Wahl von K=1 zeigt an, dass der Mischungsanteil unter der Nachweisgrenze der Methode liegt, die in In-silico-Mischexperimenten auf 4 % von 6.000 Zellen festgelegt wurde.

Genotypvergleich mit der reinen Probe

Um die Zuordnung der Genotypen zu den Individuen zu überprüfen, wurden nur gemeinsame SNVs zwischen der Genotypengruppe und der reinen Probe berücksichtigt. Anschließend wurde der durchschnittliche Genotyp aller Zellen mit dem der reinen Stichprobe verglichen. Um einen Anhaltspunkt für die prozentuale Genotypüberlappung zwischen verschiedenen Individuen zu erhalten, führten wir paarweise Vergleiche von Genotypen durch, die von denselben Individuen (11 paarweise Vergleiche) oder von verschiedenen Individuen (15 paarweise Vergleiche) stammen. Die prozentuale Genotypüberlappung zwischen denselben Individuen beträgt im Durchschnitt ∼98±0,3 %, während die prozentuale Genotypüberlappung zwischen den verschiedenen Individuen im Durchschnitt ∼73±2 % beträgt.

PCA- und tSNE-Analyse von BMMCs

Daten von sechs Proben wurden verwendet: zwei gesunde Kontrollen, AML027 vor und nach der Transplantation und AML035 vor und nach der Transplantation. Jede Probe wurde auf ∼10k sicher gemappte Reads pro Zelle heruntergestuft. Dann wurde die Gen-Zell-Barcode-Matrix aus jeder Probe verkettet. PCA, tSNE und k-means Clustering wurden mit der gepoolten Matrix durchgeführt, wobei die gleichen Schritte wie bei der PCA- und tSNE-Analyse von PBMCs befolgt wurden. Für das k-means-Clustering wurde K=10 verwendet, basierend auf der Kurve im Sum-of-Squared-Error-Scree-Plot.

Cluster-spezifische Gene wurden gemäß den Schritten identifiziert, die in „Identifizierung cluster-spezifischer Gene und markerbasierte Klassifizierung“ beschrieben sind. Die Klassifizierung erfolgte auf der Grundlage der cluster-spezifischen Gene und auf der Grundlage der Expression einiger bekannter Marker für Immunzellarten. Blasts and Immature Ery 1″ bezieht sich auf Cluster 4, der CD34, einen Marker für hämatopoetische Vorläuferzellen39, und Gata2, einen Marker für frühe Erythroiden40, exprimiert. Unreife Ery 2″ bezieht sich auf die Cluster 5 und 8, die Gata1, einen für die Erythropoese wesentlichen Transkriptionsfaktor, exprimieren41, nicht aber CD71, das häufig in stärker entwickelten erythroiden Zellen zu finden ist39. Unreife Ery 3″ bezieht sich auf Cluster 1, in dem CD71 exprimiert wird. Reife Ery“ bezieht sich auf Cluster 2. HBA1, ein Marker für reife erythroide Zellen, wird bevorzugt in Cluster 2 nachgewiesen. Cluster 3 wurde aufgrund der Expression früher Granulozytenmarker wie AZU1 und IL8 (Ref. 42) und der fehlenden Expression von CD16 als „unreife Granulozyten“ eingestuft. Cluster 7 wurde aufgrund der Expression von CD14 und FCN1 beispielsweise als „Monozyten“ eingestuft. B“ bezieht sich auf die Cluster 6 und 9 aufgrund von Markern wie CD19 und CD79A. T“ bezieht sich auf Cluster 10 aufgrund von Markern wie CD3D und CD8A.

Datenverfügbarkeit

Alle relevanten Daten sind bei den Autoren erhältlich. Einzelzell-RNA-seq-Daten wurden im Short Read Archive unter der Zugriffsnummer SRP073767 hinterlegt. Die Daten sind auch unter http://support.10xgenomics.com/single-cell/datasets verfügbar. Der Analysecode für die 68k PBMC-Analyse ist verfügbar unter https://github.com/10XGenomics/single-cell-3prime-paper.