International Journal of

1. Introduction

Infertility is an important social health condition, and is defined as the inability to conceive naturally after at least one year of unprotected intercourse. Approximately 15% couples suffer from infertility and 30%-35% of those cases relate to male factor infertility (1). Numerous factors (ranging from lifestyle to heredity) may affect semen parameters and male infertility (2). Approximately 15% of male infertility cases are attributable to genetic disorders (3). The majority of mechanisms involved in disease heritability are associated with chromosomal abnormalities and genetic variations (4). According to OMIM, autosomal or sex chromosome has more than 60 genes, leading to male infertility. CFTR mutation has been implicated in 60%-90% of cases involving infertile men, with congenital bilateral absence of the vas deferens or epididymal obstruction. Mutations in X-linked genes, such as USP26 and SOX3, have been reported in cases of severe spermatogenesis impairment (5, 6).
Approximately 30%-40% infertility cases in men with unknown etiology may be associated with genetic factors (6). The high proportion of idiopathic infertility in males is attributed to the limitations of traditional gene identification techniques, such as the Sanger sequencing. The most notable disadvantages of these approaches include length of time, low throughput data, and a broad region for candidate genes. Most of these limitations have been addressed by next-generation sequencing (NGS) (7).
The NGS approach, or whole-exome sequencing (WES), is a powerful and unbiased tool for identifying genetic variation by capturing genome coding regions (8). WES covers a region of approximately 1%-1.5% of the human genome, where approximately 85% of causative mutations are located. Overlap-, de novo-, extreme phenotypes-, and familial-based are the four main strategies for WES analysis. Of these, family-based WES is an efficient method of identifying potential causal variants (9). This approach led to the identification of causal variants that traditional methods had failed to detect in several pedigrees (10). The WES specialist can generate large amounts of data and, consequently, be of benefit to clinics. However, excluding the called variations across the whole-exome remains an issue, and it is not reasonable to verify all types of candidates via the Sanger sequencing (11, 12). Indeed, prioritizing the best candidate pathogenic variants is the main challenge. Many previous studies have focused on the family-based strategy to increase the efficiency of discovering candidate genes and reducing the number of private variants (13).
The purpose of this study was to illustrate the basic familial framework for the purpose of identifying genes that relate to male infertility. To this end, we studied seven cases in two families. In family I, we sequenced two affected members and their parents to filter out inherited variants. In family II, we focused on the most distantly related family members to reduce shared benign variations. The advantage of this approach, by sequencing non-affected siblings, is that not many private variations were missed.
 
2. Case presentation
 
2.1 Human subjects and DNA samples
In this case series study, we examined seven cases from two unrelated families (Figure 1). DNA was extracted from the peripheral blood leukocyte using salting out methods.. The quality of the DNA was checked using an agarose gel and NanoDrop analysis.
 

1. 2. Whole-exome sequencing platform

The library of all seven samples was prepared using the SureSelectXT Library Preparation Kit (Agilent Technologies, Santa Clara, CA, USA). As a next step, cluster generation, paired-end sequencing was applied on the IlluminaHiSeq 2000 platform using TruSeqv3. The binary base call (BCL) was converted to FASTQ using the Illumina bcl2fastq package.
 

1. 3. Quality evaluation of the raw data

The preprocessing and generation of raw reads (FASTQ) involved 3´ end adaptor clipping, primer removal, and the trimming of poor base sequence quality. The sequenced data were assessed using the FastQC tool (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) to analyze Phred score distribution, together with reads, GC content distribution, read length distribution, and sequence duplication level. The adaptor was removed using Trimomatic (version 1.04.636; www.github.com/ExpressionAnalysis/ea-utils/blob/wiki/FastqMcf.md).
 

1. 4. Alignment and duplicated PCR removed

Following the assessment of quality control, the raw reads were aligned to an established human reference genome (version hg19UCSC indexed in the FASTA format (14). Millions of reference genome-scale short reads have been globally streamlined using the Burrows-Wheeler Aligner (BWA) with MEM algorithm alignment in SAM format (15). The default parameters in the software were based on the Burrows-Wheeler conversion. The aligned reads were converted to a BAM file using SAMtools software (16), which was able to clean up and flag read-pairing information. The BAM file was then sorted by genomic location and indexed by SAMtools to save space and help the subsequent process. The Picard tool was combined with bamtools to filter out the mismatching and inappropriate reads for the purpose of assembled genomic data. The Picardtool identified duplicates from the PCR library. Data distribution and reads coverage (alignment statistics) were then evaluated with the CalculateHsMetrics package. Base quality was recalibrated, using GATK covariance recalibration (version 2.8; Broad Institute, Cambridge, MA, USA; www.broadinstitute.org/gatk/).
 

1. 5. Variant calling and annotation

Variant identification is the key stage in NGS data analysis. In this step, two programs, GATK and SAMtools, were applied to create a VCF file containing all the sites with potential variants (16). VCF files were filtered based on two criteria: depth of coverage and Fred score quality (DP > 9 -'QUAL > 25). After calling variants list,, the ANNOVAR software tool was used to annotate screened variations and connect the three annotation modes according to the type of gene, region, and filter (17).
 

1. 6. Variant filtering

The number of candidate variants was reduced through four steps. The first step was to filter the variants out of the coding or splicing region. The second step was to prioritize the variants with relatively low minor allele frequency (MAF <1%), or novel variants. The third step was to exclude synonymous variants and retain splice, nonsense, Indel, and nonsynonymous mutations. To reduce the number of candidate genes, we focused on the identification of genes related to two strategies (homozygosity and linkage-based). Family I underwent homozygosity strategy because the DNA samples of both parents were available. Family II was selected for linkage-based strategy because DNA samples from all affected cases were available. Five infertile men in this family were diagnosed with asthenozoospermia (ASZ); the DNA samples of two brothers and one of the most distant affected relatives were sequenced to identify the shared variations. In addition, the homozygosity region for each family was obtained using exome data. A homozygous variant in the homozygous region was prioritized.
 

1. 7. Variant analysis

Specific variants were assessed for evolutionary conservation using the GREP and PhastCons tools. Variants with a GREP score of >2.0 and a PhastCons score of >0.3 were considered a conserved variant. The prediction of candidate variants that affected protein or phenotype function was performed using five tools: Mutation Taster, PolyPhen, Sift, CADD, and Proven (18). The variants predicted to be disease-causing by more than two of these tools were further analyzed. A deleterious effect of variants on protein structure was evaluated by the online web service HOPE. Pathogenic variation is able to disrupt intron-exon splice sites, exonic splicing enhancers (ESE) and exonic splicing silencers (ESS), which can affect gene expression and cause genetic diseases by aberrant pre-mRNA splicing. The annotation of intronic and exonic mutations, leading to splicing defects, was performed using the ESEfinder and Human Splicing Finder web resource tools (19, 20).

 2.8. Ethical consideration
As a first step, written informed consent for genetic studies was obtained from subjects. This case series study was conducted in accordance with the protocols approved by the Institutional Review Board of the Royan Institute Research Center and the Royan Ethics Committee, Tehran, Iran.

 2.9. Bioinformatics tools

• FastQC tool
• Trimomatic
• Burrows-Wheeler Aligner
• SAMtools
• BamTools
• Picard tool
• CalculateHsMetrics package
• GATK
• ANNOVAR

 .3Results
3.1. WES Run Statistics and alignment results
Approximately 2.5 gigabytes of exome sequence were obtained for each of the seven samples. All sequence sets yielded a comparable total number of bases, reads, GC (%), Q20 (%), and Q30 (%). Table I shows the raw stats data of each sample. FastQC was performed to check the quality of the raw data, focusing on base quality score, GC content, N content, and sequence duplication level. Mapped results, previously generated by BWA, were computed to create the number of total reads, paired reads, properly-paired reads, duplicate reads, duplication rate percentage (number of duplicates/total reads), number of reads mapped/aligned, and unique reads mapped (Table II).

 3.2. Variant detection
Variants were called using two different tools, GATK (version 2.6) and SAMtools (version 0.1.18), which are the most widely used, and the Unified Genotyper algorithm was applied. SSCS and MSCS were separately performed for family I and family II. The number of SNPs and Indels was determined by different callers, and we evaluated the effects of SSCS and MSCS on the variant for each family (Table III). Through SAMtool, SSCS increased the number of raw SNPs by 1.6% and 1.1% in family I and family II, respectively. In contrast, MSCS increased the number of raw SNPs called by 1.2% with GATK in both pedigrees. Variants were filtered based on the depth of coverage and quality. In SSCS and MSCS, variant filtering removed 15% and 21% of raw variants called by SAMtools, respectively, whereas 16% and 9.2% of those called by GATK were filtered. These results reflect the importance of filtering in MSCS by SAMtools. We estimated the Ts/Tv ratio of SNP sets using two SSCS and MSCS pipelines. As shown in Figure 2a, in SSCS with SAMtools, the raw SNP sets had a Ts/Tv ratio of between 1.72 and 1.88, while the Ts/Tv ratio increased in the filtered data sets. The Ts/Tv ratio was constant in both unfiltered and filtered MSCS data sets with SAMtools. The Ts/Tv ratio of the raw SNP sets was 1.34-2.07 when the GATK was used, and increased in both SSCS and MSCS with GATK (Figure 2b). The results demonstrate that 83%-91% of SNVs are called by both GATK and SAMtools (Figure 3).

 3.3. Variant filtering
In this stage, the called variants were filtered in four stages. In the first stage, we excluded approximately two-thirds of variants, which were located in the intronic and intergenic regions. Less than 50% of variants were likely to be pathogenic. Of these, new and rare variants with MAF <1% required further analysis (Table IV). As the affected cases in family II were from a consanguineous marriage, we focused on homozygote variants. The number of private variants in pedigree II is shown in Table V. Approximately 0.05% (SAMtools) and 0.03% (GATK) of the called variations remained after initial filtering in family I (Table VI).
 







 3. Discussion
NGS-based technology is a powerful tool for identifying the genetic basis of the human phenotype. As one of the leading molecular techniques in the field of reproductive medicine, WES has a major impact on our understanding of the genetic causes of male infertility. WES generates a large amount of genetic data; however, the processing of WES data can be complex (11). As a standard pipeline, BWA aligns the sequencing reads against the reference genome and variant callers detect the SNVs and Indels (7). In this study, SAMtools and GATK were used as variant callers with SSCS and MSCS. Both unfiltered and filtered variants indicated that SAMtools called more SNVs than GATK did owing to its lower internal filtering criteria. Our results demonstrated that filtered and called Indel variants with SAMtools decreased compared with GATK. In parallel, other studies confirmed the high potential of GATK in the identification of true Indel variants (21). In addition, the number of called raw SNVs increased in SSCS with SAMtools and in MSCS with GATK, in both families. The results showed that SAMtools and GATK are capable of calling more raw SNVs in SSCS and MSCS, respectively, in both analyzed pedigrees. Therefore, we suggest that GATK and SAMtools are appropriate callers in MSCS and SSCS, respectively. These results, along with other reports that illustrate the role of GATK in MSCS, such as those by Liu et al., demonstrate that many variants have been lost by SAMtools in MSCS (21). Cornish and Guda compared 30 different pipelines and found that Novoalign plus the GATK Unified Genotyper showed the highest sensitivity with a low number of false positives (22).
Variants were annotated by ANNOVAR, and variants outside the coding regions, as well as those of synonymous coding, were filtered out. Subsequently, known variants with a frequency greater than 1% in ExAC, dbSNP, and 1,000 Genomes were excluded, to reduce the number of potential disease-causing variants (23-25). Homozygosity mapping and linkage approaches prioritized private variants for further analysis based on the two family-based filtering strategies. As shown in Tables V and VI, SAMtools called more private variants that were either missed by GATK or flagged as low-quality.
We have shown that WES is a powerful, efficient, and cost-effective technique that significantly reduces the number of candidate genes in a small number of infertility cases in families with multiple affected individuals. Our research demonstrates that sequencing a small number of samples, while using appropriate filters against public SNVs and in-house databases, is a sufficient approach to detecting private variants. The genomic study of familial male infertility is limited in several respects: the first challenge is the identification of families that contain more than two infertile men who are willing to participate in the study. The selection of appropriate samples for sequencing is also very important. There is a level of uncertainty when detecting genuine fertile men as a result of assisted reproductive technologies, age of incidence, and cultural barriers. Therefore, great care should be taken in the selection of fertile samples.

 4. Conclusion
This study demonstrates two strategies for the WES analysis of familial male infertility cases to suggest a convenient approach to identify potentially functional variants. This, in turn, may further our understanding of the underlying mechanisms behind male infertility.
 
Acknowledgments
The present study was supported in part by a grant from the University of Sistan and Bluchestan (USB), Zahedan, Iran, a grant from the Royan Research Institute, Tehran, Iran, and a grant from the Iran National Science Foundation (INSF).
 
Conflict of Interest
The authors declare that they have no conflicts of interest.