Gènes humains échappant à l’inactivation du X révélés par les données d’expression de cellules uniques

Un cadre pour mesurer l’échappement à l’inactivation du X dans les cellules uniques

Nous identifions les échappés en analysant l’expression des gènes à partir de cellules somatiques uniques en utilisant la méthodologie scRNA-Seq (voir Méthodes). Pour évaluer la sensibilité de la méthode, nous comparons l’expression du chromosome X (ChrX) aux autres chromosomes autosomiques. Plus précisément, nous nous sommes concentrés sur le chromosome 17 (Chr17), riche en gènes, comme prototype de chromosome autosomique. Le Chr17 a été choisi car il représente un chromosome avec un nombre minimal de gènes imprimés spécifiques aux parents. Les propriétés quantitatives de ChrX et Chr17 sont énumérées dans la Fig. 1a.

Cette étude est basée sur l’analyse de deux ressources d’origine féminine : (i) les fibroblastes primaires UCF1014 (avec 104 cellules, voir Méthodes). Cet ensemble est spécifié par une couverture plus élevée de données transcriptomiques, mais manque d’informations sur la mise en phase des haplotypes (Fig. 1b) ; (ii) un ensemble plus petit de données de lymphoblastes clonaux (n = 25) de la lignée cellulaire GM12878 avec des génomes diploïdes parentaux entièrement mis en phase et séquencés (Fig. 1c). Dans les deux ensembles de données, la transcription au niveau des SNP hétérozygotes (hSNP) est la source d’information permettant de déterminer l’expression monoallélique ou bialélique. Chaque hSNP, dans chaque cellule, qui est soutenu par des preuves d’expression au-dessus d’un seuil prédéterminé est considéré comme un SNP informatif (iSNP) (voir Méthodes, Fichier supplémentaire 1 : Texte). La somme des iSNP par gène définit son étiquette unique en tant que gène inactivé ou échappé (voir Méthodes, Fig. 1b-c, Fichier supplémentaire 1 : Texte).

Quantification de l’expression bialélique à partir de fibroblastes primaires unicellulaires

Nous avons analysé les données scRNA-Seq publiées à partir de fibroblastes humains primaires femelles . Avant l’analyse, nous avons pris soin d’un écueil expérimental pertinent pour de nombreuses études unicellulaires. Ce piège concerne les doubles de cellules dans lesquelles plus d’une cellule est séquencée par une bibliothèque. Dans un tel scénario, les différents chromosomes X actifs (Xa) de deux cellules différentes qui sont incluses dans l’échantillon produiront un signal bialélique le long du chromosome X entier. Même si la fraction de doublets devrait être faible, elle peut conduire à une interprétation erronée. Par conséquent, avant d’analyser les données provenant des fibroblastes, nous avons revu les 104 fibroblastes et testé leur rapport bialélique par rapport à ChrX (voir Méthodes). Trois cellules présentaient un degré exceptionnellement élevé d’expression bialélique qui pourrait indiquer un mélange de deux chromosomes X parentaux (fichier supplémentaire 1 : texte et figure S3). Nous avons retiré les trois cellules suspectes de toutes les analyses.

Puis, pour chaque cellule unique, nous avons compté le nombre de lectures qui ont été mappées de manière unique aux allèles hSNP. Le rapport allélique (AR) pour chaque iSNP est défini comme la fraction des lectures cartographiées à l’allèle alternatif (Alt) sur le total des lectures (voir Méthodes, Fichier supplémentaire 2 : Tableau S1). La figure 2a-c résume l’AR de ChrX, Chr17 et de l’ensemble des chromosomes autosomiques selon la collection de fibroblastes primaires (101 cellules sur 104). De plus, la figure 2d montre la distribution de l’AR d’un ensemble annoté de gènes imprimés provenant de tissus de la peau (selon ). Comme indiqué précédemment, un biais dans la cartographie vers le génome de référence (AR = 0) est évident (Fig. 2a-d). De plus, une fraction substantielle d’expression monoallélique a été observée pour tous les ensembles testés (Fig. 2a-d). Cette apparition dominante de l’expression monoallélique dans les cellules uniques est causée par une combinaison à la fois d’un sous-échantillonnage des transcrits, et d’un phénomène connu sous le nom de « rafale transcriptionnelle » .

Nous nous sommes concentrés uniquement sur les iSNP qui présentent une signature non monoallélique (c’est-à-dire excluant AR = 0 et AR = 1). Nous avons observé une différence marquée dans la distribution des AR de ChrX et des gènes imprimés par rapport à Chr17 et à tous les chromosomes autosomiques (comparer les Fig. 2a à b et les Fig. 2c à d). En conséquence, plusieurs observations peuvent être tirées des résultats présentés dans la Fig. 2 : (i) Chr17 et tous les autosomes partagent un profil AR similaire. (ii) Une tendance claire vers une expression équilibrée (AR = 0,5) est apparente pour tous les chromosomes autosomiques (Fig. 2a-b), mais pas pour ChrX ou les gènes imprimés (Fig. 2c-d). (iii) La fraction d’expression non-monoallélique dans les chromosomes autosomiques est plus élevée (~ 18%) par rapport à ChrX (~ 9%). (iv) La fraction d’expression non-monoallélique dans les gènes imprimés montre un niveau intermédiaire (13%). Ce niveau intermédiaire est probablement le reflet de l’incohérence inhérente à l’identité des gènes imprimés. Fichier supplémentaire 3 : Le tableau S2 énumère les iSNP de soutien pour tous les chromosomes analysés dans les fibroblastes, y compris l’ensemble des gènes imprimés.

Identification des évasions dans les fibroblastes primaires unicellulaires

Dans l’ensemble de données des fibroblastes primaires, il y a 232 et 485 gènes qui sont soutenus par des preuves iSNP pour ChrX et Chr17, respectivement. Comme ces cellules manquent d’informations sur la mise en phase du génome (Fig. 1b), les informations sur l’échappement du chromosome X sont limitées à l’ensemble des iSNP bialéliques (voir Fig. 1b). Nous avons agrégé les iSNP en fonction de leurs gènes correspondants (Fig. 1b). L’agrégation est effectuée sur différentes cellules uniques et sur plusieurs iSNP au sein d’une paire cellule-gène spécifique. Un gène sera étiqueté escapee candidate lorsqu’il est associé à plusieurs iSNP bialéliques. Au total, nous avons identifié 24 gènes de ce type (tableau 1) qui représentent 10,3 % de tous les gènes exprimés dans ChrX. Comme prévu, la fraction de gènes sur Chr17 présentant une expression bialélique est sensiblement plus élevée (49,3%, fichier additionnel 4 : tableau S3).

Tableau 1 énumère les candidats à l’évasion à divers degrés de soutien. Par exemple, les gènes ZFX (Zinc finger X-chromosomal protein) et SMC1A (Structural maintenance of chromosomes protein 1A) sont fortement soutenus avec 103 et 19 iSNP bialéliques, respectivement. Une autre augmentation de la fiabilité de l’identification des évadés est basée sur le fait d’avoir au moins 2 cellules indépendantes qui ont fourni des informations sur l’expression bialélique. Nous montrons que 21 gènes sur 24 répondent à ce critère strict (tableau 1). Notamment, parmi les escapees identifiés, nous n’avons détecté que trois gènes PAR (SLC25A6, CD99, et DHRSX, Tableau 1). L’attribution de ces gènes en tant que gènes échappés correspond à l’expression attendue de la RAP. D’après le nombre de gènes PAR bialéliques parmi les gènes PAR exprimés, nous avons estimé que le taux de découverte de faux négatifs pour les gènes échappés pouvait atteindre 70 % (c’est-à-dire que nous avons manqué 7 des 10 gènes PAR exprimés). Fichier additionnel 4 : Le tableau S3 montre le support pour le tableau 1.

Quantification de l’expression allélique à partir de lymphoblastes à phase clonale

Une limitation majeure du protocole décrit ci-dessus concerne l’absence de mise en phase des haplotypes parentaux. Dans ce cadre, les iSNP ne peuvent pas être assignés à Xa ou Xi. Par conséquent, le choix aléatoire de Xi qui caractérise les cellules primaires limite le taux de découverte des échappées. Nous avons étendu l’analyse de scRNA-Seq aux lymphoblastes d’origine féminine de la lignée cellulaire clonale GM12878. Voir la figure 1c et le fichier supplémentaire 1 : Figure S4.

La figure 3a montre le profil d’expression des cellules uniques de lymphoblastes clonaux (n = 25) (fichier supplémentaire 2 : tableau S1, fichier supplémentaire 5 : tableau S4). Dans toute cellule unique, l’expression monoallélique reflète la combinaison d’un sous-échantillonnage des transcrits et du phénomène d’éclatement transcriptionnel. Il est clair que l’expression maternelle provenant du Xa domine (Fig. 3a, en haut). Une observation qui concorde avec l’origine maternelle Xa rapportée de la lignée cellulaire GM12878 . Dans la plupart des cellules analysées, une fraction faible mais substantielle de l’expression totale mesurée provient du chromosome paternel Xi (Fig. 3a, haut). En revanche, le Chr17 et les chromosomes autosomiques montrent une expression égale des deux allèles (Fig. 3a, panneaux du milieu et du bas).

D’après la figure, il est évident que le phénomène d’éclatement transcriptionnel affecte tous les chromosomes y compris ChrX. Pour évaluer l’impact de ce phénomène sur l’identification des gènes échappés, nous avons comparé des cellules individuelles par rapport à un pool de cellules (Pool100, Fig. 3a, barre de droite). Alors que la majorité des iSNP de Chr17 présentent des profils bialéliques, ChrX reste dominé par une expression monoallélique maternelle.

La figure 3b est une vue agrégée de ChrX, Chr17 et des chromosomes autosomiques. Les données sont basées sur 375, 808 et 20 212 hSNP exprimés, respectivement. La figure 3b (panneaux du milieu et de droite) montre une répartition égale des allèles parentaux de Chr17 et des chromosomes autosomiques (figure 3b, en haut). En effectuant la même analyse sur les données recueillies à partir du pool 100 (Fig. 3c), on constate que la partition des allèles parentaux reste pratiquement inchangée (comparer la fraction occupée par les couleurs rose et bleue, Fig. 3b-c). De plus, nous avons observé un passage d’une expression monoallélique (Fig. 3b-c, couleur remplie) à une expression bialélique (Fig. 3b-c, couleur dépouillée). La fraction de l’expression bialélique pour Ch17 a augmenté de 19% dans les cellules uniques à 80% dans le Pool100, et pour les chromosomes autosomiques de 18 à 79% (Fig. 3c, panneaux central et droit). Les résultats de la Pool100 indiquent que l’expression monoallélique observée dans les cellules uniques est pratiquement abolie en faisant la moyenne du signal.

Les résultats de ChrX (Fig. 3b (gauche) sont fondamentalement différents par rapport à Chr17 ou aux chromosomes autosomiques (Fig. 3b-c). La différence la plus notable est que seulement 21% des iSNP exprimées sont associées à l’allèle Xi paternel dans ChrX (Fig. 3b, en haut). En outre, lors de l’analyse du Pool100, la fraction d’expression bialélique reste limitée (passage de 9 % dans les cellules individuelles à 34 % dans le Pool100). Le modèle observé de ChrX à partir du Pool100 (Fig. 3c, à gauche) s’explique mieux par une moyenne du signal monoallélique stochastique (au même degré que les autres chromosomes) tout en maintenant un fort signal de l’expression monoallélique de Xa. Voir le fichier additionnel 5 : Tableau S4 pour le ratio allélique des lymphoblastes de tous les chromosomes testés et le Pool100.

Identification des évasions des lymphoblastes unicellulaires

La figure 4a est une vue centrée sur les gènes qui montre la partition allélique iSNP des lymphoblastes (colorée selon leur origine comme expression maternelle, paternelle ou mixte, voir Méthodes). Seul le sous-ensemble de gènes qui sont soutenus par plusieurs iSNP est listé selon leur ordre le long des chromosomes. Au total, nous signalons 93 gènes annotés sur ChrX (Fig. 4a, 30 gènes échappés et 63 gènes inactivés). Notez que les gènes inactivés par le X représentent les gènes qui sont exprimés principalement à partir du Xa maternel. Un groupe de gènes avec une expression paternelle à l’extrémité du bras p de ChrX représente l’expression bialélique attendue des gènes PAR (Fig. 4a). Des preuves supplémentaires de l’expression paternelle sont localisées au XIC avec des gènes tels que XIST, JPX et FTX. Alors que la plupart des échappées sont soutenues par un nombre limité d’iSNP, quelques-unes d’entre elles, comme ZFX, CD99 et SLC25A6, sont soutenues par un nombre relativement important d’iSNP de soutien (48, 38 et 34, respectivement).

Une autre méthode pour évaluer l’ampleur du phénomène d’inactivation du X consiste à quantifier la preuve directement à partir de la somme de toutes les lectures séquencées (abrégé en protocole basé sur les lectures). Les figures 4b-c comparent le nombre de lectures de Chr17 (figures 4b) et de ChrX (figure 4c) selon l’origine paternelle ou maternelle. Nous comparons les données d’expression des cellules uniques et du Pool100. La régression linéaire pour l’expression des gènes de Chr17 montre une ligne d’ajustement de corrélation élevée (r2 = 0,823, Fig. 4b). Comme prévu, la corrélation est plus forte dans les données provenant du Pool100 (r2 = 0,946, Fig. 4b). Nous avons conclu que malgré l’expression monoallélique due au phénomène d’éclatement transcriptionnel, l’expression allélique équilibrée de tous les gènes est fortement soutenue. Pour ChrX, cependant, la régression linéaire résultante des cellules individuelles est faible (r2 = 0,238, Fig. 4c), et n’a pas été améliorée par les données du Pool100 (r2 = 0,222, Fig. 4d). L’inspection des données d’expression pour ChrX montre que les lignes de régression penchent en fait vers l’expression de la Xa maternelle (axe des x). Les données d’expression sont cohérentes avec deux lignes de régression distinctes pour ChrX. L’une qui correspond aux gènes inactivés (parallèle à l’axe des x), et l’autre correspond à une expression bialélique de confiance.

L’application du protocole conservateur basé sur les iSNP conduit à l’identification de 30 gènes comme candidats à l’évasion qui sont également soutenus par le protocole basé sur la lecture (tableau 2). Le protocole basé sur la lecture (c’est-à-dire étiquetage d’un gène comme échappé en ayant un nombre minimal de lectures paternelles, voir Méthodes) a élargi la liste des candidats échappés pour inclure globalement 49 gènes (fichier supplémentaire 4 : tableau S3).

Tester l’origine parentale des allèles le long d’un gène dans la même cellule est un test strict pour la fiabilité des iSNP. Ce test n’est valable que pour les gènes avec plusieurs hSNP. Ces gènes qui sont soutenus par deux hSNP ou plus exprimés représentent 44% des gènes. Nous considérons qu’un gène est cohérent si l’expression le long du gène dans une cellule spécifique n’est pas monoallélique aux deux allèles. Au total, nous avons identifié 3 gènes incohérents – TEX11, FTX et ZCCHC16. Pour 6 autres gènes, l’incohérence n’est que partielle car il existe d’autres observations d’expression bialélique. L’estimation de l’incohérence totale (3 gènes sur les 29 qui étaient éligibles pour ce test) suggère que la limite supérieure d’une interprétation erronée est de 10%. En outre, le protocole basé sur les iSNP a identifié 9 des 11 gènes PAR exprimés. Ainsi, nous extrapolons que le taux de détection des échappées est de 82 %. Il est intéressant de noter que l’analyse du Chr17, en supposant qu’il n’y a pas de biais allélique systématique, a montré que 7,3 et 9,6 % des gènes étaient associés à une expression monoallélique maternelle et paternelle, respectivement. Ces résultats fournissent une limite supérieure de 17,9% à la probabilité de faux étiquetage de gènes dans Chr17 et peuvent être utilisés pour estimer la limitation de la méthode.

Exiger une preuve paternelle d’au moins deux cellules réduit le nombre d’évasions de 49 à 18 (dont 5 gènes PAR). De nombreux gènes de preuve unique résultant d’une cellule aberrante spécifique (SRR764803) qui a contribué à 17 des 49 gènes signalés (tableau 2, fichier supplémentaire 5 : tableau S4). Cette cellule présente les niveaux d’expression les plus élevés sur ChrX et Chr17 (fichier supplémentaire 1 : texte et figure S4). Le tableau 2 énumère les candidats évadés découverts ainsi que les preuves à l’appui (tableau 2).

Comparaison des évadés identifiés avec les connaissances actuelles

Nous avons testé la correspondance entre les évadés identifiés de notre étude et un catalogue basé sur la littérature . Ce catalogue unifié a été compilé à partir de l’intégration de quatre études indépendantes qui couvre 1144 gènes de ChrX. Les gènes de ce catalogue sont répartis manuellement en neuf catégories définies (voir Méthodes). La plus grande catégorie représente les gènes qui manquent d’information (45%) . Environ 15% des gènes (168/1144) sont considérés comme ‘associés à des évadés’ (voir Méthodes). Nous considérons l’ensemble compilé de gènes  » associés aux évadés  » comme un étalon-or pour tester le taux de découverte des évadés dans notre étude (total de 124 gènes, gènes PAR exclus, collectivement appelés Balaton-Esc).

Nous avons appliqué un test statistique hypergéométrique (voir Méthodes) pour évaluer le chevauchement des différentes listes de gènes évadés (figure 5). La figure 5a montre le nombre d’évadés identifiés dans les fibroblastes et les lymphoblastes (à l’exclusion des gènes PAR). Notez que seuls les gènes inclus dans le repère Balaton-Esc sont inclus dans cette analyse (Fig. 5). La figure 5b montre la signification statistique du chevauchement entre les listes de gènes de la figure 5a et le repère Balaton-Esc . Comme on peut le voir, il existe un chevauchement significatif entre les gènes échappés des lymphoblastes (tableau 2 et fichier additionnel 4 : tableau S3) et la liste Balaton-Esc (Fig. 5b, valeur p = 7,43E-8). L’application du même test pour les fibroblastes primaires (tableau 1) a donné lieu à une signification plus faible (valeur p = 4,07E-2).

La figure 5c représente les gènes de chevauchement entre les escapees identifiés dans notre étude et le catalogue Balaton-Esc (168 gènes incluant les gènes PAR) . Nous avons également inclus une ressource complémentaire basée sur 940 transcriptomes provenant de scRNA-Seq (25 gènes d’évasion, Tukiainen-Esc). Le diagramme de Venn montre que chacune des études ci-dessus contribue aux connaissances actuelles sur les évadés. Les gènes échappés des deux ressources externes se chevauchent pour 18 des 25 gènes signalés (72 %). Comme le montre la figure 5c, 62 % des évasions rapportées par les lymphoblastes se chevauchent avec les listes d’évasions externes, tandis que les fibroblastes ne sont soutenus que par 38 % de chevauchement. Notamment, la plupart de nos gènes candidats d’échappement découverts à partir des fibroblastes (62%) n’ont aucune correspondance avec les autres listes testées (fichier supplémentaire 6 : tableau S5).

LncRNAs étendent la liste des candidats d’échappement

Nous avons étendu notre analyse aux longs ARN non codants (lncRNAs) en utilisant les mêmes critères d’échappement que ceux utilisés à partir des gènes codants (Fig. 6). Au total, nous avons identifié 15 lncRNAs comme candidats à l’évasion, parmi lesquels seuls quelques-uns avaient été étudiés précédemment. La localisation des lncRNAs et des gènes codants échappés le long de ChrX est présentée (Fig. 6a). Nous avons testé les positions des échappées le long du ChrX par rapport à tous les gènes du ChrX. Alors que la distribution positionnelle des échappées de lncRNA est similaire (test de Kolmogorov-Smirnov, p-value = 0,57), elle est différente pour les échappées de codage (test de Kolmogorov-Smirnov, p-value = 0,004, Fig. 6a).

La figure 6b énumère tous les 15 gènes lncRNA échappés identifiés, parmi eux, sont des gènes lncRNA du XIC qui coordonnent l’activation et le maintien de l’inactivation du X. De nombreux lncRNA sont localisés dans des segments actifs sur le plan transcriptionnel (par exemple, dans le PAR ou le XIC), tandis que d’autres sont localisés dans des régions non conservées qui sont enrichies en lncRNA longs et courts. Pour les lncRNA supplémentaires, y compris les gènes inactivés, voir le fichier additionnel 7 : tableau S6.

Partition basée sur les preuves des gènes d’échappement

La figure 6c résume la partition selon les preuves pour tous les gènes qui sont signalés avec n’importe quel niveau de confiance comme des échappatoires (fichier additionnel 6 : tableau S5). Cette liste comprend 75 candidats qui sont signalés dans cette étude, y compris les preuves de Pool100, une collection de nouveaux lncRNAs. La figure 6d fournit une liste détaillée des résultats obtenus à partir de ces données au regard de nombreuses ressources externes. Le nombre de cellules fournissant des preuves est également indiqué pour chaque gène.

En tenant compte de ces ressources externes, nous avons apparié chaque gène en fonction de la qualité du support indépendant qui lui est associé (Fig. 6c et d). Plus précisément, nous avons divisé la liste des 75 gènes en quatre groupes : (i) Les gènes sont étiquetés « confirmés » s’ils sont signalés comme échappés par les deux ressources externes précédemment discutées . Il y a 31 gènes de ce type, pour lesquels cette étude fournit une confirmation supplémentaire de leur identité d’évadés. (ii) Sept autres gènes sont étiquetés « approuvés ». Ces gènes sont marqués comme échappés par une seule des deux ressources externes. Pour ces gènes, la preuve indépendante de cette étude approuve leur identité. (iii) 12 gènes supplémentaires sont marqués comme des évasions « probables ». Ces gènes sont attribués en fonction de l’accord avec un rapport externe supplémentaire portant sur 114 gènes échappés, qui n’était pas inclus dans le catalogue basé sur la littérature, et qui peut donc être considéré comme une ressource indépendante. (iv) 25 autres gènes sont marqués comme des évasions « possibles ». Ces gènes n’ont pas d’évidence de littérature pour soutenir leur identité, ainsi leur assignation en tant qu’évadés reste moins soutenue. Parmi eux, 5 gènes ne sont soutenus que par la cellule aberrante du lymphoblaste et sont donc probablement faux. Dans cet ensemble, il y a 10 autres gènes qui n’ont pas été signalés par l’une des trois ressources externes discutées, La majorité des lncRNAs négligés appartiennent à ce groupe. Deux de ces lncRNAs échappés (TCONS_00017125, et TCONS_00017281, Fig. 6b) sont situés dans la région PAR qui est exceptionnellement active dans la transcription bialélique. Les données sources de la figure 6d sont dans le fichier additionnel 6 : tableau S5.