Massivt parallel digital transkriptionsprofilering af enkeltceller

Højhastighedsafbildning af gelperler og celler i GEM’er

Et mikroskop (Nikon Ti-E, × 10 objektiv) og et højhastigheds-videokamera (Photron SA5, billedhastighed = 4.000 s-1) blev brugt til at afbilde hver enkelt GEM, efterhånden som de blev genereret i den mikrofluidiske chip. Der blev anvendt en tilpasset analysesoftware til at tælle antallet af GEM’er, der blev genereret, og antallet af perler i hver GEM, baseret på kantdetektion og kontrasten mellem perlekanterne og GEM-kanterne og den tilstødende væske. Resultaterne af analysen er opsummeret i fig. 1c. For at estimere fordelingen af celler i GEM’er blev der anvendt manuel optælling for ∼28k frames af en video på en delmængde af GEM’er. Resultaterne viser en tilnærmelsesvis tilslutning til en Poisson-fordeling. Procentdelen af flere celleindkapslinger var imidlertid 16 % højere end den forventede værdi, muligvis på grund af fejl ved subsampling eller på grund af interaktioner mellem celler og celler (nogle klumper af to celler blev observeret under den manuelle optælling) (Supplerende fig. 1b).

Cellelinjer og transplantationspatientprøver

Jurkat (ATCC TIB-152), 293T (ATCC CRL-11268) og 3T3 (ATCC CRL-1658) celler blev erhvervet fra ATCC og dyrket i henhold til ATCC-retningslinjerne. Friske PBMC’er, frosne PBMC’er og BMMC’er blev indkøbt hos ALLCELLS. Frosne PBMC’er fra donor A blev fremstillet af friske PBMC’er fra donor A ved forsigtigt at blande 1e6 celler i frysemedium (15 % dimethylsulfoxid (DMSO) i Iscove’s modificerede Dulbecco-medium indeholdende 20 % FBS) og nedkøles i CoolCell FTS30 (BioCision) ved -80 °C i mindst 4 timer, inden de overføres til flydende nitrogen til opbevaring i 3 uger.

Det institutionelle bedømmelsesudvalg ved Fred Hutchinson Cancer Research Center godkendte undersøgelsen af transplantationsprøver. De fulgte procedurer var i overensstemmelse med Helsinki-erklæringen fra 1975 og den fælles regel. Prøverne blev udtaget, efter at patienterne havde givet skriftligt informeret samtykke til molekylære analyser. Vi identificerede patienter med AML, der gennemgik allogen hæmatopoietisk stamcelletransplantation på Fred Hutchinson Cancer Research Center. Diagnosen AML blev stillet i henhold til Verdenssundhedsorganisationens reviderede kriterier33.

Benmarvsaspirater blev udtaget til klinisk standardtest 20-30 dage før transplantationen og serielt efter transplantationen i henhold til behandlingsprotokollen. Benmarvsaspirataliquots blev behandlet inden for 2 timer efter udtagningen. BMMC’erne blev isoleret ved hjælp af centrifugering gennem en Ficoll-gradient (Histopaque-1077; Sigma Life Science, St Louis, MO, USA). BMMC’erne blev opsamlet fra serum-Ficoll-grænsefladen med en engangs Pasteur-pipette og overført til et 50 ml konisk rør med 2 % patientserum i 1 × PBS. BMMC’erne blev talt ved hjælp af et hæmacytometer, og levedygtigheden blev vurderet ved hjælp af Trypanblåt. BMMC’erne blev resuspenderet i et frysemedie med 90 % FBS og 10 % DMSO og frosset ned i en Thermo Scientific Nalgene Mr Frosty (Thermo Scientific) i en fryser ved -80 °C i 24 timer, inden de blev overført til flydende nitrogen med henblik på langtidsopbevaring.

Vurdering af RNA-indhold pr. celle

Mængden af RNA pr. celletype blev bestemt ved at kvantificere (Qubit; Invitrogen) RNA ekstraheret (Maxwell RSC simplyRNA Cells Kit) fra flere forskellige kendte antal celler.

Celleforberedelse

Friske celler blev høstet, vasket med 1 × PBS og resuspenderet ved 1 × 106 celler pr. ml i 1 × PBS og 0,04 % bovin serumalbumin. Friske PBMC’er blev frosset ned til 10 × ved at resuspendere PBMC’er i DMEM+40% FBS+10% DMSO, fryse til -by °C i en CoolCell® FTS30 (BioCision) og derefter anbringes i flydende nitrogen til opbevaring.

Frosne celleflasker fra ALLCELLS og transplantationsundersøgelser blev hurtigt optøet i et vandbad ved 37 °C i ∼2 min. Beholderne blev fjernet, når der var en lille iskrystal tilbage. Optøede PBMC’er blev vasket to gange i mediet og derefter resuspenderet i 1 × PBS og 0,04 % bovin serumalbumin ved stuetemperatur. Cellerne blev centrifugeret ved 300 r.c.f. i 5 min. hver gang. Optøede BMMC’er blev vasket og resuspenderet i 1 × PBS og 20 % FBS. Den endelige koncentration af optøede celler var 1 × 106 celler pr. ml.

Sekventeringsbibliotekskonstruktion ved hjælp af GemCode-platformen

Cellulære suspensioner blev indlæst på et GemCode Single-Cell Instrument (10x Genomics, Pleasanton, CA, USA) for at generere enkeltcellede GEM’er. RNA-Seq-biblioteker for enkeltceller blev fremstillet ved hjælp af GemCode Single-Cell 3′ Gel Bead and Library Kit (nu solgt som P/N 120230, 120231, 120232, 10x Genomics). GEM-RT blev udført i en C1000 Touch Thermal cycler med 96-Deep Well Reaction Module (Bio-Rad; P/N 1851197): 55 °C i 2 timer, 85 °C i 5 minutter; holdt ved 4 °C. Efter RT blev GEM’erne brudt, og enkeltstrenget cDNA blev renset op med DynaBeads MyOne Silane Beads (Thermo Fisher Scientific; P/N 37002D) og SPRIselect Reagent Kit (0,6 × SPRI; Beckman Coulter; P/N B23318). cDNA blev amplificeret ved hjælp af C1000 Touch Thermal cycler med 96-Deep Well Reaction Module: 98 °C i 3 min; cyklet 14 ×: 98 °C i 15 s, 67 °C i 20 s og 72 °C i 1 min; 72 °C i 1 min; holdt ved 4 °C. Forstærket cDNA-produkt blev renset op med SPRIselect Reagent Kit (0,6 × SPRI). cDNA’et blev efterfølgende skåret til ∼200 bp ved hjælp af et Covaris M220-system (Covaris; P/N 500295). Indekserede sekventeringsbiblioteker blev konstrueret ved hjælp af reagenserne i GemCode Single-Cell 3′ Library Kit, idet følgende trin blev fulgt: (1) reparation af ender og A-tailing; (2) adapterligering; (3) postligeringsoprensning med SPRIselect; (4) prøveindeks-PCR og oprydning. Stregkode-sekventeringsbibliotekerne blev kvantificeret ved kvantitativ PCR (KAPA Biosystems Library Quantification Kit for Illumina-platforme P/N KK4824). Sekventeringsbibliotekerne blev indlæst ved 2,1 pM på en Illumina NextSeq500 med 2 × 75 paired-end-kits med følgende læselængde: 98 bp Read1, 14 bp I7 Index, 8 bp I5 Index og 10 bp Read2. Nogle tidligere biblioteker blev fremstillet med 5 nt UMI, og der blev i stedet opnået 5 bp Read2. Disse biblioteker er dokumenteret i supplerende tabel 1.

ERCC assay

ERCC syntetiske spike-in RNA’er (Thermo Fisher Scientific; P/N 4456740) blev fortyndet (1:10 eller 1:50) og indlæst i et GemCode Single-Cell Instrument, der erstattede celler, der normalt anvendes til at generere GEM’er. Spike-in Mix1 og Mix2 blev begge testet. Der blev anvendt en let ændret protokol, da kun en lille del af GEM’erne blev indsamlet til RT- og cDNA-amplifikation. Efter afslutningen af GEM-RT blev 1,25 μl af emulsionen fjernet og tilsat til en bifasisk blanding af Recovery Agent (125 μl) (P/N 220016) og 25 mM additiv 1 (30 μl) (P/N 220074, 10x Genomics). Genoprettelsesmidlet blev derefter fjernet, og den resterende vandige opløsning blev renset med SPRISelect Reagent Kit (0,8 × SPRI). cDNA blev amplificeret ved hjælp af C1000 Touch Thermal cycler med 96-Deep Well Reaction Module: 98 °C i 3 min; cyklet 14 ×: 98 °C i 15 s, 67 °C i 20 s og 72 °C i 1 min; 72 °C i 1 min; holdt ved 4 °C. Amplificeret cDNA-produkt blev renset op med SPRIselect Reagent Kit (0,8 × ) cDNA blev efterfølgende skåret til ∼200 bp ved hjælp af et Covaris M220-system for at konstruere prøveindekserede biblioteker med 10x Genomics-adaptere. Det forventede antal ERCC-molekyler blev beregnet på grundlag af den anvendte mængde ERCC-molekyler og prøvefortyndingsfaktorer. Antallet blev sammenlignet med detekterede molekyletællinger (UMI-tællinger) for at beregne konverteringseffektiviteten.

ddPCR-assay

Jurkat-celler blev anvendt i ddPCR-assays for at estimere konverteringseffektiviteten på følgende måde: (1) mængden af RNA pr. Jurkat-celle blev bestemt ved at kvantificere (Qubit, Invitrogen) RNA ekstraheret (Maxwell RNA Purification Kits) fra flere forskellige kendte antal Jurkat-celler. (2) Bulk RT-ddPCR (Bio-Rad One-Step RT-ddPCR Advanced Kit for Probes 1864021) blev udført på det ekstraherede RNA for at bestemme antallet af kopier pr. celle af otte udvalgte gener. (3) Ca. 5 000 Jurkat-celler blev behandlet ved hjælp af GemCode Single-Cell 3′-platformen, og enkeltstrenget cDNA blev indsamlet efter RT i GEM’er efter de protokoller, der er anført i afsnittet “Sequencing library construction using the GemCode platform”. cDNA-kopier af de otte gener blev bestemt ved hjælp af ddPCR (Bio-Rad ddPCR Supermix for Probes (no dUTP) P/N 1863024). Det faktiske antal Jurkat-celler blev fundet ved sekventering af en delmængde af GEM-RT-reaktionerne på en MiSeq. Konverteringseffektiviteten er forholdet mellem cDNA-kopier pr. celle (trin 3) og RNA-kopier pr. celle fra bulk-RT-ddPCR (trin 2), idet der antages en effektivitet på 50 % i RT-ddPCR34.

Sondersekvenserne til ddPCR-assayet er som følger: SERAC1_f, 5′-CACGAGAGCCGCCGCCAGC-3′ og SERAC1_r, 5′-TCTGCAACACAGAGATGACGCAATAAG-3′; SERAC1_p: /56-FAM/CGCCCCTGCCG/ZEN/GCAGAATGTC/3IABkFQ/. AP1S3_f, 5′-GAAGCAGCAGCCATATGGTCTCTAAGC-3′ og AP1S3_r, 5′-CCTTGTCGACTGAAGAGAGCAATATATG-3′; AP1S3_p: /56-FAM/CGGCCCCCAGC/ZEN/CACGATGATGATACACAT/3IABkFQ/OR. AOV1_f, 5′-CCGGAAGAGTGGGGGGTCTCTCGTOR-3′ og AOV1_r, 5′-TTCTTCTTCATAGCCTTCTCCCCCGATACCOR-3′; AOV1_p: /56-FAM/TCGTGTGATGG/ZEN/CGGATGAGAGAGAGGTTTCA/3IABkFQ/. DOLPP1_f, 5′-ATGGCAGCAGCGGACGGA-3′ og DOLPP1_r, 5′-GGCTCTCAGGTAGGCAAGGA-3′; DOLPP1_p: /56-FAM/CCACACGTCGA/ZEN/ATATATCCTGCAGCAGGTGTGATCT/3IABkFQ/. KPNA6_f, 5′-TGAAAGAGCTGCCGCCGCTGAAG-3′ og KPNA6_r, 5′-CCCTGGGCTCGTCGCCAT-3′; KPNA6_p: /56-FAM/CGGACCCCCGC/ZEN/GATGGAGAGACC/3IABkFQ/. ITSN2_f, 5′-GTGACACAGAGGCTACGCAACACAG-3′ og ITSN2_r, 5′-TCCTGAGAGTTTTCCTCTTGCTAGCT-3′; ITSN2_p: /56-FAM/AGGGCGCCCA/ZEN/GATGGCTGA/3IABkFQ/. LCMT1_f, 5′-GTCGACCCCCCGCTTCTCCA-3′ og LCMT1_r, 5′-GGTCATGCCAGTAGCCAATG-3′; LCMT1_p: /56-FAM/ATGCTTCTCCCCC/ZEN/TGTGCAAGAGAGGTTTGTGC/3IABkFQ/. AP2M1_f, 5′-GCAGCAGCGGCGGGCAGACG-3′ og AP2M1_r, 5′-ATGGCGGCGGCAGATCAGCAGTCT-3′; AP2M1_p: /56-FAM/CATCGCCTCT/ZEN/GAGAACAGAGACACCTGGGGTG/3IABkFQ/.

Beregning af effektiviteten af celleindfangning

Effektiviteten beregnes ved at tage forholdet mellem antallet af celler, der detekteres ved sekventering, og antallet af celler, der indlæses i chippen. Sidstnævnte bestemmes ud fra (tilsat volumen × indgangskoncentration af celler). Indgangskoncentrationen af celler blev bestemt ved hjælp af en Countess II Automated Cell Counter (Thermo Fisher Scientific). Det er værd at bemærke, at der er en fejl på 15-20 % i celletællinger, hvilket kan forklare i det mindste en del af variabiliteten i de beregnede effektiviteter.

Chimerismeassay

PowerPlex 16 System (Promega) blev anvendt sammen med en Applied Biosystems (Life Technologies) 3130xl Genetic Analyzer. Donor BMMC’er blev anvendt som referencebasislinje.

Alignment, stregkodetildeling og UMI-tælling

Cell Ranger Single-Cell Software Suite blev anvendt til at udføre demultiplexing af prøver, stregkodebehandling og 3′-genoptælling af enkeltceller (http://software.10xgenomics.com/single-cell/overview/welcome). Først blev der udført demultiplexing af prøven på grundlag af prøvindekset på 8 bp for at generere FASTQ’er for Read1- og Read2-parret-end-læsninger samt GemCode-stribe-koden på 14 bp. UMI-tags med 10 basepar blev ekstraheret fra Read2 (14 biblioteker blev fremstillet med UMI-tags med 5 bp, som anført i supplerende tabel 1, på grund af en tidligere iteration af metoderne. For disse prøver blev UMI-tags på 5 bp ekstraheret fra Read2.). Derefter blev Read1, som indeholder cDNA-indsatsen, justeret til et passende referencegenom ved hjælp af STAR35. For museceller blev mm10 anvendt, og for humane celler blev hg19 anvendt. For prøver med blandinger af museceller og humane celler blev foreningen af hg19 og mm10 anvendt. For ERCC-prøver blev ERCC-referencen (https://tools.thermofisher.com/content/sfs/manuals/cms_095047.txt) anvendt.

Dernæst blev GemCode-stregkoder og UMI’er filtreret. Alle de kendte anførte stregkoder, der er 1-Hamming-afstand fra en observeret stregkode, blev taget i betragtning. Derefter beregnes den efterfølgende sandsynlighed for, at den observerede stregkode er fremkommet som følge af en sekventeringsfejl, på baggrund af den observerede stregkodes basekvaliteter og den forudgående sandsynlighed for at observere kandidatstregkoden (taget fra den samlede stregkodeantalfordeling). Hvis den efterfølgende sandsynlighed for en kandidatstregkode er mindst 0,975, korrigeres stregkoden til den kandidatstregkode, der har den højeste efterfølgende sandsynlighed. Hvis alle kandidatsekvenser er lige sandsynlige, vælges den, der optræder først i leksikalsk rækkefølge.

UMI’er med sekventeringskvalitetsscore >10 blev betragtet som gyldige, hvis de ikke var homopolymerer. Kval=10 indebærer en nøjagtighed på 90 % af basekaldet. En UMI, der er 1-Hamming-afstand fra en anden UMI (med flere læsninger) for den samme cellestregkode og det samme gen, korrigeres til UMI’en med flere læsninger. Denne fremgangsmåde er næsten identisk med den i Jaitin et al.4 og ligner den i Klein et al.8 (selv om Klein et al.8 også brugte UMI’er til at løse multimapped læsninger, hvilket ikke blev implementeret her).

Sidst blev PCR-duplikater markeret, hvis to sæt læsningspar delte en stregkodesekvens, et UMI-tag og et gen-id (Ensembl GTF’er GRCh37.82, ftp://ftp.ensembl.org/pub/grch37/release-84/gtf/homo_sapiens/Homo_sapiens.GRCh37.82.gtf.gz og GRCm38.84, ftp://ftp.ensembl.org/pub/release-84/gtf/mus_musculus/Mus_musculus.GRCm38.84.gtf.gz, blev brugt). Kun sikkert kortlagte (MAPQ=255), ikke-PCR-duplikater med gyldige stregkoder og UMI’er blev anvendt til at generere gen-stregkodematrix.

Stregkoder for celler blev bestemt på grundlag af fordelingen af UMI-tællinger. Alle de øverste stregkoder inden for samme størrelsesorden (>10% af den øverste n-te stregkode, hvor n er 1% af det forventede genvundne celletal) blev betragtet som celle-stregkoder. Antallet af læsninger, der giver meningsfulde oplysninger, beregnes som produktet af fire metrikker: (1) gyldige stregkoder; (2) gyldigt UMI; (3) forbundet med en celle-stregkode; og (4) sikkert kortlagt til exoner.

I blandingsforsøgene med mus og mennesker blev multipletraten defineret som to gange hastigheden af celle-stregkoder med signifikante UMI-tællinger fra både mus og mennesker, hvor den øverste 1% af UMI-tællingerne blev betragtet som signifikant. Omfanget af stregkode-crosstalk blev vurderet ved brøkdelen af muselæsninger i humane stregkoder eller omvendt.

Prover, der er behandlet fra flere kanaler, kan kombineres ved at sammenkæde gen-celle-stribekode-matricer. Denne funktionalitet findes i Cell Ranger R-sættet (http://support.10xgenomics.com/single-cell/software/pipelines/latest/rkit). Sekventeringsdata fra flere sekventeringskørsler af et bibliotek kan kombineres ved at tælle ikke-duplikerede læsninger. Denne funktionalitet findes i Cell Ranger-pipelinen. Desuden kan sekventeringsdata underprøves for at opnå et givet antal UMI-tællinger pr. celle. Denne funktionalitet leveres også i Cell Ranger R-kittet og er nyttig, når data fra flere prøver kombineres med henblik på sammenligning.

PCA-analyse af blanding af Jurkat- og 293T-celler

Gen-celle-barcode-matrixen fra hver af de fire prøver blev sammenkædet. Kun gener med mindst én UMI-tælling, der er påvist i mindst én celle, er anvendt. UMI-normalisering blev udført ved først at dividere UMI-tællingerne med de samlede UMI-tællinger i hver celle, efterfulgt af multiplikation med medianen af de samlede UMI-tællinger på tværs af cellerne. Derefter tog vi den naturlige logaritme af UMI-tallene. Endelig blev hvert gen normaliseret således, at det gennemsnitlige signal for hvert gen er 0, og standardafvigelsen er 1. PCA blev udført på den normaliserede matrix af gen-stribe-koder. De normaliserede UMI-tællinger for hvert gen bruges til at vise ekspression af en markør i et tSNE-plot.

SNV-analyse af Jurkat- og 293T scRNA-seq-data

SNV’er blev kaldt ved at køre Freebayes 1.0.2 (ref. 36) på genom BAM produceret af Cell Ranger. SNV’er af høj kvalitet (SNV calling Qual>=100 med mindst 10 UMI-tællinger fra mindst to celler; indels ignoreres), som kun blev observeret i Jurkat- eller 293T-celler (men ikke begge), blev udvalgt. Cellerne blev mærket som Jurkat eller 293T på grundlag af Jurkat- og 293T-specifikke SNV-tællinger, hvor brøkdelen af tællinger fra den anden art er <0,2. Celler med en fraktion af SNV fra en af arterne på mellem 0,2 og 0,8 betragtes som multiplets. Den udledte multipletrate er 2* den observerede multipletrate (for at tage højde for Jurkat:Jurkat- og 293T:293T-multiplets).

PCA- og tSNE-analyse af PBMC’er

Gener med mindst én UMI-tælling, der er påvist i mindst én celle, er anvendt. De øverste 1.000 mest variable gener blev identificeret på grundlag af deres gennemsnit og spredning (varians/middelværdi), hvilket svarer til den fremgangsmåde, der blev anvendt af Macoscko et al.7 Generne blev placeret i 20 bins baseret på deres gennemsnitlige udtryk. Normaliseret spredning er beregnet som den absolutte forskel mellem spredning og median spredning af ekspressionsgennemsnittet, normaliseret med median absolut afvigelse inden for hver bin.

PCA blev kørt på den normaliserede gen-barcode matrix af de øverste 1.000 mest variable gener for at reducere antallet af feature (gen) dimensioner. UMI-normalisering blev udført ved først at dividere UMI-tællinger med de samlede UMI-tællinger i hver celle, efterfulgt af multiplikation med medianen af de samlede UMI-tællinger på tværs af cellerne. Derefter tog vi den naturlige logaritme af UMI-tallene. Endelig blev hvert gen normaliseret således, at det gennemsnitlige signal for hvert gen er 0, og standardafvigelsen er 1. PCA blev udført på den normaliserede gen-barcode-matrix. Efter at have kørt PCA blev Barnes-hut37-approximation til t-SNE16 udført på de første 50 PC’er for at visualisere cellerne i et todimensionalt rum. Der blev anvendt 50 PC’er, fordi: (1) det ville tage meget lang tid at bruge alle PC’er med tSNE-analyse; (2) de forklarede ∼25% af den samlede varians. K-means15 clustering blev kørt for at gruppere celler til clusteringanalysen. k=10 blev valgt baseret på summen af kvadreret fejl scree plot (Supplerende fig. 5d).

Identificering af klyngespecifikke gener og markørbaseret klassifikation

For at identificere gener, der er beriget i en specifik klynge, blev den gennemsnitlige ekspression af hvert gen beregnet på tværs af alle celler i klyngen. Derefter blev hvert gen fra klyngen sammenlignet med medianekspressionen af det samme gen fra celler i alle andre klynger. Generne blev rangeret på baggrund af deres ekspressionsforskel, og de 10 mest berigede gener fra hver klynge blev udvalgt. Til hierarkisk clustering blev parvis korrelation mellem hver klynge beregnet, og centreret ekspression af hvert gen blev brugt til visualisering ved hjælp af heatmap.

Klassificering af PBMC’er blev udledt af annotationen af klyngespecifikke gener. I tilfælde af klynge 10 blev der påvist markørudtryk af flere celletyper (f.eks. B, dendritisk og T). Da den relative klyngestørrelse for B, dendritisk og T er henholdsvis 5,7 %, 6,6 % og 81 %, ville vi forvente, at klynge 10 (som kun er 0,5 %) ville indeholde multiplets, der hovedsagelig består af B:dendritisk (0.36 %) og B:dendritisk:T (0,3 %).

Selektion af rensede subpopulationer af PBMC’er

Hver population af rensede PBMC’er blev downsamplet til ∼16k læser pr. celle. PCA, tSNE og k-means clustering blev udført for hver downsamplet matrix efter de samme trin som beskrevet i PCA- og t-SNE-analyse af PBMC’er. Der blev kun påvist én klynge i de fleste prøver, hvilket er i overensstemmelse med FACS-analyserne (Supplerende figur 6). For prøver med mere end én klynge blev kun de klynger, der viste den forventede markørgenekspression, udvalgt til downstream-analyse. For CD14+ monocytter blev der observeret to klynger, som blev identificeret som CD14+ monocytter og dendritiske celler baseret på ekspression af henholdsvis markørgenerne FTL og CLEC9A.

Celleklassifikationsanalyse ved hjælp af rensede PBMC’er

Hver population af rensede PBMC’er blev downsamplet til ∼16k sikkert kortlagte læserater pr. celle. Derefter blev der beregnet en gennemsnitlig (gennemsnitlig) genekspressionsprofil på tværs af alle celler. Derefter blev genekspressionen fra hver celle i den komplekse population sammenlignet med genekspressionsprofilerne for rensede populationer af PBMC’er ved hjælp af Spearmans korrelation. Cellen blev tildelt ID’et for den rensede population, hvis den havde den højeste korrelation med den pågældende population. Bemærk, at forskellen mellem den højeste og den næsthøjeste korrelation var lille for nogle celler (f.eks. forskellen mellem cytotoksiske T- og NK-celler), hvilket tyder på, at celletildelingen ikke var så sikker for disse celler. Nogle få af de rensede PBMC-populationer overlappede med hinanden. F.eks. omfatter CD4+ T-hjælperceller alle CD4+ celler. Det betyder, at celler fra denne prøve vil overlappe med celler fra prøver, der indeholder CD4+ celler, herunder CD4+/CD25+ T reg, CD4+/CD45RO+ T memory, CD4+/CD45RA+/CD25- naive T. Når en celle blev tildelt ID’et CD4+ T-hjælpercelle på grundlag af korrelationsscoren, blev den næsthøjeste korrelation således kontrolleret for at se, om den var en af CD4+ prøverne. Hvis det var tilfældet, blev cellens ID opdateret til den celletype med den næsthøjeste korrelation. Den samme procedure blev udført for CD8+ cytotoksisk T og CD8+/CD45RA+ naiv cytotoksisk T (som er en undergruppe af CD8+ cytotoksisk T).

Den R-kode, der blev anvendt til at analysere 68k PBMC’er og rensede PBMC’er, kan findes her: https://github.com/10XGenomics/single-cell-3prime-paper.

Cell clustering og klassificering med Seurat

Gen-celle-barcode-matrixen af 68k PBMC’er blev log-transformeret som input til Seurat. De 469 mest variable gener, der blev udvalgt af Seurat, blev anvendt til at beregne PC’erne. De første 22 PC’er var signifikante (P<0,01) baseret på den indbyggede jackstraw-analyse og blev brugt til tSNE-visualisering. Celleklassificering blev taget fra celleklassificeringsanalyse ved hjælp af rensede PBMC’er.

Sammenligning mellem friske og frosne PBMC’er

Sekventeringsdataene fra 68k friske PBMC’er og 3k frosne PBMC’er blev down-samplet således, at hver prøve har ∼14k sikkert kortlagte læsninger pr. celle. Kun gener, der er detekteret i mindst én celle, blev inkluderet i sammenligningen, som anvender gennemsnittet af hvert gen på tværs af alle celler.

Til sammenligning af celleklassificering mellem rensede og frosne PBMC’er samlede vi alle celler, der var mærket som T- eller naturlige dræberceller, i en pulje. Dette skyldes, at subpopulationerne inden for T og mellem T- og naturlige dræberceller nogle gange er vanskelige at gruppere separat. Vi ønskede ikke, at sammenligningen mellem friske og frosne celler skulle påvirkes af de anvendte clusteringmetoder.

SNV-baseret genotypetildeling

SNV’er blev kaldt ved at køre Freebayes 1.0.2 (ref. 36) på genom BAM produceret af Cell Ranger. Kun SNV’er med støtte fra mindst to cellestregkoder, med en minimal SNV Qual score >=30, minimal SNV base Qual>=1 blev medtaget. Antallet af referencealleler (R) og alternative alleler (A) blev beregnet for hver SNV, hvilket resulterede i en matrix af UMI-tal for celle-referencealleler og UMI-tal for celle-alternativalleler. Disse matricer blev modelleret som en blanding af to genomer, hvor sandsynligheden for en af de tre genotyper (R/R, R/A eller A/A) på et sted blev antaget at være binomialt fordelt med en fast fejlprocent på 0,1 %. For hver prøve blev der parallelt udledt to modeller, en hvor der kun er ét genom til stede (K=1) og en anden hvor der er to genomer til stede (K=2). Inferens af modelparametrene (celle-til-genom-tilknytninger og K-sættene af genotyper) blev udført ved hjælp af en Gibbs sampler til at tilnærme sig deres posteriorfordelinger. For at afhjælpe label-switching-problemet i Monte Carlo-inferens af blandingsmodeller blev der foretaget ommærkning af de udtagne celle-til-genom-tildelinger i overensstemmelse med Stephens et al.38

I in silico-celleblandingseksperimenter, hvor K=2-modellen ikke kunne adskille de to genomer tilstrækkeligt, rapporterede den en fordeling af posterior-sandsynligheder nær 0,5 for celle-genom-tildelingerne, hvilket indikerer manglende tillid til disse tildelinger. Vi anvendte et krav om, at 90 % af cellerne skal have en posterior sandsynlighed >75 % for at vælge K=2-modellen frem for K=1-modellen. Valg af K=1 indikerer, at blandingsfraktionen er under metodens detektionsniveau, som i in silico-blandingseksperimenter blev bestemt til at være 4 % af 6.000 celler.

Genotypesammenligning med den rene prøve

For at sikre tilknytningen af genotyper til individer blev der kun taget hensyn til fælles SNV’er mellem genotypegruppen og den rene prøve. Derefter blev den gennemsnitlige genotype for alle cellerne sammenlignet med genotypen for den rene prøve. For at opnå en vis basislinje for den procentvise genotypeoverlapning mellem forskellige individer foretog vi parvis sammenligning af genotyper, der blev kaldt fra de samme individer (11 parvise sammenligninger) eller fra forskellige individer (15 parvise sammenligninger). Det procentvise genotypeoverlap mellem de samme individer er i gennemsnit ∼98±0,3 %, mens det procentvise genotypeoverlap mellem de forskellige individer er i gennemsnit ∼73±2 %.

PCA- og tSNE-analyse af BMMC’er

Data fra seks prøver blev anvendt: to sunde kontroller, AML027 før og efter transplantation samt AML035 før og efter transplantation. Hver prøve blev downsamplet til ∼10k sikkert kortlagte læsninger pr. celle. Derefter blev gen-celle stregkodematrixen fra hver prøve sammenkædet. PCA, tSNE og k-means clustering blev udført på den sammenlagte matrix efter de samme trin som beskrevet i PCA- og tSNE-analysen af PBMC’er. Til k-means-gruppering blev K=10 anvendt på grundlag af bøjningen i summen af kvadreret fejl scree-plot.

Klyngespecifikke gener blev identificeret efter de trin, der er beskrevet i “Identifikation af klyngespecifikke gener og markørbaseret klassificering”. Klassificering blev tildelt på grundlag af klyngespecifikke gener og på grundlag af ekspression af nogle velkendte markører for immuncelletyper. ‘Blasts and Immature Ery 1’ henviser til klynge 4, som udtrykker CD34, en markør for hæmatopoietiske progenitorer39 , og Gata2, en markør for tidlige erytroider40. ‘Immature Ery 2’ henviser til klynge 5 og 8, som viser udtryk for Gata1, en transkriptionsfaktor, der er afgørende for erythropoiese41 , men ikke CD71, som ofte findes i mere engagerede erytroide celler39. “Immature Ery 3” henviser til klynge 1, som viser udtryk for CD71. “Mature Ery” henviser til klynge 2. HBA1, en markør for modne erytroide celler, er fortrinsvis påvist i klynge 2. Klynge 3 blev tildelt som “umodne granulocytter” på grund af ekspressionen af tidlige granulocytmarkører såsom AZU1 og IL8 (ref. 42) og manglen på ekspression af CD16. Klynge 7 blev tildelt som “Monocytter” på grund af ekspressionen af bl.a. CD14 og FCN1. ‘B’ henviser til klynge 6 og 9 på grund af markører som CD19 og CD79A. ‘T’ henviser til klynge 10 på grund af markører som CD3D og CD8A.

Databesiddelse

Alle relevante data er tilgængelige hos forfatterne. Single-cell RNA-seq-data er blevet deponeret i Short Read Archive under tiltrædelsesnummer SRP073767. Data er også tilgængelige på http://support.10xgenomics.com/single-cell/datasets. Analysekoden for 68k PBMC-analysen er tilgængelig på https://github.com/10XGenomics/single-cell-3prime-paper.