DASAF: An R Package for Deep Sequencing-Based Detection of Fetal Autosomal Abnormalities from Maternal Cell-Free DNA

Background. With the development of massively parallel sequencing (MPS), noninvasive prenatal diagnosis using maternal cell-free DNA is fast becoming the preferred method of fetal chromosomal abnormality detection, due to its inherent high accuracy and low risk. Typically, MPS data is parsed to calculate a risk score, which is used to predict whether a fetal chromosome is normal or not. Although there are several highly sensitive and specific MPS data-parsing algorithms, there are currently no tools that implement these methods. Results. We developed an R package, detection of autosomal abnormalities for fetus (DASAF), that implements the three most popular trisomy detection methods—the standard Z-score (STDZ) method, the GC correction Z-score (GCCZ) method, and the internal reference Z-score (IRZ) method—together with one subchromosome abnormality identification method (SCAZ). Conclusions. With the cost of DNA sequencing declining and with advances in personalized medicine, the demand for noninvasive prenatal testing will undoubtedly increase, which will in turn trigger an increase in the tools available for subsequent analysis. DASAF is a user-friendly tool, implemented in R, that supports identification of whole-chromosome as well as subchromosome abnormalities, based on maternal cell-free DNA sequencing data after genome mapping.


Introduction
Fetal autosomal aneuploidies are one type of abnormalities for chromosome number with a death rate of 6%-11% in newborns. And the most common autosomal aneuploidies are Down's syndrome (trisomy 21) with the incidence of 1 in every 160 newborns causing mental retardation and hypoplasia [1]. Besides whole-chromosome aneuploidies, a considerable number of fetuses are at high risk for subchromosomal abnormalities [2][3][4] that also result in mental illnesses and other abnormalities [5].
The traditional screening for chromosomal abnormalities combined the maternal age, ultrasonographic examination of the fetus, and levels of various proteins or hormones in the maternal blood which refers to traditional noninvasive detection [6]. However, the traditional noninvasive methods are lacking accuracy because they are indirect measures of the underlying chromosomal defect [7,8]. So pregnant women have to choose the invasive methods including chorionic villus sampling and cordocentesis, coupled with fetal cell karyotyping which yield definitive answers. But there is a 0.5% risk of miscarriage which adds additional concern to the pregnant women and their families [9,10].
The discovery of cell-free fetal DNA in maternal serum [11] and recent advances in massive parallel sequencing (MPS) technologies [12][13][14] now enable noninvasive prenatal testing (NIPT) of fetal chromosomal aneuploidies [15][16][17][18], with very high specificity and sensitivity [19][20][21]. In addition to being noninvasive, NIPT requires only 5 mL of maternal peripheral blood for sequencing. Sequences are analyzed using bioinformatics methods to calculate a hazard score, which is then used to determine whether fetal chromosomes are normal or not. Although the standard -score (STDZ) method was originally used, it was later discovered that the accuracy of this method varied depending on the GC content of the chromosomes in question [15]. More specifically, the coefficients of variance (of measuring the percentage of representation of each chromosome) were much larger for chromosomes 18 and 13 than for chromosome 21 [15,16]. This variation in accuracy is linked to the difference in sequencing efficacy as a function of chromosome size and GC content. In recent years, many methods have emerged to solve the aforementioned problem, including a GC correctionscore (GCCZ) method [21], internal reference -score (IRZ) method [20], and the noninvasive fetal trisomy (NIFTY) test [19], as well as the method of Srinivasan et al., henceforth referred to as the subchromosome abnormality -score (SCAZ) method [22]. The first three methods are similar to the standard method ( -score) for identifying abnormalities in whole chromosomes, while the last method is used to identify subchromosomal (i.e., chromosomal regions) losses and gains. Lau et al. indicated that the standard -score (STDZ) method accurately detects trisomy 21 early in pregnancy of 11 weeks with low accuracy for other aneuploidies, being 0% for trisomy 13 and 40% for trisomy 18, while the GC correction with LOESS regression method (GCCZ) is more accurate than STDZ but still with low detection rate for trisomy 18. And the adjusted method using -scores with an internal reference (IRZ), which corrects for GC bias and sequencing efficiency, substantially improved the performance of the test [20]. On the other hand, Verweij et al. investigated the attitudes among pregnant women regarding NIPT for the detection of trisomy 21 (T21): they had a positive attitude regarding NIPT for detection of T21, and more than 50% of them who rejected the traditional screening would accept NIPT if available [23]. However, although NIPT has become increasingly popular and acceptable and subsequent data analysis algorithms have emerged, there are no tools currently available to implement these data analysis methods. In the present study, we developed an R package, DASAF, that implements the three most popular trisomy detection methods (STDZ, GCCZ, and IRZ) and one subchromosome abnormality identification method (SCAZ). We have also included a fetal gender prediction module in the DASAF package. With the cost of DNA sequencing declining and with advances in personalized medicine, we believe that the demand for NIPT will increase, which will undoubtedly trigger an increase in the tools available for subsequent analysis.

Materials and Methods
This study was approved by the Independent Ethics Committee of Shanghai Clinical Research Center. The reference data used here consists of DNA sequencing data from one hundred and twenty pregnant women from Huzhou Maternity & Child Care Hospital located in Huzhou, Zhejiang, China. All data were produced by Illumina HiSeq2000 for 100 bp pair-end with 7 × 10 6 to 17 × 10 6 sequence read pairs per sample.
The sequencing reads were aligned to the human genome assembly hg19 with Bowtie short read aligner (version 1.1.2), allowing for two base mismatches at most when aligning [24]. Only uniquely mapped reads were kept.
Before using DASAF, sequencing data should be aligned using the above method, which is independent of the DASAF software and needs to be completed by the users themselves. The results file from Bowtie is used as input for DASAF. All the reference datasets are described in Table 1. A typical DASAF workflow involves two procedures: mapping read statistics and autosomal aneuploidy prediction ( Figure 1).

Read Mapping Statistics.
Read mapping statistics produce two files: one contains the unique mapping read counts for every chromosome and the other contains the mapping location for every unique read. The normalized chromosome ratio (NCR) is generated according to the following equation for every chromosome in each sample: NCR is the ratio of number of reads uniquely mapped to the specific chromosome divided by the total number of reads uniquely mapped to all autosomal chromosomes [15,25]. If the GCCZ method is used, the GC content for every chromosome is calculated from the mapping results.  The results from (a) are used to calculate the risk score using any of the four methods implemented in (b). STDZ: standard -score; GCCZ: GC correction -score; IRZ: internal reference -score; and SCAZ: subchromosome abnormalities -score.

Autosomal Aneuploidy Prediction
In the standard -score (STDZ) theory method, a hazard ratio of the -score is calculated to determine whether the fetal chromosome is normal or not: where NCR is the ratio of the sequence counts uniquely mapped to the specific chromosome and the total number of the sequences uniquely mapped to all of the autosomal chromosomes, NCR is the average NCR of chromosome in the reference samples, SD is the standard deviation for NCRs of chromosome in the reference samples, and is the specific chromosome number, that is, 13, 18, and 21 [15]. For the average value and standard deviation values for the NCRs, one can use the reference files (NCR ref.txt) contained in the DASAF package or calculate them based on one's own samples. The -score is a number indicating how far an observation deviates from the average in a population [26]. Usually, a -score of 3 is selected as threshold to determine whether the fetus is normal or not [22].

GC Correction
where NCR GC is the NCR value after GC correction, NCR is the original value, and GC average ref and Slope ref are the mean values of references' chromosomal GC content and the slope of linear regression from the reference samples [21]. Then, the mean and SD of the GC-corrected NCR were calculated for the reference dataset and the -score was calculated for the chromosome of the sample tested using (1) with a -score cutoff of 3.

Internal Reference -Score Method.
To minimize the sequencing bias (stemming from differences in GC content), Lau et al. presented a -score method that relies on an internal reference chromosome [20]. They showed that using chromosomes 4, 8, and 14 as internal reference chromosomes provided the most accurate results for the detection of trisomy 13, trisomy 18, and trisomy 21, respectively. The method is as follows: the comparative NCR is calculated using the value from the internal reference as NCR /NCR IR , where IR is the internal reference chromosome for chromosome . The -score is also calculated by (1) that the IR adjusted NCR value for the test sample subtracts the averaged IR adjusted NCR values from the reference samples and the difference is then divided by the standard deviation from the IR adjusted NCR values for the reference samples. A -score of 3 was selected as threshold for the diagnosis of trisomy in chromosome of the testing sample [20].

Subchromosome Abnormality -Score Method.
In addition to whole-chromosome abnormalities, subchromosome losses and gains are also important components of chromosomal diseases [4]. The subchromosome abnormalityscore (SCAZ) is a method used to identify abnormalities for chromosomal regions with lengths between 100 kb and 1 Mb [22]. In the first step, positions uniquely mapped to the genome are retrieved and counted as tags. And the whole genome was divided into continuous bins with length of 1 Mb and 100 kb and tags were assigned to individual bins for the following analysis. Then GC content percentage of each bin was calculated to rank the bins across the entire genome. And then every bin was normalized using the ratio of tags within the bin to the sum of the tag counts in bins with the nearest GC content percentages. Bins with nearest GC content percentages include 10 bins of 1 Mb length and 40 bins of 100 kb length. The equation is as follows: where BRV is the ratio for the th bin for chromosome and Tags is the count of tags in the th bin for chromosome . represents the bins with length of 100 kb and 1 Mb. Further, every BRV was examined for deviation from the median values collected across all the reference samples which is similar to the standard -score method, while the median absolute deviations (MAD) were adjusted to MAD (i.e., MAD was multiplied by 1.4826); here is 1.4826. Consider The absolute values of -score larger than 3 indicate that there were CNVs in fetal chromosome for the specific genomic regions [22].

Results and Discussions
We built the DASAF R package, which supports three existing methods for identifying whole chromosome abnormalities and one for identifying subchromosome abnormalities from MPS data. We then compared the running time and identification accuracy of the four methods.

Comparison of Chromosomal Abnormality Detection
Methods. All the detection methods used here were derived from existing algorithms and their accuracy has been tested previously [19,20]. Here, we therefore only list the previously reported results for these algorithms. Lau et al. provided detection rates and false-positive rates for the three wholechromosome trisomy detection methods. Their research revealed that the false-positive rates were 0 for all the three methods and the method of IRZ was the most sensitive, with a 100% detection rate for all trisomies examined (13, 18, and 21). For the method of STDZ, the detection rate was 100% for detecting trisomy 21 but only 40% for trisomy 18 and almost 0% for trisomy 13, while the GCCZ method with a detection rate of 100% for trisomy 21, 90% for trisomy 18, and 100% for trisomy 13 was better than the standard method but worse than the IRZ method [20]. Jiang et al. also evaluated the performance of these three methods for 903 cases and found that the Coefficient of Variation (CV) for the STDZ method was larger than that for the other two approaches among clinically relevant chromosomes (13, 18, and 21). Thus, the STDZ method has poor sensitivity for the detection of trisomies 13 and 18. However, the performance of the GCCZ approach demonstrated over 99% sensitivity and specificity for the detection of trisomies 13, 18, and 21, while the IRZ approach displayed CV larger than GCCZ but smaller than STDZ for chromosomal trisomies 13, 18, and 21 [19]. In summary, the adjusted methods (GCCZ and IRZ) more accurately identify trisomies than the STDZ method. It was also reported that the SCAZ method, which identifies chromosome CNVs, can accurately detect losses and gains for chromosomal regions [22].

Evaluation of Diagnostic Accuracy as a Function of Sequencing Depth.
In order to evaluate the effect of sequencing depth on diagnostic accuracy, we randomly subsampled the 100 bp pair-end (PE) sequencing data at read counts of at least 3 M, 5 M, 7 M, and 12 M. Using the STDZ method, cases were diagnosed as T21-positive or T21-negative. Importantly, we found that, even at a read count of 3 M, T21 was accurately   diagnosed, which suggests that the cost of sequencing can be considerably reduced by decreasing the sequence coverage. The results shown in Figure 2 demonstrate that the -scores for all the positive samples are larger that 3 (above the horizontal line of = 3). STDZ and IRZ ran faster than the other methods if the NCR values for references were prepared beforehand. The GCCZ method requires the user to calculate the GC content for every chromosome, which consumes a considerable amount of time. The SCAZ method had the longest runtime because the BRV needs to be calculated for every bin by counting the tags. While all running times were acceptable, these times can be dramatically reduced by decreasing the sequence read counts to 3-5 M ( Table 2).

Conclusions
We developed an R package that supports chromosomal abnormality detection. For chromosomal abnormality detection, users can select one of four supported methods or, for whole chromosomal abnormality detection, summarize the results of the three available methods (i.e., average the three -scores) for detection of trisomies 13, 18, and 21. We chose a -score threshold of 3 to predict fetal chromosome abnormalities. The reference datasets under the directory of data in the package can be updated or replaced by users as the samples increase, which can promote the accuracy of these methods. A detailed vignette is included with the DASAF package to assist nonexperts in the field (http://lifecenter.sgst.cn/dasaf/).
The cost of high-throughput sequencing has decreased dramatically over the past few years, thus increasing its utility in clinical practice [27,28]. Noninvasive prenatal diagnosis is the most widely used method for detecting trisomic abnormalities or the loss or gain of chromosomal regions, and an increasing number of pregnant women are benefitting from this technology. In August 2014, noninvasive prenatal DNA diagnosis finally obtained legal status in China following the approval of the registration of second-generation genesequencing diagnostic products. This represents a major advance in the field of prenatal screening that will undoubtedly benefit numerous pregnant women and their families.