http://orcid.org/0000-0002-0038-461X
Figures
Abstract
Over the last decade, a large number of nucleotide sequences have been generated by next-generation sequencing technologies and deposited to public databases. However, most of these datasets do not specify the sex of individuals sampled because researchers typically ignore or hide this information. Male and female genomes in many species have distinctive sex chromosomes, XX/XY and ZW/ZZ, and expression levels of many sex-related genes differ between the sexes. Herein, we describe how to develop sex marker sequences from syntenic regions of sex chromosomes and use them to quickly identify the sex of individuals being analyzed. Array-based technologies routinely use either known sex markers or the B-allele frequency of X or Z chromosomes to deduce the sex of an individual. The same strategy has been used with whole-exome/genome sequence data; however, all reads must be aligned onto a reference genome to determine the B-allele frequency of the X or Z chromosomes. SEXCMD is a pipeline that can extract sex marker sequences from reference sex chromosomes and rapidly identify the sex of individuals from whole-exome/genome and RNA sequencing after training with a known dataset through a simple machine learning approach. The pipeline counts total numbers of hits from sex-specific marker sequences and identifies the sex of the individuals sampled based on the fact that XX/ZZ samples do not have Y or W chromosome hits. We have successfully validated our pipeline with mammalian (Homo sapiens; XY) and avian (Gallus gallus; ZW) genomes. Typical calculation time when applying SEXCMD to human whole-exome or RNA sequencing datasets is a few minutes, and analyzing human whole-genome datasets takes about 10 minutes. Another important application of SEXCMD is as a quality control measure to avoid mixing samples before bioinformatics analysis. SEXCMD comprises simple Python and R scripts and is freely available at https://github.com/lovemun/SEXCMD.
Citation: Jeong S, Kim J, Park W, Jeon H, Kim N (2017) SEXCMD: Development and validation of sex marker sequences for whole-exome/genome and RNA sequencing. PLoS ONE 12(9): e0184087. https://doi.org/10.1371/journal.pone.0184087
Editor: Daniel Doucet, Natural Resources Canada, CANADA
Received: May 15, 2017; Accepted: August 17, 2017; Published: September 8, 2017
Copyright: © 2017 Jeong et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the paper and its Supporting Information files.
Funding: This work was supported by grants from the National Research Foundation of Korea, Grant NRF-2011-0030049, NRF-2014M3C9A3064552 (http://www.nrf.re.kr/) and the Korea Research Institute of Bioscience and Biotechnology, NK., KRIBB Initiative program (http://www.kribb.re.kr/). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Next-generation sequencing (NGS) has jump-started a wide array of novel genomics research in the last decade. NGS is replacing two previous technologies—Sanger sequencing and array-based methods. Now, researchers can directly examine genomes, transcriptomes, and epigenomes through whole-exome/genome sequencing, RNA-Seq, DNA methylation, ChIP-Seq, and other methods. In a population-scale experiment, researchers may create hundreds of NGS datasets in a single project, and public databases such as 1000 Genome project [1], Sequence Read Archive (SRA) [2], and The Cancer Genome Atlas (TCGA) [3] have already archived petabytes of information from genome and transcriptome sequencing. However, the sex of individuals analyzed to create these datasets are not always explicitly denoted in the associated reports. Males and females in species with sex chromosomes can have two different combinations of sex chromosomes, XX/XY for mammals and ZW/ZZ for avian, which could lead to different gene expression and DNA methylation patterns. We believe some datasets are assigned the incorrect sex as a result of errors in annotation. Finally, large sequencing centers have sex markers for array-based technology, but not for NGS. If the sex of individuals in a given NGS dataset was known, as part of quality control before bioinformatics analysis, we could prevent errors, including label error, caused by analysis of mixed samples.
Sex-specific genetic markers can be obtained using a simple molecular identification method based on polymerase chain reaction (PCR) [4]. Chaveerach et al. analyzed amplification products from PCR targeting four loci to deduce the sex of individuals sampled [5]. These two methods require primer sequence information a priori for a known polymorphic sex marker. Array-based technologies typically use single nucleotide polymorphism (SNP) markers from the human amelogenin genes, AMELX and AMELY, because they are single copy genes and can be used to extract sex-specific sequence markers [6]. Some software packages, such as PLINK [7], PLATO [4], seXY [8], and Golden Helix (Bozeman, MT, USA, http://www.goldenhelix.com), use whole X chromosome heterozygosity or intensities to identify sex. These analyses are based on genotypes, intensities, and heterozygosity thresholds of the X chromosome. Qu et al. took into account both intensities and genotypes of SNPs on the X chromosome simultaneously and calculated the probability of errors in sex identification by logistic regression [9]. The above methods all require SNP polymorphic markers and have biases from copy number variations or large-scale structural variations in individual genomes.
With many pipelines, it is necessary to align NGS reads onto a reference genome in order to identify the genotypic sex of an individual, XX/XY or ZZ/ZW. This process can take hours or days, depending on dataset size, and consume hundreds of gigabytes of disk space. The sex of a given individual can typically be deduced by measuring the B-allele frequency of the X or Z chromosome. Instead of using all NGS reads and the whole reference genome, we developed a simple strategy to use tens of sex-specific marker sequences from syntenic regions of the sex chromosomes. Alignment software chooses whether a given read can be mapped onto one of the sex-specific markers. If the alignment score is not high enough, they are reported as un-mapped. We can thereby collect mapped reads and count the total numbers of hits for each sex-specific marker sequence.
We have implemented our analysis pipeline as (1) designing sex-specific marker sequences, (2) training using a known dataset, and (3) optimal sex marker sequence selection. Because the algorithm uses differences within the sex chromosomes, it can be applied to genome and transcriptome sequencing. Small datasets, such as whole-exome and RNA sequencing, take a few minutes to analyze, while whole-genome sequencing takes more than 10 minutes, depending on computing platform. We have validated our pipeline for two organisms, humans (XX/XY) and chickens (Gallus gallus, ZZ/ZW). SEXCMD is freely available at https://github.com/lovemun/SEXCMD.
Materials and methods
Extracting sex-specific marker sequences from a reference genome
Sex-specific marker sequences were extracted through six simple steps. First, sex chromosome sequences were aligned with each other using LASTZ [10]. Second, we found a syntenic region between sex chromosomes that did not have long blocks of exactly matching sequences. Blocks without polymorphisms should be less than 35 base pairs (Python script: 2.sex_marker.py)—no longer than read length. This internal parameter can be adjusted if no candidate syntenic regions can be found because the sex chromosome reference genome sequences are incomplete. We used human and chicken reference genomes in this study, and were able to identify appropriate syntenic regions in their sex chromosomes. However, the total number of sex-specific polymorphic regions in the pig genome (susScr3) were not sufficient. Considering recent improvements in long read lengths, i.e., up to 150 bp, minimum length of marker sequences was set to 151 bp with 5 bp mismatches (Python script: 3.sex_marker_filtered.py). This length should be longer than the input NGS read lengths. Finally, we used BLAST [11, 12] to check whether those marker sequences were unique by comparing against all autosomal sequences with minimum identity of 90%. Most of the candidate sequences were removed in this filtering step. If the total number of marker pairs is less than 10, two parameters can be adjusted—–“minimum exact match blocks” in 2.sex_marker.py and “minimum length of marker sequences” in 3.sex_marker_filtered.py. If more than 30 marker sequences are identified, only the sequences within genic regions can be selected to improve efficiency when analyzing data from RNA sequencing. We can assume that most syntenic regions will be within genic regions of sex chromosomes. Non-genic regions can be excluded from the analysis if there are no hits from RNA sequencing alignments. FASTA sequences from final candidate sex-specific marker sequences were indexed by BWA [13] with the “–a is” option. Note that we used the “is” mode because the size of the FASTA file is small.
The workflow for extracting sex-specific marker sequences is depicted in Fig 1. The main goal of this pipeline is to find polymorphic syntenic regions with minimal exact match blocks. Once the sequences of the sex chromosomes are known, we can design sex-specific marker sequences.
Two sex chromosomes were aligned with each other using LASTZ, and syntenic regions with polymorphisms were extracted. Final sex-specific marker sequences were selected after removal of similar sequences (90% identity by BLAST).
Software and URLs are described on the SEXCMD webpage along with examples of command lines.
Mapping and counting sex-specific hits and interpretation
Each sex-specific marker pair consists of two sex-specific marker sequences originated from the two sex chromosomes. The main concept is to align small numbers of input reads onto those marker sequences and count the number of hits for each pair. Individuals with heterozygous sex chromosomes will have equal numbers of hits on both sex-specific marker sequences in each pair, whereas homozygous individuals will only have hits for one marker sequence per pair (markers on the X or Z chromosome). Because we use only small target sequences instead of alignment with the whole genome, calculation time is reduced. It typically takes hours or days to align all whole-exome/genome or RNA sequencing reads onto a reference genome, but less than tens of minutes to align onto sex-specific marker sequences. We used the BWA-MEM algorithm to align reads onto sex-specific marker sequences and SAMTOOLS [14] to count hits for each marker pair.
Individuals with XY or ZW chromosomes will have a ratio of 1:1 hits for each sex chromosome, and individuals with XX or ZZ chromosomes will have a ratio of 1:0 hits for the two sex chromosomes. Essentially, we used two numbers–the number of hits on the dominant chromosome (Y or W) and the number of hits on the recessive chromosome (X or Z)–to identify the sex of a given individual. Owing to false hits from wrong alignments, we set the minimum fraction of Y or W marker hits on a homozygous individual to be 0.2. If the fraction of hits on the dominant chromosome (Y or W) markers was less than 0.2, we considered the individual to be heterozygous.
The proportion of sex chromosome reads in whole genome sequencing is frequently smaller than in whole exome and RNA sequencing, because whole exome and RNA sequencing are focused on the coding regions of a genome. It takes tens of minutes to align entire datasets onto sex-specific marker sequences. In order to make calculation time shorter, we implemented a ‘test mode’ to measure the minimum number of input reads. It sums the number of hits on each marker and suggests a minimum number of input reads to be analyzed. We recommend that the sum of read counts for all markers be greater than 100 hits. Once we set the minimum requirements, the program uses that number of reads instead of analyzing the entire dataset. This significantly reduces calculation time.
To estimate the minimum number of input reads, we selected a hundred samples for each sequencing type (whole-exome, genome, and RNA sequencing) and measured the average mapped read counts according to input number of reads. The minimum number of input sequence reads necessary for identifying the sex of an individual differs depending on the sequencing type (Fig 2). For human data, we set default values of input read sizes to be five million for whole exome or RNA sequencing datasets and one hundred million for whole genome sequencing datasets.
Red arrows indicate minimum read counts: 5 million (5×106) reads for whole-exome sequencing and RNA sequencing and 100 million (1×108) reads for whole-genome sequencing. The red horizontal line denotes the minimum average read counts of sex-specific marker sequences.
SEXCMD is a simple and fast method of identifying the sex of given individuals without mapping whole NGS datasets onto a reference genome.
Results
Input data
A reference genome with two sex chromosomes is essential for the identification of sex-specific marker sequences. We used human (hg38) and chicken (galGal4) reference genomes from UCSC genome browser (http://genome.ucsc.edu). We downloaded 400 human whole exome sequencing datasets (182 males and 218 females) from the 1000 Genomes Project (http://www.1000genomes.org/data/) and 48 human whole genome sequencing datasets (27 males and 21 females) from the Sequence Read Archive (SRA, http://ncbi.nlm.nih.nih.gov/sra/). In addition, we used 202 human whole exome sequencing datasets and 378 human whole genome sequencing datasets from our in-house database. A total of 131 human RNA sequencing datasets were downloaded (59 males and 72 females) from the ArrayExpress archive (http://www.ebi.ac.uk/arrayexpress/). For chicken, we used 120 whole genome sequencing datasets and 36 chicken RNA sequencing datasets from our in-house database. A summary of these datasets including the sex of individuals sampled is shown in Table 1 and the accession numbers in public databases are provided in S1 Table.
Sex-specific marker sequences
Sex-specific marker sequences were derived using Python scripts as described in the Methods section. Genomic coordinates of human and chicken marker sequences are summarized in Tables 2 and 3. In the human genome, SEXCMD selected AMELX, USP9X, NLGN4X, TBL1X, and KDM6A from the X chromosome, and AMELY, USP9Y, UTY, NLGN4Y, TBL1Y, and KDM5C from the Y chromosome. In the chicken genome, SEXCMD selected SMAD2 from the Z chromosome, and FET1 and UBAP2 from the W chromosome. The average length and number of mismatches for marker sequences were 216 bp and 28 in humans and 196 bp and 31 in chickens, respectively.
Interpretation and identification of sex
We used Python and R scripts to implement SEXCMD. These scripts can be applied to genome and transcriptome sequencing datasets that contain sex-specific fragments. We measured the accuracy of sex identification for humans and chickens using the datasets we downloaded and achieved 100% for all three data types–whole exome, whole genome, and RNA sequencing (Table 4).
To obtain accurate results from SEXCMD, we recommend the number of input read sequences should be at least five million. The more input reads used, the better the accuracy of the results. However, calculation time must also be considered. Current implementation of SEXCMD takes approximately 10 minutes for whole genome sequencing, 5 minutes for whole exome sequencing and 1 minute for RNA sequencing datasets.
Implementation
SEXCMD is available at https://github.com/lovemun/SEXCMD and consists of R and Python scripts. We used a Linux environment for development and validation. Python, R, BWA, and SAMTOOLS should be installed and their PATH included as an environmental variable.
Conclusions
SEXCMD is a novel pipeline to design sex-specific marker sequences from a reference genome and to identify the sex of an individual analyzed in a given NGS dataset. This pipeline successfully extracts the polymorphic syntenic region of two sex chromosomes, and the NGS dataset of interest will be aligned onto these small sex-specific marker sequences. Then, the number of hits on each sequence is counted. With a few simple criteria, we have achieved 100% accuracies for human and chicken whole exome, whole genome, and RNA sequencing datasets. Perhaps most importantly, the calculation time is typically less than ten minutes because the pipeline only uses a fraction of the input dataset–just enough to identify sex.
This pipeline can be very useful if there is no sex information provided with the NGS reads. Most large genome projects do not give the sex of individuals analyzed because one may think it is trivial information. Even when the sex of individuals is provided, we have found that there is some level of inaccuracy in public databases because of incorrect annotations or mixed samples. Large genome sequencing facilities could use the SEXCMD pipeline prior to bioinformatic analysis to eliminate the possibility of errors due to incorrect annotation or mixed samples.
One limitation of SEXCMD is that it depends on the availability of high quality reference genomes with two sex chromosomes. Many well-studied species have a high quality reference genome, but this is not always the case. SEXCMD may not be able to extract sex-specific marker sequences if the differences between sex chromosomes are not that big enough. If you have a bunch of scaffolds, not in chromosomal level assembly, and if you can decide which the scaffolds are on sex chromosomes of X and Y or Z and W chromosomes, then you can use SEXCMD. One can adjust some parameters in SEXCMD to design proper sex-specific marker sequences. This does not consume CPU time or storage space because SEXCMD is not saving any temporary files.
Most software packages measure the B-allele frequency of X and Z chromosomes and align all input reads onto the reference genome prior to identifying the sex of an individual. For the simple purpose of sex identification, our SEXCMD pipeline is unique.
Finally, we can assume that genotypes of X or Z chromosomes will change if one aligns NGS reads onto heterozygous sex chromosomes regardless of the sex of the given individual. Because the human Y chromosome is inherited from the male parent and is very diverse, researchers tend to ignore this and simply use the XY reference genome for all datasets. This can lead to incorrect genotypes for syntenic regions of sex chromosomes.
We have upload all datasets to SEXCMD github page. Supporting information provides accession number for each sequencing type from 1000 genome project, Sequence Read Archive, and ArrayExpress archive. Because some dataset is ongoing project, the uploaded data contains only mapped sequences to sex markers.
Supporting information
S1 Table. Accession numbers of public database in this analysis.
https://doi.org/10.1371/journal.pone.0184087.s001
(XLSX)
Acknowledgments
This work was supported by grants from the National Research Foundation of Korea (NRF-2011-0030049, NRF-2014M3C9A3064552) and the KRIBB Initiative program.
References
- 1. The 1000 Genomes Project Consortium (2015). "A global reference for human genetic variation." Nature 526(7571): 68–74. pmid:26432245
- 2. Leinonen R., Sugawara H., and Shumway M. (2011). "The sequence read archive." Nucleic Acids Res 39(Database issue): D19–21. pmid:21062823
- 3. Cancer Genome Atlas Research, N. (2008). "Comprehensive genomic characterization defines human glioblastoma genes and core pathways." Nature 455(7216): 1061–1068. pmid:18772890
- 4. Grady B. J., Torstenson E, Dudek S.M., Giles J., Sexton D., and Ritchie M.D. (2010). "Finding unique filter sets in PLATO: a precursor to efficient interaction analysis in GWAS data." Pac Symp Biocomput: 315–326. pmid:19908384
- 5. Arunrat Chaveerach R. S., Tanee Tawatchai, Sanubol Arisa, Thooptianrat Tikumporn, Faijaidee Waraporn, and Yaipool Kittibodee (2015). "Genetic markers for sex identification in Thai population." Forensic Science International: Genetics Supplement Series 5: e390–e391.
- 6. Tzvetkov M. V., Meineke I., Sehrt D., Vormfelde S.V., and Brockmoller J. (2010). "Amelogenin-based sex identification as a strategy to control the identity of DNA samples in genetic association studies." Pharmacogenomics 11(3): 449–457. pmid:20235797
- 7. Purcell S., Neale B., Todd-Brown K., Thomas L., Ferreira M.A., Bender D., et al. (2007). "PLINK: a tool set for whole-genome association and population-based linkage analyses." Am J Hum Genet 81(3): 559–575. pmid:17701901
- 8. Qian David C., Busam Jonathan A., Xiao Xiangjun, O’Mara Tracy A., Eeles Rosalind A., Schumacher Frederick R., Phelan Catherine M., Amos Christopher I.; seXY: a tool for sex inference from genotype arrays. Bioinformatics 2017; 33 (4): 561–563. pmid:28035028
- 9. Qu C., Schuetz J.M., Min J.E., Leach S., Daley D., Spinelli J.J., et al. (2011). "Cost-effective prediction of gender-labeling errors and estimation of gender-labeling error rates in candidate-gene association studies." Front Genet 2: 31. pmid:22303327
- 10.
Harris R. S. (2007). Improved pairwise alignment of genomic DNA, The Pennsylvania State University.
- 11. Boratyn G. M., Schaffer A.A., Agarwala R, Altschul S.F., Lipman D.J., and Madden T.L. (2012). "Domain enhanced lookup time accelerated BLAST." Biol Direct 7: 12. pmid:22510480
- 12. Altschul S. F., Madden T.L., Schaffer A.A., Zhang J., Zhang Z., Miller W., and Lipman D.J. (1997). "Gapped BLAST and PSI-BLAST: a new generation of protein database search programs." Nucleic Acids Res 25(17): 3389–3402. pmid:9254694
- 13.
Li, H. (2013). "Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM." arXiv 1303.3997v1.
- 14. Li H., Handsaker B., Wysoker A;, Fennell T., Ruan J., Homer N., et al. (2009). "The Sequence Alignment/Map format and SAMtools." Bioinformatics 25(16): 2078–2079. pmid:19505943