Primary fibroblast cultures

Our protocol for generating primary skin fibroblast cultures from a skin biopsy is based on Vangipuram et al.59 From adults (18+ years of age) at the NIH Clinical Center we obtained two 4-mm skin punch biopsies from the upper arm, which were immediately placed into a 15 mL conical tube with 10 mL of media (DMEM/F12 (Gibco), 20% FBS, and 100 IU/mL Penicillin-Streptomycin. Tubes were shipped to the Page lab overnight on ice for processing. Each biopsy was used to generate a separate skin fibroblast culture. Biopsies were cut into 18 pieces of equal size and placed 3/well in gelatinized 6-well plates with 1 mL media (High Glucose DMEM (Gibco), 20% FBS, L-Glutamine (MP Biomedicals), MEM Non-Essential Amino Acids (Gibco), 100 IU/mL Penicillin/Streptomycin (Lonza)). Plates were gelatinized by incubating 1 mL sterile 0.1% gelatin (Sigma) solution per well for 30 min at room temperature.

Plates were incubated for 1 week at 37°C without disturbance to allow biopsies to attach to the plate and begin to grow out. During week 2, we added 200 μL of fresh media per well every 2-3 days, being careful not to disturb the biopsies. The following week (week 3), we aspirated the media and replaced with 1 mL fresh media per well every 2-3 days. During week 4, we aspirated the media and replaced with 2 mL fresh media per well every 2-3 days. At this point, the fibroblasts generally reached the edges of the wells and were expanded to two T75 gelatinized flasks per 6 well plate. After two days, we combined the cells from the two T75 flasks and split them to three T175 gelatinized flasks. After two days, cells were viably frozen with 1 million cells per vial in freezing media (fibroblast culture media + 5% DMSO). Cells were also preserved for RNA extraction (see below). During optimization of the protocol, cell culture purity was confirmed by immunofluorescence of SERPINH1, a fibroblast marker.

Cell collection for subsequent analysis

Cells were collected when LCL cultures reached 30mL, and fibroblasts were ~80% confluent in three T175 plates. All cell counting was performed using the Countess II cell counter (Life Technologies) and Trypan Blue exclusion. Cultures with >85% cell viability were used in subsequent experiments. To preserve cells for subsequent RNA extraction, 1 million cells were washed in PBS, pelleted, and resuspended in 500 μL TRIzol (Invitrogen) or 200 μL RNAprotect Cell Reagent (Qiagen). Cell suspensions were then frozen at −80°C. Cell cultures were maintained at low passage number; RNA-sequencing experiments were performed on samples at or below passage 4.

Periodically, and on each passage used for experiments, cell cultures were confirmed negative for mycoplasma contamination using either the MycoAlert Kit (Lonza) following the manufacturer’s instructions, or PCR using SapphireAmp Fast PCR Master Mix (Takara) and the following primers:

Myco2(cb): 5′ CTTCWTCGACTTYCAGACCCAAGGCAT-3′

Myco11(cb): 5′ ACACCATGGGAGYTGGTAAT-3′

PCR for GAPDH was performed on the same sample, using the following primers:

hGAPDH-F: TGT CGC TGT TGA AGT CAG AGG AGA

hGAPDH-R: AGA ACA TCA TCC CTG CCT CTA CTG.

Known mycoplasma positive and negative samples were used as a reference.

RNA extraction, library preparation, and sequencing

RNA was extracted from 1 million cells per experiment using the RNeasy Plus Mini Kit (Qiagen) following the manufacturer’s instructions, with the following modifications: Cells in RNAprotect Cell Reagent were thawed on ice, pelleted, and lysed in buffer RLT supplemented with 10μL β-mercaptoethanol per mL. For most samples, ERCC control RNAs were added to the lysate based on the number of cells: 10μL of 1:100 dilution of ERCC control RNAs was added per 1 million cells. The lysate was then homogenized using QIAshredder columns (Qiagen), and transferred to a gDNA eliminator column. All subsequent optional steps in the protocol were performed, and RNA was eluted in 30 μL RNase-free water. RNA levels were measured using a Qubit fluorometer and the Qubit RNA HS Assay Kit (ThermoFisher). Before we switched to the per-cell spike-in protocol, we prepared 18 samples in which ERCC control RNAs were added based on amount of RNA after isolation: 2 μL of a 1:100 dilution of ERCC control RNAs was added per 1 μg of RNA. These samples are: #2237, 2245, 6312, 711, 4032, 706, 3429, 3430, 3442, 2690, 2703, 3107, 5297, 5566, 5755, 6029, 2547, and 525. RNA quality control was performed using the 5200 Fragment Analyzer System (Agilent); we consistently purified high-quality RNA with RNA integrity numbers (RIN) near 10. We randomized the samples by karyotype into batches for RNA extraction, library preparation, and sequencing.

RNA sequencing libraries were prepared using the TruSeq RNA Library Preparation Kit v2 (Illumina) with modifications as detailed in Naqvi et al,60 or using the KAPA mRNA Hyper-Prep Kit V2 (Roche). In both cases, libraries were size selected using the PippinHT system (Sage Science) and 2% agarose gels with a capture window of 300-600 bp. Paired-end 100x100 bp sequencing was performed on a HiSeq 2500 or NovaSeq 6000 (Illumina). Table S1 lists the library preparation kit and sequencing platform for each sample.

RNA-seq data processing and analysis

All analyses were performed using human genome build hg38, and a custom version of the comprehensive GENCODE v24 transcriptome annotation.52 This annotation represents the union of the “GENCODE Basic” annotation and transcripts recognized by the Consensus Coding Sequence project.61 Importantly, the GENCODE annotation lists the PAR gene annotations twice – once on Chr X and once on Chr Y– which complicates analysis. We removed these annotations from Chr Y so the PAR genes are only listed once in our annotation, on Chr X. To analyze samples in which ERCC spike-ins were added, we merged our custom transcript annotation with the ERCC Control annotation.

Reads were pseudoaligned to the transcriptome annotation, and expression levels of each transcript were estimated using kallisto software v0.42.5.54 We included the “--bias” flag to correct for sequence bias. The resulting count data (abundance.tsv file) were imported into R with the tximport package v1.14.062 for normalization using DESeq2 v1.26.0.55 For downstream analysis, we used only protein-coding genes (as annotated in ensembl v104) with the following exceptions: we included genes annotated as pseudogenes on Chr Y that are members of X-Y pairs (TXLNGY, PRKY) and well-characterized long non-coding RNAs (lncRNAs) involved in X-inactivation or other processes (XIST, JPX, FTX, XACT, FIRRE, TSIX). We annotated genes distal to XG, which spans the pseudoautosomal boundary on Xp and is truncated on Chr Y, as part of PAR1 - 15 genes in total. PAR2 comprised the four most distal genes on Xq and Yq. Annotations of non-pseudoautosomal region of the X (NPX) genes with homologs on the non-pseudoautosomal region of the Y (NPY) were derived from Bellott et al.63 224 protein-coding genes on Chr 21 (ensembl v104) were used as a starting point for our analyses. We excluded 21 annotated genes in several regions with high homology between the long and short arms of Chr 21 because the assembly was not fully validated in these regions (https://www.ncbi.nlm.nih.gov/grc/human/issues?filters=chr:21).

Identifying genes affected by changes in chr X, Y, or 21 copy number

We first defined lists of expressed NPX, NPY, PAR, or Chr 21 genes as those with median TPM of at least 1 in 46,XX or 46,XY samples. To ensure that no genes with robust expression were excluded, we also analyzed LCL and fibroblast expression data from GTEx,53 and included several genes that were just below our TPM cutoff but had median TPM of at least 1 in those datasets.

For each expressed NPX, NPY, or PAR gene we performed linear modeling using the lm() function in R. These calculations suppose that each additional chromosome adds a consistent and equal increment to the total expression level of the gene in question.

For NPX and PAR genes we used the following equation:

E represents the expression (read counts) per gene, ${\beta }_0$ represents the intercept, ${\beta }_X$ and ${\beta }_Y$ are the coefficients of the effect of additional copies of Chr Xi or Y, respectively, and $\in$ is an error term. For this equation, the intercept represents the 45,X samples.

For NPY genes we employed the following equation, analyzing only those samples with one or more copies of Chr Y:

For this equation, the intercept represents the 46,XY samples.

For Chr 21 genes we employed the following equation, analyzing only those samples with 46,XX; 46,XX; 47,XY,+21; or 47,XX,+21 karyotypes:

${\beta }_{21}$ and ${\beta }_{Sex}$ are the coefficients of the effect of an additional copy of Chr 21 and sex (XY vs XX), respectively. For this equation, the intercept represents the 46,XX samples.

The resulting p values were adjusted for multiple hypothesis testing using the p.adjust() function in R, specifying the Benjamini Hochberg method. Genes with a false discovery rate (FDR) < 0.05 were considered significant. To compute the normalized expression change per Chr Xi (ΔE_X) or Y (ΔE_Y), we divided the coefficient of interest (${\beta }_X$ or ${\beta }_Y$) by the average intercept across batches, which corresponds to the baseline expression of the gene in samples with only one X chromosome (for NPX and PAR genes) or one Y chromosome (in the case of NPY genes). For Chr 21, we computed ΔE₂₁ by dividing the coefficient (${\beta }_{21}$) by the average intercept across batches divided by two to obtain the average expression from one copy of Chr 21.

In the case of XIST, which is only expressed when two or more copies of Chr X are present, we used the following equations:

We calculated the standard error (SE) of ΔE_X, ΔE_Y, and ΔE₂₁ using the following equations:

To confirm the validity of our approach, we used bootstrapping to sample our dataset with replacement 1000 times and obtained similar results. BEX1 was removed from downstream analyses in fibroblasts because two samples (one 45,X and one 49,XXXXY) had high expression values for this gene resulting in >25 times higher error values for ΔE_X and ΔE_Y compared to all other genes.

Saturation analysis for sex chromosome-encoded genes

For LCLs and fibroblasts, size-n subsets of available RNA-seq libraries were sampled randomly without replacement, 100 times for each sample size, n. After confirming that the model matrix would be full rank in each sampling (for example, that samples would not all be of the same karyotype or batch), we performed linear modeling on NPX, PAR, NPY, and Chr 21 genes as described above to identify genes whose expression changes significantly (FDR < 0.05) with copy number of Chr X, Y or 21.

Assessing linearity of sex-chromosome gene expression changes

To assess whether sex-chromosome gene expression changed linearly (i.e., by a fixed amount) with additional X or Y chromosomes, their expression levels across the LCL or fibroblast samples were fit by non-linear least squares to the power curves shown below, using the “nlsLM” function from the R package “minpack.lm”.NPX genes:

PAR genes:

NPY genes:

In each of the equations above, $y_j^{\ast }$ is the normalized RNA-seq read count for a given gene in sample $j$, given by the raw read count in sample $j$ divided by the average read count in the set of samples $S_{\{·\}}$ with only one chromosome of the relevant type: for NPX genes, 1 copy of Chr X (and any number of Y chromosomes); for PAR genes, 45,X samples; for NPY genes, 1 copy of Chr Y (and any number of X chromosomes). $b=0.5$ and $a=1$ were used as initial parameter values. Fitted values of $a\approx 1$ indicate a linear relationship between expression and sex-chromosome count. Fitted values of $a\approx 0$ or $b\approx 0$ indicate no change in expression with X or Y count.

ΔE_X calculations in samples with 0 Y chromosomes (females) and 1 Y chromosome (males)

We took subsets of the samples with either zero Y chromosomes (females) and with one Y chromosome (males) and performed the same linear modeling and ΔE_X calculations as above. We removed MAP7D2 in female LCLs, IL13RA2 in female fibroblasts, and FHL1 in male fibroblasts because their error values (likely due to smaller sample size) were much higher than those of other genes. To compare the linear modeling results, we performed Pearson correlations between the results using all samples, and those from male-only or female-only samples.

Reanalysis of array data and comparison to RNA-seq data

A previous study performed gene expression analysis, using Illumina oligonucleotide BeadArrays, of LCLs from 68 individuals of the following karyotypes: 45,X; 46,XX; 46,XY; 47,XXX; 47,XXY; 47,XYY; and 48,XXYY.24 Since this microarray dataset was generated from an independent set of samples, we sought to validate our results through a reanalysis of the data.

The raw data from the microarrays was not publicly available, but the authors provided us pre-processed data upon request, which we used to perform our analysis. To identify genes that cleared a minimum signal threshold to be considered expressed in the microarray data, we assessed the median signal in 46,XY samples for all Chr Y genes annotated on the microarray. We focused on Chr Y genes in this analysis because many are known to be expressed exclusively in testes, and therefore could provide us with an appropriate sense of the background signal expected for genes not expressed in LCLs. From this analysis, we concluded that a signal threshold of 111 would be appropriate for identifying expressed genes (Figure S9A). We used this threshold to identify 278 expressed Chr X genes (including PAR and NPX genes) in the microarray dataset. This was fewer than the 341 expressed Chr X genes identified in our LCL RNA-seq data, but more than double the 121 expressed Chr X genes reported in Raznahan et al. (Table S4 in Raznahan et al.).24 This discrepancy could not be resolved by simply increasing the signal threshold in our analysis, as TMSB4X, one of the most highly expressed genes in LCLs, was excluded from the previously reported list of expressed genes.

Using our list of 278 expressed genes from the Raznahan et al. dataset, we analyzed the microarray signal values (in place of RNA-seq read counts) using linear models as a function of Xi copy number, controlling for Chr Y copy number. We calculated ΔE_X values from the microarray data and compared these to our RNA-seq dataset using a Pearson correlation, which revealed that the results were generally concordant (Figure S9B). For genes that were lowly expressed, however, ΔE_X values tended to be much lower in the microarray dataset, consistent with the higher sensitivity of RNA-seq data.

Isoform-specific analysis of RNA-seq data

After estimating counts for each transcript using kallisto software (as described above, with 100 bootstraps), we used sleuth v0.30.0 to normalize those transcript counts.56 X chromosome transcripts were called expressed if their corresponding gene was on the list of expressed genes (above) and median transcript counts were >200. Linear regressions and ΔE_X calculations for transcripts were performed as for genes (above) to identify transcripts whose abundance changes significantly with additional copies of Chr X.

The following ENCODE datasets were used for visualization in IGV software58 at the UBA1 locus.

Gene constraint analysis

To investigate sensitivity to a reduction in gene dosage, we used three metrics: LOEUF, RVIS, and pHI. We downloaded LOEUF (loss-of-function observed/expected upper fraction) scores from gnomAD (v2.1.1.lof_metris.by_gene.txt; https://gnomad.broadinstitute.org/), and only used scores with a minimum of 10 expected LoF variants. Updated RVIS (residual variation intolerance scores)36 including the ExAC dataset were downloaded from http://genic-intolerance.org/data/RVIS_Unpublished_ExAC_May2015.txt. Updated probability of haploinsufficiency (pHI) scores37 were downloaded from https://www.deciphergenomics.org/files/downloads/HI_Predictions_Version3.bed.gz. To complement these data, we obtained a list of genes with observed homozygous loss-of-function variants.35 For sensitivity to an increase in gene dosage, we used the per-gene average probability of conserved miRNA targeting scores (P_CT).38

For each metric, we computed a percentile rank score, ranking from most-to least-constrained. Because several of the metrics calculated scores separately for autosomal (including PAR genes) and NPX genes, we ranked autosomal (and PAR) genes separately from NPX genes. All annotated genes, regardless of expression status in LCLs or fibroblasts, were included in the rankings, with the following exceptions: 1) NPX genes previously annotated as “ampliconic”,64^,65 since constraint metrics cannot be accurately applied to these highly similar genes, and 2) genes with <2 annotations across all metrics.

To obtain an aggregate sense of a gene’s expression constraint across multiple metrics, we calculated the average. Among NPX genes, we considered those with |ΔE_X ≥ 0.1| (FDR < 0.05) to be most likely to contribute to phenotypes mediated by Xi copy number, prioritizing the top ten genes by the average gene-constraint metric. For PAR1 genes, we prioritized genes with an average gene constraint percentile ranking of at least 50%. To assess the phenotypic roles of highly constrained genes, we annotated them for disease phenotypes with known molecular basis from Online Mendelian Inheritance in Man (OMIM).66

Comparisons to published annotations of X-inactivation status

We re-compiled XCI status annotations of individual genes from four studies of Chr X allelic ratios.14^,15^,16^,18 Previous XCI status compilations9 incorporated DNA methylation data, which we excluded because it does not directly measure Xi transcription. Previous compilations also incorporated information about expression in human-rodent hybrid cell lines carrying a human Xi13; we incorporated this information only where allelic ratios in human cells were not available. Our final XCI status annotations are listed in Table S6, with the workflow for generating these annotations explained below.

The first dataset that we incorporated was derived from paired genomic and cDNA SNP-chips in skewed LCL and fibroblast cell cultures (Additional file 7 from Cotton et al.).14 We used the AR values provided (average Xi expression column) for genes informative in at least 5 samples, resulting in AR values for 424 genes. Using the provided numbers of informative samples and standard deviations of AR values, we computed 95% confidence intervals for the AR values. We considered a gene “Subject” to XCI if the AR 95% confidence interval included zero or the AR value was <0.1; otherwise we considered the gene to “Escape”.

The second dataset that we incorporated was derived from bulk or single cell RNA-seq of LCLs.15 The bulk RNA-seq was from an individual in the GTEx dataset with 100% skewed XCI across the body (Table S5 from Tukiainen et al.).15 The single-cell RNA-seq in LCLs was from three individuals (Table S8 from Tukiainen et al.; we excluded data from one dendritic cell sample).15 For each dataset, we calculated an AR for each gene using read counts from the more lowly and highly expressed alleles in each sample, and used the provided adjusted p-values to identify genes with significant Xi expression (padj < 0.05). For a gene to be considered informative, we required data from at least two individuals in the single cell dataset, or one individual in the single cell dataset and informative in the bulk RNA-seq dataset, resulting in 82 informative genes. We called a gene as “Subject” to XCI if there was no significant expression from Xi in either the bulk or single-cell datasets, and “Escape” if one or both of the datasets showed evidence of Xi expression.

The third dataset that we incorporated was derived from single-cell allelic expression in fibroblasts.16 The dataset includes five individuals (Dataset 3 from Garieri et al.)16 and we required data from at least two samples to be considered informative for a given gene, resulting in 203 genes. We converted their reported values (Xa reads/total reads) to AR values using the following formula: $AR=\frac{1}{\frac{Xareads}{totalreads}}-1$. We used the AR threshold calculated in the previous study16 to consider a gene significantly expressed from Xi in each sample (AR > 0.0526). If a gene had no samples with significant expression from Xi or a mean AR value < 0.1 across samples, we considered it “Subject” to XCI; otherwise, it was judged to “Escape” XCI.

The fourth dataset that we incorporated was derived from allele-specific bulk RNA-seq performed on 136 samples with skewed XCI from the set of GEUVADIS LCLs (Tables S4 and S5 from Sauteraud et al.).18 For a gene to be considered scorable, we required that it be informative in at least 10 samples, resulting in 215 genes. We calculated an AR for each gene in each sample using the read counts from the more lowly and highly expressed alleles in each sample, adjusting for the level of skewing in each sample. To identify genes that were significantly expressed from Xi across samples, we performed paired, two-sample, one-sided t tests using the t.test function in R, asking whether the raw (pre-adjusted for skewing) AR values were greater than the baseline AR given the level of skewing in each sample ($baselineAR=\frac{1-skewingcoefficient}{skewingcoefficient}$); we corrected the resulting p values for multiple comparisons with the p.adjust function in R using the Benjamini-Hochberg method. Genes were considered “Escape” if they had padj < 0.01, and “Subject” otherwise.

Next, we synthesized the calls from these four datasets. We assigned a gene as “Subject” if: 1) all studies were “Subject” or 2) most (>50%) studies were “Subject” and the average AR across all studies was <0.1. We assigned “Escape” if 1) most (>50%) studies were “Escape” or 2) 50% or fewer (but more than 0) studies were “Escape” and either i) there was more than one study with evidence of escape or ii) the average AR across all studies was ≥1. Finally, we assigned “No call” if the gene was not informative in any of the four datasets. For these genes, we investigated whether there were any calls using hybrids from Carrel et al.13 as compiled in Balaton et al.9 If a gene had no call in any of the four AR datasets, but had a proportion of expression in Xi hybrids <0.22, we considered the gene “Subject”; genes with a greater proportion were called “Escape.”

To compare our calls with previous XCI consensus calls, we made the following modifications to the Balaton list: XG had been listed as a PAR gene, but we excluded it from our list of PAR genes because it is located at the PAR boundary and truncated on the Y chromosome. We updated its annotation to escape (“E”) since the Balaton table lists evidence for escape. The Balaton table lists XIST as “mostly subject” to XCI, but given its exclusive expression from Xi, we updated its status to escape (“E”). We manually examined all genes on our list that were not found in the Balaton list to make sure that genes were not misclassified due to differences in official gene names. For those genes still not present in the Balaton list after this correction, we list “No call”. To compare with our annotations, we grouped the Balaton calls into “Escape” if they were annotated as “PAR”, “Escape”, “Mostly escape”, “Variable Escape”, “Mostly Variable Escape”, or “Discordant”. We grouped Balaton calls into “Subject” if they were annotated as “Mostly subject” or “Subject”.

We compared our new calls with the Balaton calls for the 423 genes expressed in fibroblasts or LCLs, finding 48 where they differed. Of these, nine had a call in Balaton, but no call in the newer datasets. For two genes (TCEAL3, TMSB4X), there was no call in Balaton, but newer data enabled a call to be made (both “Subject”). Nineteen genes were called “Subject” in Balaton, but new data indicates that they have expression from Xi and we categorize them as “Escape.” The final eighteen genes were called “Escape” in Balaton, but new data suggested they have no expression from Xi. In total, our classification found 86 genes that “Escape”, 315 genes that are “Subject” to XCI, and 22 genes with “No call.”

Author contributions

Conceptualization, A.K.S.R. and D.C.P.; data curation, A.K.S.R. and H.S.; formal analysis, A.K.S.R., A.K.G., H.S., D.W.B., A.F.G., H.L.H., L.V.B., and J.F.H.; funding acquisition, A.K.S.R., A.F.G., L.V.B., S.M.D., N.R.T., M.M., and D.C.P.; investigation, A.K.S.R., S.P., L.B., A.D., and E.P.; methodology, A.K.S.R., A.K.G., H.S., A.F.G., and H.L.H.; project administration, A.K.S.R., J.F.H., L.B., A.B., P.K., N.B., P.C.L., C.K., and S.M.D.; resources, A.B., P.K., N.B., P.C.L., C.K., S.M.D., N.R.T., C.S.-S., and M.M.; software, A.K.S.R., A.K.G., H.S., A.F.G., H.L.H., and L.V.B.; supervision, A.K.S.R., J.F.H., N.R.T., C.S.-S., M.M., and D.C.P.; validation, A.K.G., H.S., D.W.B., and L.V.B.; visualization, A.K.S.R., A.K.G., and H.L.H.; writing – original draft preparation, A.K.S.R. and D.C.P.; writing – review and editing, A.K.S.R., A.K.G., D.W.B., A.F.G., H.L.H., L.V.B., J.F.H., S.M.D., and D.C.P.

Declaration of interests

The authors declare no competing interests.