Profilage transcriptionnel numérique massivement parallèle de cellules uniques

Imagerie à grande vitesse des billes de gel et des cellules dans les GEM

Un microscope (Nikon Ti-E, objectif × 10) et une caméra vidéo à grande vitesse (Photron SA5, fréquence d’images=4 000 s-1) ont été utilisés pour imager chaque GEM au fur et à mesure qu’ils étaient générés dans la puce microfluidique. Un logiciel d’analyse personnalisé a été utilisé pour compter le nombre de GEM générés et le nombre de billes présentes dans chaque GEM, en se basant sur la détection des bords et le contraste entre les bords des billes et des GEM et le liquide adjacent. Les résultats de l’analyse sont résumés dans la figure 1c. Pour estimer la distribution des cellules dans les GEM, un comptage manuel a été utilisé pour ∼28k images d’une vidéo sur un sous-ensemble de GEMs. Les résultats indiquent une adhésion approximative à une distribution de Poisson. Cependant, le pourcentage d’encapsulations cellulaires multiples était 16% plus élevé que la valeur attendue, peut-être en raison d’une erreur de sous-échantillonnage ou d’interactions cellule-cellule (certains amas de deux cellules ont été observés lors du comptage manuel) (Supplementary Fig. 1b).

Lignées cellulaires et échantillons de patients transplantés

Les cellules Jurkat (ATCC TIB-152), 293T (ATCC CRL-11268) et 3T3 (ATCC CRL-1658) ont été acquises auprès de l’ATCC et cultivées conformément aux directives de l’ATCC. Les PBMCs fraîches, les PBMCs congelées et les BMMCs ont été achetées chez ALLCELLS. Les PBMC congelées du donneur A ont été fabriquées à partir des PBMC fraîches du donneur A en mélangeant délicatement 1e6 cellules dans un milieu de congélation (15 % de diméthylsulfoxyde (DMSO) dans le milieu Dulbecco modifié d’Iscove contenant 20 % de FBS), puis en les refroidissant dans le CoolCell FTS30 (BioCision) à -80 °C pendant au moins 4 h avant de les transférer dans l’azote liquide pour les conserver pendant 3 semaines.

Le comité d’examen institutionnel du Fred Hutchinson Cancer Research Center a approuvé l’étude sur les échantillons de transplantation. Les procédures suivies étaient conformes à la déclaration d’Helsinki de 1975 et à la règle commune. Les échantillons ont été obtenus après que les patients aient donné leur consentement éclairé écrit sur les analyses moléculaires. Nous avons identifié les patients atteints de LAM ayant subi une greffe allogénique de cellules souches hématopoïétiques au Fred Hutchinson Cancer Research Center. Le diagnostic de LAM a été établi selon les critères révisés de l’Organisation mondiale de la santé33.

Des aspirats de moelle osseuse ont été obtenus pour les tests cliniques standard 20-30 jours avant la transplantation et en série après la transplantation selon le protocole de traitement. Les aliquotes d’aspiration de moelle osseuse ont été traitées dans les 2 h suivant le prélèvement. Les CMMB ont été isolées par centrifugation à travers un gradient Ficoll (Histopaque-1077 ; Sigma Life Science, St Louis, MO, USA). Les CMMB ont été recueillies à l’interface sérum-Ficoll avec une pipette Pasteur jetable et transférées dans un tube conique de 50 ml contenant 2 % de sérum de patient dans 1 × PBS. Les CMMB ont été comptées à l’aide d’un hémacytomètre et leur viabilité a été évaluée à l’aide de bleu Trypan. Les BMMC ont été remises en suspension dans un milieu de congélation à 90 % de FBS, 10 % de DMSO et congelées à l’aide d’un Nalgene Mr Frosty (Thermo Scientific) dans un congélateur à -80 °C pendant 24 h avant d’être transférées dans l’azote liquide pour un stockage à long terme.

Estimation du contenu en ARN par cellule

La quantité d’ARN par type de cellule a été déterminée en quantifiant (Qubit ; Invitrogen) l’ARN extrait (Maxwell RSC simplyRNA Cells Kit) de plusieurs nombres différents et connus de cellules.

Préparation des cellules

Les cellules fraîches ont été récoltées, lavées avec 1 × PBS et remises en suspension à 1 × 106 cellules par ml dans 1 × PBS et 0,04% d’albumine de sérum bovin. Les PBMC fraîches ont été congelées à 10 × en remettant en suspension les PBMC dans du DMEM+40% FBS+10% DMSO, en les congelant à -by °C dans un CoolCell® FTS30 (BioCision) puis en les plaçant dans de l’azote liquide pour les conserver.

Les flacons de cellules congelées provenant de l’étude ALLCELLS et des études de transplantation ont été rapidement décongelés dans un bain-marie à 37 °C pendant ∼2 min. Les flacons ont été retirés lorsqu’un minuscule cristal de glace était laissé. Les PBMC décongelées ont été lavées deux fois dans le milieu, puis remises en suspension dans 1 × PBS et 0,04 % d’albumine de sérum bovin à température ambiante. Les cellules ont été centrifugées à 300 r.c.f. pendant 5 minutes à chaque fois. Les BMMC décongelées ont été lavées et remises en suspension dans 1 × PBS et 20 % de FBS. La concentration finale des cellules décongelées était de 1 × 106 cellules par ml.

Construction de bibliothèques de séquençage à l’aide de la plateforme GemCode

Les suspensions cellulaires ont été chargées sur un instrument GemCode Single-Cell (10x Genomics, Pleasanton, CA, USA) pour générer des GEM unicellulaires. Les librairies RNA-Seq unicellulaires ont été préparées à l’aide du kit GemCode Single-Cell 3′ Gel Bead and Library (désormais vendu sous les références 120230, 120231, 120232, 10x Genomics). La RT des GEM a été effectuée dans un thermocycleur C1000 Touch avec module de réaction 96 puits profonds (Bio-Rad ; P/N 1851197) : 55 °C pendant 2 h, 85 °C pendant 5 min ; maintenu à 4 °C. Après la réaction rapide, les GEM ont été brisés et l’ADNc simple brin a été nettoyé avec les billes de silane DynaBeads MyOne (Thermo Fisher Scientific ; P/N 37002D) et le kit de réactifs SPRIselect (0,6 × SPRI ; Beckman Coulter ; P/N B23318). L’ADNc a été amplifié à l’aide du thermocycleur C1000 Touch avec module de réaction 96 puits profonds : 98 °C pendant 3 min ; cyclage 14 × : 98 °C pendant 15 s, 67 °C pendant 20 s et 72 °C pendant 1 min ; 72 °C pendant 1 min ; maintien à 4 °C. Le produit ADNc amplifié a été nettoyé avec le kit de réactifs SPRIselect (0,6 × SPRI). L’ADNc a ensuite été cisaillé à ∼200 pb à l’aide d’un système Covaris M220 (Covaris ; P/N 500295). Des bibliothèques de séquençage indexées ont été construites à l’aide des réactifs du kit de bibliothèque 3′ GemCode Single-Cell, en suivant les étapes suivantes : (1) réparation des extrémités et A-tail ; (2) ligature des adaptateurs ; (3) nettoyage postligature avec SPRIselect ; (4) PCR et nettoyage de l’indexation des échantillons. Les bibliothèques de séquençage de code-barres ont été quantifiées par PCR quantitative (KAPA Biosystems Library Quantification Kit for Illumina platforms P/N KK4824). Les bibliothèques de séquençage ont été chargées à 2,1 pM sur un NextSeq500 d’Illumina avec des kits 2 × 75 paired-end utilisant la longueur de lecture suivante : 98 pb Read1, 14 pb I7 Index, 8 pb I5 Index et 10 pb Read2. Certaines bibliothèques antérieures ont été réalisées avec des UMI de 5 nt, et on a obtenu une Read2 de 5 pb à la place. Ces bibliothèques ont été documentées dans le tableau supplémentaire 1.

ErCC assay

ERCC synthetic spike-in RNAs (Thermo Fisher Scientific ; P/N 4456740) ont été dilués (1:10 ou 1:50) et chargés dans un instrument GemCode Single-Cell, remplaçant les cellules normalement utilisées pour générer des GEMs. Les mélanges Spike-in Mix1 et Mix2 ont tous deux été testés. Un protocole légèrement modifié a été utilisé car seule une petite fraction des GEMs a été collectée pour la RT et l’amplification de l’ADNc. Une fois la RT des GEM terminée, 1,25 μl de l’émulsion a été retiré et ajouté à un mélange biphasique d’agent de récupération (125 μl) (P/N 220016) et d’additif 1 25 mM (30 μl) (P/N 220074, 10x Genomics). L’agent de récupération a ensuite été retiré et la solution aqueuse restante a été nettoyée avec le kit de réactifs SPRISelect (0,8 × SPRI). L’ADNc a été amplifié à l’aide du thermocycleur C1000 Touch avec module de réaction 96 puits profonds : 98 °C pendant 3 min ; cyclage 14 × : 98 °C pendant 15 s, 67 °C pendant 20 s et 72 °C pendant 1 min ; 72 °C pendant 1 min ; maintien à 4 °C. Le produit d’ADNc amplifié a été nettoyé avec le kit de réactifs SPRIselect (0,8 × ). L’ADNc a ensuite été cisaillé à ∼200 pb à l’aide d’un système Covaris M220 pour construire des bibliothèques indexées par échantillon avec des adaptateurs Genomics 10x. Le nombre attendu de molécules ERCC a été calculé en fonction de la quantité de molécules ERCC utilisées et des facteurs de dilution de l’échantillon. Les comptages ont été comparés aux comptages de molécules détectées (comptages UMI) pour calculer l’efficacité de conversion.

test ddPCR

Les cellules Turkat ont été utilisées dans des tests ddPCR pour estimer l’efficacité de conversion comme suit : (1) la quantité d’ARN par cellule Jurkat a été déterminée en quantifiant (Qubit, Invitrogen) l’ARN extrait (kits de purification d’ARN Maxwell) de plusieurs cellules Jurkat en nombre connu. (2) Une RT-ddPCR en masse (Bio-Rad One-Step RT-ddPCR Advanced Kit for Probes 1864021) a été réalisée sur l’ARN extrait pour déterminer le nombre de copies par cellule de huit gènes sélectionnés. (3) Environ 5 000 cellules Jurkat ont été traitées à l’aide de la plateforme GemCode Single-Cell 3′, et de l’ADNc simple brin a été recueilli après RT dans des GEM selon les protocoles indiqués dans la section  » Construction de librairies de séquençage à l’aide de la plateforme GemCode « . Les copies d’ADNc des huit gènes ont été déterminées par ddPCR (Bio-Rad ddPCR Supermix for Probes (no dUTP) P/N 1863024). Le nombre réel de cellules Jurkat a été déterminé en séquençant un sous-ensemble de réactions GEM-RT sur un MiSeq. L’efficacité de conversion est le rapport entre les copies d’ADNc par cellule (étape 3) et les copies d’ARN par cellule provenant de la RT-ddPCR en vrac (étape 2), en supposant une efficacité de 50 % dans la RT-ddPCR34.

Les séquences de sonde pour le test ddPCR sont les suivantes : SERAC1_f, 5′-CACGAGCCGCCAGC-3′ et SERAC1_r, 5′-TCTGCAACAGATGACGCAATAAG-3′ ; SERAC1_p : /56-FAM/CGCCTGCCG/ZEN/GCAGAATGTC/3IABkFQ/. AP1S3_f, 5′-GAAGCAGCCATGGTCTAAGC-3′ et AP1S3_r, 5′-CCTTGTCGACTGAAGAGCAATATG-3′ ; AP1S3_p : /56-FAM/CGGCCCAGC/ZEN/CACGATGATACAT/3IABkFQ/OR. AOV1_f, 5′-CCGGAAGTGGGTCTCGTOR-3′ et AOV1_r, 5′-TTCTTCATAGCCTTCCCGATACCOR-3′ ; AOV1_p : /56-FAM/TCGTGATGG/ZEN/CGGATGAGAGGTTTCA/3IABkFQ/. DOLPP1_f, 5′-ATGGCAGCGGACGGA-3′ et DOLPP1_r, 5′-GGCTCAGGTAGGCAAGGA-3′ ; DOLPP1_p : /56-FAM/CCACGTCGA/ZEN/ATATCCTGCAGGTGATCT/3IABkFQ/. KPNA6_f, 5′-TGAAAGCTGCCGCTGAAG-3′ et KPNA6_r, 5′-CCCTGGGCTCGCCAT-3′ ; KPNA6_p : /56-FAM/CGGACCCGC/ZEN/GATGGAGACC/3IABkFQ/. ITSN2_f, 5′-GTGACAGGCTACGCAACAG-3′ et ITSN2_r, 5′-TCCTGAGTTTTCCTTGCTAGCT-3′ ; ITSN2_p : /56-FAM/AGGGCGCCA/ZEN/GATGGCTGA/3IABkFQ/. LCMT1_f, 5′-GTCGACCCCGCTTCCA-3′ et LCMT1_r, 5′-GGTCATGCCAGTAGCCAATG-3′ ; LCMT1_p : /56-FAM/ATGCTTCCC/ZEN/TGTGCAAGGTTTGC/3IABkFQ/. AP2M1_f, 5′-GCAGCGGGCAGACG-3′ et AP2M1_r, 5′-ATGGCGGCAGATCAGTCT-3′ ; AP2M1_p : /56-FAM/CATCGCTCT/ZEN/GAGAACAGACCTGGTG/3IABkFQ/.

Calcul de l’efficacité de capture des cellules

L’efficacité est calculée en prenant le rapport du nombre de cellules détectées par séquençage par rapport au nombre de cellules chargées dans la puce. Ce dernier est déterminé à partir de (volume ajouté × concentration d’entrée des cellules). La concentration d’entrée des cellules a été déterminée à l’aide d’un compteur cellulaire automatisé Countess II (Thermo Fisher Scientific). Il convient de noter qu’il existe une erreur de 15 à 20 % dans le comptage des cellules, ce qui pourrait expliquer au moins une partie de la variabilité des efficacités calculées.

Dosage du chimérisme

Le système PowerPlex 16 (Promega) a été utilisé conjointement avec un analyseur génétique Applied Biosystems (Life Technologies) 3130xl. Les BMMC de donneurs ont été utilisées comme ligne de base de référence.

Alignement, attribution de codes-barres et comptage des UMI

La suite logicielle monocellulaire Cell Ranger a été utilisée pour effectuer le démultiplexage des échantillons, le traitement des codes-barres et le comptage des gènes 3′ monocellulaires (http://software.10xgenomics.com/single-cell/overview/welcome). Tout d’abord, le démultiplexage de l’échantillon a été effectué sur la base de la lecture de l’index de l’échantillon de 8 pb afin de générer des FASTQ pour les lectures paires-end Read1 et Read2, ainsi que le code-barres GemCode de 14 pb. Des étiquettes UMI de 10 paires de bases ont été extraites de Read2 (14 bibliothèques ont été réalisées avec des étiquettes UMI de 5 paires de bases, comme indiqué dans le tableau supplémentaire 1, en raison d’une itération antérieure des méthodes. Pour ces échantillons, les balises UMI de 5 pb ont été extraites de Read2). Ensuite, Read1, qui contient l’insert d’ADNc, a été aligné sur un génome de référence approprié à l’aide de STAR35. Pour les cellules de souris, on a utilisé mm10 et pour les cellules humaines, on a utilisé hg19. Pour les échantillons contenant des mélanges de cellules de souris et humaines, l’union de hg19 et mm10 a été utilisée. Pour les échantillons ERCC, la référence ERCC (https://tools.thermofisher.com/content/sfs/manuals/cms_095047.txt) a été utilisée.

Puis, les codes-barres GemCode et les UMI ont été filtrés. Toutes les listes connues de codes à barres qui sont à une distance de 1-Hamming d’un code à barres observé sont considérées. Ensuite, la probabilité postérieure que le code-barres observé ait été produit par une erreur de séquençage est calculée, étant donné les qualités de base du code-barres observé et la probabilité antérieure d’observer le code-barres candidat (tirée de la distribution globale du nombre de codes-barres). Si la probabilité postérieure de tout code-barres candidat est d’au moins 0,975, alors le code-barres est corrigé en fonction du code-barres candidat ayant la probabilité postérieure la plus élevée. Si toutes les séquences candidates sont également probables, alors celle qui apparaît en premier par ordre lexical est choisie.

Les UMI avec un score de qualité de séquençage >10 ont été considérées comme valides si elles n’étaient pas des homopolymères. Qual=10 implique une précision d’appel de base de 90%. Un UMI qui est à une distance de 1-Hamming d’un autre UMI (avec plus de lectures) pour le même code-barres cellulaire et le même gène est corrigé vers l’UMI avec plus de lectures. Cette approche est presque identique à celle de Jaitin et al.4, et similaire à celle de Klein et al.8 (bien que Klein et al.8 ait également utilisé les UMI pour résoudre les lectures multimappées, ce qui n’a pas été mis en œuvre ici).

Enfin, les doublons PCR ont été marqués si deux ensembles de paires de lectures partageaient une séquence de code-barres, une étiquette UMI et un ID de gène (les GTF Ensembl GRCh37.82, ftp://ftp.ensembl.org/pub/grch37/release-84/gtf/homo_sapiens/Homo_sapiens.GRCh37.82.gtf.gz et GRCm38.84, ftp://ftp.ensembl.org/pub/release-84/gtf/mus_musculus/Mus_musculus.GRCm38.84.gtf.gz, ont été utilisés). Seuls les duplicatas non PCR, cartographiés avec confiance (MAPQ=255), avec des codes-barres et des UMI valides, ont été utilisés pour générer la matrice gène-code-barres.

Les codes-barres des cellules ont été déterminés sur la base de la distribution du nombre d’UMI. Tous les codes à barres supérieurs dans le même ordre de grandeur (>10% du nième code à barres supérieur, où n est 1% du nombre de cellules récupérées prévu) ont été considérés comme des codes à barres de cellules. Le nombre de lectures qui fournissent des informations significatives est calculé comme le produit de quatre métriques : (1) codes-barres valides ; (2) UMI valides ; (3) associés à un code-barres cellulaire ; et (4) cartographiés avec confiance à des exons.

Dans les expériences de mélange de souris et d’humains, le taux de multiplets a été défini comme étant le double du taux de codes-barres cellulaires avec des comptages UMI significatifs à la fois chez la souris et chez l’humain, où le top 1 % des comptages UMI était considéré comme significatif. L’étendue de la diaphonie des codes-barres a été évaluée par la fraction de lectures de souris dans les codes-barres humains, ou vice versa.

Les échantillons traités à partir de plusieurs canaux peuvent être combinés en concaténant les matrices gène-cellule-code-barres. Cette fonctionnalité est fournie dans le kit Cell Ranger R (http://support.10xgenomics.com/single-cell/software/pipelines/latest/rkit). Les données de séquençage provenant de plusieurs cycles de séquençage d’une bibliothèque peuvent être combinées en comptant les lectures non dupliquées. Cette fonctionnalité est fournie dans le pipeline Cell Ranger. En outre, les données de séquençage peuvent être sous-échantillonnées pour obtenir un nombre donné de comptages UMI par cellule. Cette fonctionnalité est également fournie dans le kit Cell Ranger R, et est utile pour combiner les données de plusieurs échantillons à des fins de comparaison.

Analyse PCA du mélange de cellules Jurkat et 293T

La matrice des codes-barres des gènes de chacun des quatre échantillons a été concaténée. Seuls les gènes avec au moins un comptage UMI détecté dans au moins une cellule sont utilisés. La normalisation des UMI a été effectuée en divisant d’abord les comptes UMI par les comptes UMI totaux dans chaque cellule, puis en les multipliant par la médiane des comptes UMI totaux dans toutes les cellules. Ensuite, nous avons pris le log naturel des comptes UMI. Enfin, chaque gène a été normalisé de telle sorte que le signal moyen de chaque gène soit égal à 0 et que l’écart-type soit égal à 1. L’ACP a été exécutée sur la matrice gène-code-barres normalisée. Les comptes UMI normalisés de chaque gène sont utilisés pour montrer l’expression d’un marqueur dans un graphique tSNE.

Analyse SNV des données scRNA-seq Jurkat et 293T

Les SNV ont été appelés en exécutant Freebayes 1.0.2 (réf. 36) sur le génome BAM produit par Cell Ranger. Les SNV de haute qualité (SNV calling Qual>=100 avec au moins 10 comptes d’UMI provenant d’au moins deux cellules ; les indels sont ignorés) qui ont été observés uniquement dans les cellules Jurkat ou 293T (mais pas les deux) ont été sélectionnés. Les cellules ont été étiquetées comme étant Jurkat ou 293T sur la base des comptages de SNV spécifiques à Jurkat et 293T, où la fraction des comptages provenant de l’autre espèce est <0,2. Les cellules dont la fraction de SNV de l’une ou l’autre espèce est comprise entre 0,2 et 0,8 sont considérées comme des multiplets. Le taux de multiplets inférés est 2* le taux de multiplets observés (pour tenir compte des multiplets Jurkat:Jurkat et 293T:293T).

Analyse PCA et tSNE des PBMC

Les gènes avec au moins un comptage UMI détecté dans au moins une cellule sont utilisés. Les 1 000 gènes les plus variables ont été identifiés en fonction de leur moyenne et de leur dispersion (variance/moyenne), ce qui est similaire à l’approche utilisée par Macoscko et al.7 Les gènes ont été placés dans 20 bacs en fonction de leur expression moyenne. La dispersion normalisée est calculée comme la différence absolue entre la dispersion et la dispersion médiane de la moyenne d’expression, normalisée par l’écart absolu médian au sein de chaque bin.

L’ACP a été exécutée sur la matrice gène-code à barres normalisée des 1 000 gènes les plus variables afin de réduire le nombre de dimensions des caractéristiques (gènes). La normalisation de l’UMI a été effectuée en divisant d’abord le nombre d’UMI par le nombre total d’UMI dans chaque cellule, puis en le multipliant par la médiane du nombre total d’UMI dans toutes les cellules. Ensuite, nous avons pris le log naturel des comptes d’UMI. Enfin, chaque gène a été normalisé de telle sorte que le signal moyen pour chaque gène soit égal à 0 et que l’écart-type soit égal à 1. L’ACP a été exécutée sur la matrice normalisée gène-code-barres. Après l’exécution de l’ACP, l’approximation de Barnes-hut37 à t-SNE16 a été effectuée sur les 50 premiers PC pour visualiser les cellules dans un espace bidimensionnel. Cinquante PC ont été utilisés car : (1) l’utilisation de tous les PC prendrait beaucoup de temps avec l’analyse tSNE ; (2) ils expliquaient ∼25 % de la variance totale. Le clustering K-means15 a été exécuté pour regrouper les cellules pour l’analyse de clustering. k=10 a été sélectionné sur la base de la somme des erreurs au carré scree plot (Supplementary Fig. 5d).

Identification des gènes spécifiques au cluster et classification basée sur les marqueurs

Pour identifier les gènes qui sont enrichis dans un cluster spécifique, l’expression moyenne de chaque gène a été calculée sur toutes les cellules du cluster. Ensuite, chaque gène du cluster a été comparé à l’expression médiane du même gène à partir des cellules de tous les autres clusters. Les gènes ont été classés en fonction de leur différence d’expression, et les 10 gènes les plus enrichis de chaque groupe ont été sélectionnés. Pour le clustering hiérarchique, la corrélation par paire entre chaque cluster a été calculée, et l’expression centrée de chaque gène a été utilisée pour la visualisation par carte thermique.

La classification des PBMC a été déduite de l’annotation des gènes spécifiques aux clusters. Dans le cas du cluster 10, l’expression des marqueurs de plusieurs types de cellules (par exemple, B, dendritiques et T) a été détectée. Puisque la taille relative des clusters B, dendritique et T est de 5,7 %, 6,6 % et 81 %, respectivement, nous nous attendons à ce que le cluster 10 (qui n’est que de 0,5 %) contienne des multiplets composés principalement de B:dendritique (0.36%) et B:dendritique:T (0,3%).

Sélection de sous-populations purifiées de PBMCs

Chaque population de PBMCs purifiées a été sous-échantillonnée à ∼16k reads par cellule. L’ACP, le tSNE et le clustering k-means ont été réalisés pour chaque matrice sous-échantillonnée, en suivant les mêmes étapes décrites dans l’analyse ACP et t-SNE des PBMCs. Un seul cluster a été détecté dans la plupart des échantillons, ce qui correspond aux analyses FACS (figure supplémentaire 6). Pour les échantillons présentant plus d’un cluster, seuls les clusters affichant l’expression génique marqueur attendue ont été sélectionnés pour une analyse en aval. Pour les monocytes CD14+, deux clusters ont été observés et identifiés comme monocytes CD14+ et cellules dendritiques sur la base de l’expression des gènes marqueurs FTL et CLEC9A, respectivement.

Analyse de classification des cellules à l’aide de PBMC purifiées

Chaque population de PBMC purifiées a été sous-échantillonnée à ∼16k lectures cartographiées avec confiance par cellule. Ensuite, un profil d’expression génique moyen (moyenne) de toutes les cellules a été calculé. Ensuite, l’expression génique de chaque cellule de la population complexe a été comparée aux profils d’expression génique des populations purifiées de PBMC par corrélation de Spearman. La cellule s’est vue attribuer l’ID de la population purifiée si elle présentait la corrélation la plus élevée avec cette population. Il est à noter que la différence entre la corrélation la plus élevée et la deuxième plus élevée était faible pour certaines cellules (par exemple, la différence entre les cellules T cytotoxiques et les cellules NK), ce qui suggère que l’attribution de la cellule n’était pas aussi fiable pour ces cellules. Quelques-unes des populations de PBMC purifiées se chevauchaient. Par exemple, les cellules T-helper CD4+ comprennent toutes les cellules CD4+. Cela signifie que les cellules de cet échantillon se chevaucheront avec les cellules des échantillons qui contiennent des cellules CD4+, y compris les T reg CD4+/CD25+, les T mémoire CD4+/CD45RO+, les T naïfs CD4+/CD45RA+/CD25-. Ainsi, lorsqu’une cellule se voyait attribuer l’ID de cellule T auxiliaire CD4+ sur la base du score de corrélation, la corrélation la plus élevée suivante était vérifiée pour voir s’il s’agissait d’un des échantillons CD4+. Si c’était le cas, l’ID de la cellule était mis à jour avec le type de cellule présentant la corrélation la plus élevée suivante. La même procédure a été effectuée pour les T cytotoxiques CD8+ et les T cytotoxiques naïfs CD8+/CD45RA+ (qui est un sous-ensemble des T cytotoxiques CD8+).

Le code R utilisé pour analyser les 68k PBMCs et les PBMCs purifiés peut être trouvé ici : https://github.com/10XGenomics/single-cell-3prime-paper.

Classement de cellules et classification avec Seurat

La matrice gène-cellule-code-barres de 68k PBMCs a été transformée en logarithme comme entrée pour Seurat. Les 469 gènes les plus variables sélectionnés par Seurat ont été utilisés pour calculer les PCs. Les 22 premiers PC étaient significatifs (P<0,01) selon l’analyse jackstraw intégrée, et ont été utilisés pour la visualisation tSNE. La classification des cellules a été reprise de l’analyse de classification des cellules utilisant des PBMC purifiées.

Comparaison entre les PBMC fraîches et congelées

Les données de séquençage de 68k PBMC fraîches et de 3k PBMC congelées ont été sous-échantillonnées de telle sorte que chaque échantillon comporte ∼14k lectures cartographiées avec confiance par cellule. Seuls les gènes détectés dans au moins une cellule ont été inclus pour la comparaison, qui utilise la moyenne de chaque gène sur toutes les cellules.

Pour la comparaison de la classification cellulaire entre les PBMC purifiées et congelées, nous avons regroupé toutes les cellules étiquetées comme T ou cellules tueuses naturelles. En effet, les sous-populations au sein de T et entre T et cellules tueuses naturelles sont parfois difficiles à regrouper séparément. Nous ne voulions pas que la comparaison entre les cellules fraîches et congelées soit affectée par les méthodes de regroupement utilisées.

Assignation de génotype basée sur les SNV

Les SNV ont été appelés en exécutant Freebayes 1.0.2 (réf. 36) sur le génome BAM produit par Cell Ranger. Seuls les SNVs supportés par au moins deux codes barres cellulaires, avec un score SNV Qual minimal >=30, SNV base Qual>=1 ont été inclus. Les nombres d’allèles de référence (R) et alternatifs (A) ont été calculés pour chaque SNV, produisant une matrice de nombres d’UMI de référence de cellule et de nombres d’UMI d’allèles alternatifs de cellule. Ces matrices ont été modélisées comme un mélange de deux génomes où la probabilité de n’importe lequel des trois génotypes (R/R, R/A ou A/A) à un site a été considérée comme étant distribuée de façon binomiale avec un taux d’erreur fixe de 0,1 %. Pour chaque échantillon, deux modèles ont été inférés en parallèle, un où un seul génome est présent (K=1) et un autre où deux génomes sont présents (K=2). L’inférence des paramètres des modèles (affectations cellule-génome et les K ensembles de génotypes) a été réalisée en utilisant un échantillonneur de Gibbs pour approximer leurs distributions postérieures. Pour améliorer le problème de changement d’étiquette dans l’inférence de Monte Carlo des modèles de mélange, le réétiquetage des assignations cellule-génome échantillonnées a été effectué selon Stephens et al.38

Dans les expériences de mélange de cellules in silico, lorsque le modèle K=2 n’a pas réussi à séparer les deux génomes de manière adéquate, il a signalé une distribution de probabilités postérieures proche de 0,5 pour les appels cellule-génome, indiquant un manque de confiance dans ces appels. Nous avons appliqué une exigence selon laquelle 90% des cellules ont une probabilité postérieure >75% pour sélectionner le modèle K=2 plutôt que le modèle K=1. La sélection de K=1 indique que la fraction de mélange est inférieure au niveau de détection de la méthode, qui, dans les expériences de mélange in silico, a été déterminée comme étant de 4 % de 6 000 cellules.

Comparaison des génotypes avec l’échantillon pur

Pour vérifier l’attribution des génotypes aux individus, seuls les SNV partagés entre le groupe génotype et l’échantillon pur ont été considérés. Ensuite, le génotype moyen de toutes les cellules a été comparé à celui de l’échantillon pur. Pour obtenir une base de référence pour le % de chevauchement des génotypes entre différents individus, nous avons effectué une comparaison par paire des génotypes appelés des mêmes individus (11 comparaisons par paire) ou de différents individus (15 comparaisons par paire). Le pourcentage de chevauchement génotypique entre les mêmes individus est en moyenne de ∼98±0,3%, tandis que le pourcentage de chevauchement génotypique entre les différents individus est en moyenne de ∼73±2%.

Analyse PCA et tSNE des BMMC

Les données de six échantillons ont été utilisées : deux contrôles sains, AML027 pré et post-transplantation, et AML035 pré et post-transplantation. Chaque échantillon a été sous-échantillonné à ∼10k lectures cartographiées avec confiance par cellule. Ensuite, la matrice de code-barres gène-cellule de chaque échantillon a été concaténée. Les analyses PCA, tSNE et k-means ont été réalisées sur la matrice regroupée, en suivant les mêmes étapes que celles décrites dans l’analyse PCA et tSNE des PBMC. Pour le clustering k-means, K=10 a été utilisé sur la base de la courbure du scree plot de la somme des erreurs au carré.

Les gènes spécifiques aux clusters ont été identifiés en suivant les étapes décrites dans ‘Identification des gènes spécifiques aux clusters et classification basée sur les marqueurs’. La classification a été attribuée sur la base des gènes spécifiques aux clusters, et sur la base de l’expression de certains marqueurs bien connus des types de cellules immunitaires. Blasts et Immature Ery 1  » fait référence au cluster 4, qui exprime CD34, un marqueur des progéniteurs hématopoïétiques39, et Gata2, un marqueur des érythroïdes précoces40. Immature Ery 2  » fait référence aux grappes 5 et 8, qui présentent l’expression de Gata1, un facteur de transcription essentiel à l’érythropoïèse41, mais pas de CD71, que l’on trouve souvent dans les cellules érythroïdes plus engagées39. Immature Ery 3′ se réfère au cluster 1, qui montre l’expression de CD71. Mature Ery  » fait référence au groupe 2. HBA1, un marqueur des cellules érythroïdes matures, est détecté de manière préférentielle dans le cluster 2. Le cluster 3 a été désigné comme ‘Granulocytes immatures’ en raison de l’expression de marqueurs granulocytaires précoces tels que AZU1 et IL8 (réf. 42), et de l’absence d’expression de CD16. Le cluster 7 a été attribué aux ‘Monocytes’ en raison de l’expression de CD14 et FCN1, par exemple. B  » fait référence aux clusters 6 et 9 en raison de marqueurs tels que CD19 et CD79A. ‘T’ fait référence au cluster 10, en raison de marqueurs tels que CD3D et CD8A.

Data availability

Toutes les données pertinentes sont disponibles auprès des auteurs. Les données RNA-seq des cellules uniques ont été déposées dans les Short Read Archive sous le numéro d’accession SRP073767. Les données sont également disponibles à http://support.10xgenomics.com/single-cell/datasets. Le code d’analyse pour l’analyse 68k PBMC est disponible à https://github.com/10XGenomics/single-cell-3prime-paper.

.