We compared the performances of three existing gender-inferring methods and seGMM using the TGS data from the following 3 datasets: 1) Dataset 1: exon-targeted sequencing data of 1,000 genes (34 X chromosomal genes and two Y chromosomal genes) for a cohort of 110 males and 98 females from the 1,000 Genomes Project (Supplementary Table S1) (Genomes Project et al., 2010); 2) Dataset 2 (in-house): massive parallel sequencing of 785 deafness-related genes (eight genes in the X chromosome) for an in-house cohort of 8,805 males and 7,582 females; and 3) Dataset 3: targeted sequencing data of 187 autism risk genes (13 genes in the X chromosome) for a cohort of 42 females and 205 males from the Autism Clinical and Genetic Resources in China (ACGC) (Guo et al., 2018).

We also used the following two publicly available datasets (Supplementary Tables S2, S3) and one in-house dataset for analyzing the performance of seGMM in determining gender using WES and WGS data: 1) Dataset 4: exome sequencing data of 164 males and 118 females from the 1,000 Genomes Project (Genomes Project et al., 2015); 2) Dataset 5 (in-house): exome sequencing data of 1,257 males and 1,136 females; and 3) Dataset 6: high-coverage whole genome sequencing data of 11 males and 16 females from the 1,000 Genomes Project (Genomes Project et al., 2015).

The publicly available BAM files were previously mapped to the reference genome (GRCh37) and directly used for downstream analyses. For the in-house datasets, Fastp was used to remove the adapters and low-quality reads, and the quality of sequencing data was evaluated using measures such as Q20, sequence duplication levels, coverage, and GC content (Chen et al., 2018). Clean DNA sequencing reads were mapped to the human reference genome (GRCh37) using the BWA-MEM algorithm (Li and Durbin, 2009). Duplicated reads in the BAM files from the public and in-house datasets were removed using the sambamba tool (Tarasov et al., 2015). The variants were identified based on the Genome Analysis Toolkit best practices recommendations (McKenna et al., 2010) and filtered with VCFtools (Danecek et al., 2011) using parameters such as missing data in more than 50% of samples, minor allele count <3, overall SNP quality (QUAL) score <30, and read depth <5.

The model for seGMM included five gender-associated features, namely, X chromosome heterozygosity (XH), reads mapped to the X chromosome (Xmap), reads mapped to the Y chromosome (Ymap), the ratio of X/Y counts (XYratio), and the mean depth of the sex-determining region of the Y chromosome (SRY) gene (SRY_dep). The seGMM tool computed XH as the fraction of all genotypes on the X chromosome with two different allele calls, excluding the missing genotypes. Xmap/Ymap was computed as the fraction of high-quality reads (mapq > 30) that mapped to the X/Y chromosome divided by the total number of high-quality reads that mapped to the genome using the samtools algorithm (Li et al., 2009). XYratio was computed as the ratio of Xmap to Ymap (Xmap/Ymap). SRY_dep was determined using the mosdepth tool (Pedersen and Quinlan, 2018). The seGMM tool allows the users to customize feature selection for the GMM model because different TGS panel designs may only provide some features. For example, if the gene panel contains only genes located on the X chromosome, the relevant features on the X chromosome (XH and Xmap) are extracted and put into the model for gender determination.

The features extracted from the BAM and VCF files were normalized to the same level using the scale function in R 4.1.2 (R Core Team., 2021). Then, the mclust (v.5.4.9) R package was used to perform model-based clustering with the expectation-maximization (EM) algorithm and the samples were classified into two clusters (Scrucca et al., 2016). The gender was inferred based on the cluster results for a group of samples. The outliers were identified when uncertainty (probability of being assigned to two different clusters) was greater than 0.1. When a single sample was submitted, gender was inferred using the reference data that was analyzed with the same features as those in the seGMM model (Figure 1).

We defined the gates to classify individual karyotypes. The distribution of Xmap and Ymap in the females and males of the large cohort was normal. The ratio of samples with sex chromosomal abnormalities was 0.022% (Nielsen and Wohlert, 1990). This data was in agreement with the empirical rule, which states that 99.7% of normally distributed data lies within 3 standard deviations (sd) of the mean. Hence, we defined the normal gates as mean±3sd. The fold changes in Xmap or Ymap values indicated sex chromosomal aneuploidy.

To identify sex chromosomal abnormalities in the samples, we first calculated the mean value and standard deviation values of Xmap (mean_xmap and sd_xmap) and Ymap (mean_ymap, and sd_ymap) in the genetically determined male and female samples. The values for the males and females were denoted as m and f, respectively. The following six gates were then used to classify the karyotypes of individuals: •

XY Gate

The performance of seGMM was compared to PLINK 1.9, XYalign (v.1.1.6), and seXY (v.20170316). For PLINK 1.9, the pseudoautosomal region of the X chromosome was first split off with the parameter--split-x. Then, the parameter--check-sex was run without parameters. After reviewing the distribution of F estimates, the parameter--check-sex was rerun with parameters corresponding to the empirical gap. XYalign was performed following the method described in the original literature. The CHROM_STATS module was used to obtain the depths of chromosomes 1, X, and Y. The depths of X and Y chromosomes were normalized relative to the depth of chromosome 1. Then, a scatter plot of normalized X and Y chromosomes depth was plotted to assess gender in samples. Gender of the samples was inferred with seXY using the X.ped and Y.ped data that was derived from PLINK. The training dataset was provided by seXY. We expected to compare seGMM and other existing tools for all the six datasets. However, target gene panel data for datasets 2 and 3 did not contain genes on the Y chromosome. Therefore, the performance of XYalign and seXY was not available for these two datasets.

The STR analysis was performed using the customized multiplex PowerPlex® 16 System, which allowed co-amplification and four-color detection of amelogenin and other gene loci. The following primers were used for amplifying amelogenin: forward, 5′- GTT​AG​ACG​TGT​GCT​TCA​ACT​TCA​GCT​ATG​AGG​TAA​TTT​TTC—3′; reverse, 5′- ATC​CGA​CGG​TAG​TGT​CCA​ACC​ATC​AGA​GCT​TAA​ACT​GG-3′. All genetic loci were amplified simultaneously in a single tube and analyzed in a single lane. One of the primers for the amelogenin gene was labeled with carboxyrhodamine (ROX). The amplicons were separated in the ABI 3730XL Genetic Analyzer and the data was extracted using GeneMapper ID v3.2. The gender was inferred according to the peaks for the amelogenin gene. If only one peak was observed for the amelogenin locus, the gender was designated as female. If two distinct peaks differing by 6 bp were observed in the amelogenin locus, the gender was designated as male.

Genomic DNA (gDNA) was extracted using the MagMAX High Purity Free DNA Separation Kit (Magen, China). DNA concentration of the samples was measured using the NanoDrop One spectrophotometer (Thermo Fisher Scientific, United States). The working concentration of all DNA samples was 20 ng/µl. The primers targeting SRY, zinc finger protein Y-linked (ZFY), and deleted in azoospermia 1 (DAZ1) genes were designed using the Primer-BLAST online tool to determine the Y chromosome copy number (Ye et al., 2012). RPP30 was used as the internal control. All the primers used in this study are listed in Supplementary Table S4. The qPCR reaction mix included 0.6 µl of gDNA, 0.4 µl of each primer, 5 µl of iTaq™ Universal SYBR® Green Supermix (Bio–Rad, United States), and 3.6 µl of double-distilled water. Each sample was analyzed with three replicates. The quantitative real-time PCR assay (RT–qPCR) was performed in the QuantStudio 5 Real-Time PCR system (Thermo Fisher Scientific) using the following conditions: initial hot start cycle at 98°C for 2 min followed by 40 cycles consisting of denaturation at 98°C for 10 s, annealing at 60°C for 10 s, and the final extension step of 30 s at 72°C.

The distribution of XH, Xmap, Ymap, and XYratio values for dataset 1 (n = 208) shown in Figures 2A–D. The accuracy of seGMM was 99.52% and none of the samples were outliers (Figure 2E; Table 1). The accuracy of seGMM in females and males was 98.98 and 100%, respectively. The XYratio of one female sample (NA19054) resembled that of males and was incorrectly classified as male by seGMM. The gender-inference performance of seGMM for dataset 1 was superior to PLINK, seXY and XYalign. The PLINK tool analysis showed that the F coefficients for the dataset1 samples ranged from 0 to 0.9 and gap of F coefficients was not observed (Supplementary Figure S1). The accuracy of PLINK was 81.44% by running --check-sex without parameters. The accuracy of seXY for the dataset 1 was only 62.5% (Table 1). XYalign does not directly indicate predicted gender. Therefore, plotting the normalized sequence depth of the sex chromosomes and cluster samples along two ellipses using the stat_ellipse function resulted in a confidence level of 99.99%. XYalign plot showed that one female sample was located along with the male samples and three female samples were located between the two ellipses. Hence, the predicted gender of these four female samples was ambiguous (Supplementary Figure S2) and the accuracy of XYalign is 98.08%.

The performance of seGMM for the target gene panel data was validated using dataset 2 (n = 16,387) and dataset 3 (n = 247). The read counts, base quality, and GC distribution of the sequencing data of all the 16,387 subjects in dataset 2 was assessed. The average total number of sequence read per sample was 10.72 million. The average quality score for all bases was above 30 and the average GC content was 50.11% per subject. The average targeted sequence coverage was 90.43%, and unique mapping rate of each sample was 99.26%. We identified 16,988 variants in eight genes located on the X chromosome. Because the target gene panel for dataset 2 did not contain genes on the Y chromosome, we only used XH and Xmap to analyze the performance of the seGMM model. XH and Xmap plots showed distinct clusters for males and females (Supplementary Figure S3). The overall accuracy of seGMM and PLINK was 99.92 and 87.10%, respectively (Figure 2F; Table 2). The accuracy of seGMM in females and males was both 99.98%. One self-reported female sample (HL-001200) and two self-reported male samples (CTRL-002692 and CTRL-002753) were misclassified. Therefore, we performed STR analysis with the sex chromosome marker gene, amelogenin, to verify the gender of these three ambiguous samples. All these three samples were identified as males because two distinct peaks were observed with a difference of 6 bp for the amelogenin gene (Figures 3A–C). CTRL-002692 and CTRL-002753 were misclassified because all other male samples had a XH value of 0, while these two samples had a non-zero value (XH_CTRL-002692 = 0.0025; XH_CTRL-002753 = 0.0046), which could be caused by the individual variation in targeting region.

The overall accuracy of seGMM for dataset 3 was 92.31% (97.56% for females and 84.88% for males, Figure 2G and Table 2). The accuracy of PLINK was only 38.87% for dataset 3. The performances of seGMM and PLINK were significantly better for datasets 1 and 2 compared to dataset 3 because the number of X chromosome SNPs (81) were lower and sequencing data for the Y chromosome was absent in dataset 3, thereby affecting the distribution of XH values from the male and female samples (Supplementary Figure S4A). In contrast to PLINK, seGMM collected additional information for the reads mapped to the X chromosome, thereby enabling better separation between the female and male samples (Supplementary Figure S4B). Furthermore, we assessed the performance of seGMM using features only extracted from the X chromosome in dataset 1. The seGMM tool showed that 59 samples were outliers and the accuracy for the remaining samples was only 84.56%. We then evaluated the computation time of different methods using 1 core, 10 cores and 20 cores on a server with 64 Intel(R) Xeon(R) CPU E7-8895 v3 at 2.60 GHz. The analysis time for the seGMM tool was longer than PLINK and seXY because it collected additional features such as reads mapped to the X and Y chromosomes. Moreover, the analysis time for seGMM with 1 core was longer than XYalign and 10 times faster than XYalign with 20 cores (Supplementary Table S5).

We then evaluated the performance of the seGMM tool for the WES and WGS data. First, we analyzed the publicly available WES data (dataset 4). The accuracy of seGMM was 100% for the samples in dataset 4 (n = 282) (Table 3 and Supplementary Figure S5). The accuracy of PLINK and seXY was also 100%. The accuracy of XYalign was 99.65% (Supplementary Figure S6).

Next, we analyzed the in-house WES data (dataset 5, n = 2,393) using seGMM and other tools. In dataset 5, the average number of sequencing reads per sample was 114.43 million. The average Q20, Q30 and GC content of the reads per subject was 97.36%, 93.21%, and 51.23%, respectively. Furthermore, the average unique mapping rate for each individual sample was 99.92%. We identified 89,273 variants on the X chromosome and 4,866 variants on the Y chromosome. The concordance between inferred gender and self-reported gender using the seGMM tool based on the five features for the in-house WES dataset 5 was 99.75% (99.76% for males and 99.74% for females, Figure 4A, Supplementary Figure S7). Six mismatched samples (HL-005584, HL-006009, HL-006904, HL-007335, HL-007935 and HL-012246) were identified by comparing SNP-inferred gender and self-reported gender. This indicated misregistration of clinical information for some samples. Therefore, we performed STR analysis to validate the gender of these six samples. Three samples were classified as females because they showed only one peak for the amelogenin locus, whereas the remaining three samples showed two distinct peaks with a difference of 6 bp and were classified as males (Table 4). The results demonstrated that the actual accuracy of seGMM prediction was 100%. We also evaluated the correlation between age and reads mapped to the Y chromosome in the male samples (Supplementary Figure S8A) and reads mapped to the X chromosome in the female samples from the in-house WES dataset 5 (Supplementary Figure S8B). The results showed significant negative correlation (p = 6.899e-09; correlation coefficient: −0.17) between reads mapped to the Y chromosome and age, thereby indicating loss of Y chromosome during aging.

We then compared the performances of PLINK, seXY and XYalign for dataset 5 using the corrected gender information. The accuracy of PLINK was 99.79% with five mismatched samples (HL-033182, HL-020292, HL-011500, HL-019211 and HL-012554) (Table 3). The accuracy of XYalign was 99.91% with two mismatched samples (HL-009389 and HL-012554). Overall, six samples were mismatched, as predicted by PLINK and XYalign. STR analysis showed that the gender of these samples was consistent with their self-reported gender and matched the predicted results of seGMM analysis (Table 4). Furthermore, the accuracy of seXY was 49.23%. The loss of accuracy in seXY for dataset 5 was because the distribution of Y chromosome missingness in the male and female samples were confounded (Supplementary Figure S9). Finally, the performances of these tools were assessed using the WGS data (dataset 6, n = 27). The accuracy of all tools was 100% for dataset 6 (Table 3).

The seGMM tool can identify six sex chromosomal karyotypes (XX, XY, XYY, XXY, XXX, and X) using Xmap and Ymap. In a large cohort, the distribution of Xmap and Ymap was normal in females and males. The Xmap or Ymap values of samples with sex chromosome abnormalities such as XYY and XXY were significantly different and were recognized as outliers compared to samples with normal sex chromosomes. Three samples in dataset 5 (HL-029620, HL-009382 and HL-019110) were classified as the XYY karyotype. In dataset 5, the average rate of reads mapping to the X and Y chromosomes in the female samples were 0.035 ± 0.0020 and 1.61e-05 ± 3.36e-05, respectively, and 0.018 ± 0.0011 and 0.00067 ± 0.00014, respectively, for the male samples. The rate of reads mapping to the Y chromosome for the three outlier samples was twice as high as the mean value of Ymap in all the male samples (Ymap_HL-029620 = 0.0015, Ymap_HL-009382 = 0.0018, and Ymap_HL-019110 = 0.0015), thereby suggesting a XYY karyotype by seGMM (Figure 4B). Furthermore, although HL-009389 and HL-012554 samples were located close together in the middle of the plot, they were correctly predicted by seGMM as female and male, respectively (Figure 4B). This is because features such as Xmap and SRY_dep, which are not shown in Figure 4B, clearly separated all female and male samples (Supplementary Figure S10). This demonstrated the significance of incorporating key features to improve the accuracy of the gender prediction model. In the other datasets, sex chromosomal abnormalities were not identified.

Next, we evaluated the accuracy of the data-based sex chromosome karyotype of these three samples by analyzing the copy number ratios of Y chromosome-specific genes (SRY, ZFY and DAZ1) by RT–qPCR. We used HL-007935 and HL-012246 samples as controls for females and males based on the STR analysis results. The copy number ratio for normal females was 0. The copy number ratio for normal male samples was 1. We analyzed the copy number ratios of HL-029620, HL-009382 and HL-019110 samples in dataset 5 and found that the copy number ratio of HL-019110 was 2 (Figure 5). This confirmed that the karyotype for the HL-019110 sample was XYY.