WGDdetector: a pipeline for detecting whole genome duplication events using the genome or transcriptome annotations

With the availability of well-assembled genomes of a growing number of organisms, identifying the bioinformatic basis of whole genome duplication (WGD) is a growing field of genomics. The most extant software for detecting footprints of WGDs has been restricted to a well-assembled genome. However, the massive poor quality genomes and the more accessible transcriptomes have been largely ignored, and in theoretically they are also likely to contribute to detect WGD using dS based method. Here, to resolve these problems, we have designed a universal and simple technical tool WGDdetector for detecting WGDs using either genome or transcriptome annotations in different organisms based on the widely used dS based method. We have constructed WGDdetector pipeline that integrates all analyses including gene family constructing, dS estimating and phasing, and outputting the dS values of each paralogs pairs processed with only one command. We further chose four species (Arabidopsis thaliana, Juglans regia, Populus trichocarpa and Xenopus laevis) representing herb, wood and animal, to test its practicability. Our final results showed a high degree of accuracy with the previous studies using both genome and transcriptome data. WGDdetector is not only reliable and stable for genome data, but also a new way to using the transcriptome data to obtain the correct dS distribution for detecting WGD. The source code is freely available, and is implemented in Windows and Linux operation system.


Background
Polyploidy or whole genome duplication (WGD) is just like what it sounds: an event of nondisjunction during meiosis which drives species diversification and evolutionary novelties with additional copies of the entire genome [1][2][3]. As a common phenomenon in plants, all extant seed plants have experienced at least one ancient WGD, and many flowering plants have undergone multiple rounds of paleopolyploidy [4,5]. WGD has long been considered as the major force for rapid genome evolution [6][7][8], which could increase organism complexity, enhance adaptation through dosage effect and induce the speciation and biodiversifcation by imme diately producing reproductive isolation with other relatives [9][10][11]. Moreover, WGD also plays the key role in the domestication of many crops, such as maize, wheat and cotton [12]. For these reasons, there is an increasing interest in detecting the bioinformatic basis of whole genome duplication events.
There are three main methods to search for evidence of WGD [13]. The most straightforward evidence for WGDs is the presence of large syntenic regions within a genome, while these methods need a well-assembled genome, the nearer to chromosome level the more accurate of the results [14,15]. With a growing number of published draft genomes, two other methods based on phylogenetics [4,16] and distribution of pairwise paralogs synonymous substitutions per synonymous site (dS) are more suitable [17,18]. For the former, the WGDs are estimated through the gene count data where the number of gene copies in various gene families across a group of taxa along the phylogeny is counted with the gene birth and death rates in consideration [16]. And the dS based method assumes that each gene family has the constant rates of birth or loss death [19], while WGDs violate this assumption and produce peaks in cumulative distributions of pairwise dS between paralogs within a genome [18]. Recently, the dS based method has become the most common and widely used approach to inferring WGD. Theoretically, peaks in cumulative distributions of pairwise dS between paralogs within the same species should be universal in both genome and transcriptome annotations. Here, we just focused on the dS based method to develop a technical tool for detecting WGDs and trying to break its limitation for utilization of only genome annotations.
Within dS based method, the core and initial step is to identify the pairwise paralogs among the genome, and then, to estimate the distribution of fourfold synonymous third-codon transversion rate (4DTV) or dS between paralogs pairs to determine the WGDs. There are two main approaches to identifying the pairwise paralogs. One is to use the combined gene similarity and gene order information to identify syntenic pairwise paralogs, through many software including ADHoRe [20], DAGchainer [21], ColinearScan [22], MCscan [23], MCScanX [24], SyMAP [25], and so on. The gene order information is unavailable in the poor quality genome assembly or the transcrip tomes, which will limit the usage of those software. The other is to use a gene family based approach to identifying pairwise paralogs, which does not need the gene location information and can be suitable for most genomes. However, how to convert the gene family dS to represent the pairwise paralogs dS is complex, since the large gene family need to correct the redundant dS values. Those analyses or approaches are mainly achieved by in-house scripts [18,26], which are difficult to transfer the same analyses for other species or web servers [27]. Therefore, the repeated attempts are needed, which would cause the wastes of time and resources. To phase those problems, we construct this WGDdetector pipeline for WGDs detecting that integrates all analyses processed with only one command, which includes gene family constructing and dS estimating and phasing, and outputting the dS values of each paralogs pairs.

Implementation
WGDdetector is written in Perl. BioPerl must be installed and other seven easily installed software (BLAST [28], MMseqs2 [29], BlastGraphMetrics [30], MCL [31], MAFFT [32], PAL2NAL [33] and R [34]) are also needed. Their function was used in our pipeline and the major steps were exhibited in Fig. 1, and the detailed process was described as follow:

Gene family constructing
In this step, WGDdetector supplied two methods to detect the gene similarity: BLAST [28] or MMseqs2 [29] with an e-value cut-off of 1e-10. Here, we recommend selecting MMseqs2, as it can run 10,000 times faster than BLAST, and the results were similar. Then the Blas-tGraphMetrics [30] was used to phase the similarity file, and the followed MCL [31] was selected to construct the gene families based on the Markov Cluster algorithm.

dS value estimating
WGDdetector automatically aligned the protein and CDS sequence within each gene family using MAFFT [32] and PAL2NAL [33], and assigned the corresponding dS values for each pair paralogs (gap-stripped alignment length > 90 bp) within each gene family via the 'Bio::Align::DNAStatistics' Perl module based on the Nei-Gojo bori algorithm.

dS correction for redundant
As the above estimates, a gene family of n members originated from n-1 retained single gene duplications and generated the number of possible pairwise comparisons is n(n-1)/2. To correct the redundancy of dS values, we used a slightly modified strategy as described in Arabidopsis [18] and Norway spruce analysis [35]. We used the dS as a distance measure, and constructed a tentative phylogenetic tree with an average linkage clustering algorithm using the 'hclust' R module. A series of clusters (from 1 to n, n is the gene numbers within one family) were generated by the 'cutree' function for each gene family. Subsequently, they were divided into the subfamilies with the dS values less than 5 and each subfamily contained as many genes as possible. Then, a tentative phylogenetic tree was constructed again for each subfamily, and 'cutree' was used to intercept only two child clades. We summed the dS values for all combinations between the two child clades, and weighed the number of combinations to represent this subfamily, which corresponded to a duplication event. Finally, we collected all the dS values of each subfamily and supply the R script to plot the distribution.

Results
Four organisms' genome or/and transcriptome datasets were selected to evaluate the performance of WGDdetector, including three plants (Arabidopsis thaliana, Juglans regia and Populus trichocarpa) and one frog (Xenopus laevis) ( Table 1 and Additional file 1: Table S1). For the genome datasets, a total of 27,301, 32,436, 39,410 and 41,073 genes satisfied our criteria in A. thaliana, J. regia, P. trichocarpa and X. laevis, respectively: retaining the longest coding sequence (CDS) for each gene, removing CDS with premature stop codons and those protein sequences < 50 amino acids (AA). For the transcriptome datasets, the raw reads were download from NCBI SRA and assembled by The format "S (X + Y)" represent the total elapsed time (S), the protein clustering elapsed time (X) and the dS calculating elapsed time (Y) c In the software ADHoRe and MCScanX, all the homologous genes were recorded. In the WGDdetector pipeline, only the sub-gene families were recoded Trinity v2.5.1 [36] with the default parameters except "--trimmomatic" and "--normalize_reads". The constructed transcripts were filtered by the SeqClean [37] to remove contamination, and then the TransDecoder v5.3.0 [38] was used to identify candidate coding regions. The candidate alternative splicing were filtered by CD-HIT-EST with the parameter '-c 0.9′ [39], and the proteins with length less than 50 AA were further removed. A total of 23,495 and 20,354 transcripts were obtained for the following analysis in A. thaliana and P. trichocarpa, respectively.
All datasets with a gene number ranging from 20,354 to 41,073, showed a memory usage approximately 6~35G and the elapsed time around 5-19 h ( Table 1). As our pipeline could use multiple CPUs, this elapsed time would be shorter if more CPUs used. To evaluate the performance of WGDdetector, ADHoRe and MCScanX were selected for comparisons. The general similar trajectories of the density or histogram were observed in all the datasets implemented by WGDdetector, ADHoRe and MCScanX, and different software have different superiority at the recent or ancient WGD events (Fig. 2). All the first peaks were coincidence by different approaches within each species, which indicated high sensitivity and accuracy in the detection of recent WGD event using WGDdetector, based on both genome and transcriptome datasets. For A. thaliana, a major peak (the second) with a long range (0.7~2) was detected using both ADHoRe and MCScanX, which was hard to discern the ancient WGD event. While, the result of WGDdetector showed two peaks (~0.6 and~1.9), representing the 1R and 2R WGD events within A. thaliana and coincident with the previous studies [18,40]. In the other three species, WGDdetector also showed a high sensitive on the detection of ancient WGD event, as a more obvious second peak detected than the other two software. But we also found slightly larger dS values in the second peak in WGDdetector than the other software, as detected in P. trichocarpa

Discussion
As the methodological distinctness at the dS distribution obtaining, WGDdetector elapsed more time and memory than ADHoRe and MCScanX (Table 1). This was mainly caused by the most time consuming step that WGDdetector calculated the dS values between all the possible homologous gene pairs within each gene family. While ADHoRe and MCScanX needed the gene order information to identify the synteny gene pairs and thereby a small number of dS values were calculated [24]. In the results of accuracy evaluation, WGDdetector showed a high accu racy and more sensitive of detecting the recent WGD events. In the genome data of J. regia or the transcriptome data of A. thaliana and P. trichocarpa, WGDdetector was also detected noise signal peaks (near the origin), which might reflect the unmerged allelic haplotypes in the genome data [41] or the alternative splice transcripts within the transcriptome data that was still retained. Our results of the genome data of A. thaliana also proved the distinct first and second peaks, rather than a long range peak detected in MCScanX and ADHoRe, which reflecting the high performance of detecting the ancient WGD events in WGDdetector. The second peaks in each dataset showed a little difference in different software. We speculated that this difference might be caused by dS saturation when the dS value > 1 [42], and the higher sensitive performance in the detecting ancient WGD in WGDdetector than that in ADHoRe and MCScanX.

Conclusions
The WGDdetector was designed as a user-friendly pipeline with a very simple command which only needed the CDS and protein files. This pipeline integrated the gene family constructing, dS estimating and hierarchical clustering, dS correcting and distributing plotting. This methodology eliminates the limitation of gene order information and is more suitable for the well/poor quality genomes and transcriptomes. In our practice based on the genome and transcriptome datasets, WGDdetector showed a high perfor mance in the detection of recent and ancient WGD events and matched well with the previous studies and/or the software of ADHoRe and MCScanX. With the development and rapidly declining cost of next-generation sequencing (NGS) technologies and third-generation long-range DNA sequencing, more and more species would be resolved by sequencing their genomes and/or transcriptomes [43,44]. Totally, WGDdetector gives a reliable and acceptable way to infer WGD event using either genome or transcriptome data by the dS-based method, and will help to accelerate our understanding of the evolutionary history of WGDs within all organisms.

Availability and requirements
Project name: WGDdetector.
License: GNU GPL v3. Any restrictions to use by non-academics: none.