text:
Identification of an early diagnostic biomarker of lung adenocarcinoma based on co-expression similarity and construction of a diagnostic model

Background The purpose of this study was to achieve early and accurate diagnosis of lung cancer and long-term monitoring of the therapeutic response. Methods We downloaded GSE20189 from GEO database as analysis data. We also downloaded human lung adenocarcinoma RNA-seq transcriptome expression data from the TCGA database as validation data. Finally, the expression of all of the genes underwent z test normalization. We used ANOVA to identify differentially expressed genes specific to each stage, as well as the intersection between them. Two methods, correlation analysis and co-expression network analysis, were used to compare the expression patterns and topological properties of each stage. Using the functional quantification algorithm, we evaluated the functional level of each significantly enriched biological function under different stages. A machine-learning algorithm was used to screen out significant functions as features and to establish an early diagnosis model. Finally, survival analysis was used to verify the correlation between the outcome and the biomarkers that we found. Results We screened 12 significant biomarkers that could distinguish lung cancer patients with diverse risks. Patients carrying variations in these 12 genes also presented a poor outcome in terms of survival status compared with patients without variations. Conclusions We propose a new molecular-based noninvasive detection method. According to the expression of the stage-specific gene set in the peripheral blood of patients with lung cancer, the difference in the functional level is quantified to realize the early diagnosis and prediction of lung cancer.

Background Among patients with non-small cell lung cancer(NSCLC), which accounts for approximately 85% of all lung cancer cases, nearly 50% have lung adenocarcinoma, which is the most common lung cancer [1]. Lung adenocarcinoma usually begins at the outer part of the lung tissue. Early lung cancer does not yield significant clinical manifestations; only in the late stage will there be a gradual emergence of chronic coughs, bloody sputum, and other symptoms. In that case, some early symptoms (such as fatigue, shortness of breath, or upper back and chest pain) are likely to be overlooked [2].
The clinical diagnosis of lung cancer is mainly dependent on an X-ray; lung cancer lesions often appear as abnormal shadows on X-rays. However, in up to 25% of lung cancer cases, the chest X-ray does not reveal any abnormal lesions and returns a perfect "normal" diagnosis. If you still suspect cancer, you can also use other more sensitive diagnostic methods, including computed tomography (CT) or MRI scans [3,4]. According to the results, the doctor may wish to obtain a lung tissue sample by using a puncture to confirm. However, there are many controversies associated with this invasive detection method. First, nobody can guarantee that puncture sampling will always obtain tumour cells, and the invasion of the sampling process has a potential risk of cancer metastasis. Therefore, a new blood detection method called liquid biopsy has emerged, and it can accurately detect the expression of specific genes in lung adenocarcinoma under noninvasive conditions. According to the appearance of particular genes in the peripheral blood of patients, this new detection method will achieve the early diagnosis of lung cancer with 91% accuracy and longterm monitoring of therapeutic response [5]. Beyond that, this gene detection method has a high sensitivity and specificity compared with traditional detection methods while avoiding the risk of cancer cell proliferation due to invasive exposure.
Once the doctor clinically confirms that a patient has lung cancer, the doctor will classify the result into different stages based on the progress of lung cancer, such as whether it has spread, and if any other articles may be involved [6]. Multi-staging helps the direct treatment in a more appropriate manner, neither under treating a malignancy or over treating it and causing more harm than good. According to the malignant degree of carcinoma, there will be four stages: stage I, stage II, stage III and stage IV [7]. Stage I: the cancer is localized and has not spread to any lymph nodes. Stage 2: the cancer has spread to the lymph nodes, the lining of the lungs, or the significant passageways of the lungs. Stage 3: the cancer has spread to nearby tissue. Stage 4: the cancer has spread (metastasized) to farther reaches of the body. Lung adenocarcinoma is classified by the qualitative analysis of cancer tissue after diagnosis, which means that in the diagnosis stage, there is no distinction between the stages of lung adenocarcinoma [8].
This process also leads to many lung cancer patients not knowing whether they are probably in an advanced stage. By combining the genetic analyses at the cellular and functional level, experimenters identify specific genes and mutations in the four stages. These genes may be the markers of different stages of lung cancer or the target for personalized treatment. Additionally, the variability function explains the biological difference between the various stages and the malignant mechanism of the fatal development of cancer to a certain extent.
In conclusion, the traditional method of finding disease markers is mostly dependent on the expression level of single genes. The assumption is that each gene is relatively independent, so the doctor will use the selected gene's characteristics to design diagnostic kits, make the combination of chips or sequencing, and predict the risk of illness based on gene expression [9]. However, in organisms, the gene is not relatively independent, and there is a functional interaction. By considering the patient's molecular-level and functional-level changes from a higher perspective, this study converts gene expression information into a functional imbalance. On the one hand, this method overcomes the instability of cross-sample single gene markers. On the other hand, it prompts the underlying pathogenesis from the functional level, while the gene with the specific function we collected is likely to be a clinically significant therapeutic target or diagnostic marker. More importantly, our analysis found 12 genes that not only can predict lung cancer patients in the early phase but can also present dynamic changes during cancer development. Moreover, patients carrying variations in these genes tend to have a poor outcome in terms of survival status compared with patients without variations.

Data remodelling and grouping
The RNA-seq data GSE20189 were downloaded from the GEO database [10]. We first split the lung cancer data into groups based on the clinical data. The original data consistes of 22,277 genes and 162 samples. There were 81 controls, 28 patients in the early phase (stage I) and 53 patients in the late phase (stage II-stage IV).

Data normalization
Initially, we accomplish the data pre-processing and standardization, and we removed genes and samples with the ratio of missing values higher than 10%. The missing values within remaining samples are replaced with the mean values of the corresponding genes in the other samples. We calculated the mean value and standard deviation of each gene in the control group. Then, we applied Z-score normalization [11] to all samples after which the expression of the gene in the control group obeys the standard normal distribution, with a mean value of 0 and a variance of 1. Therefore, for gene I, if there is no difference between different stages of the lung adenocarcinoma samples, the case groups should be subject to a normal distribution. Otherwise, we believe the gene i in a specific stage is obviously different compared to the control group; this differentially expressed gene is likely to act as a crucial biomarker in early diagnosis.

Stage grading-specific gene extraction
In the evolutionary process, most of the genes are conservative and commonly do not have an apparent difference in expression under the disease signal stimulation. If the gene is significantly associated with lung adenocarcinoma, then different expression should be observed in at least one stage group. Therefore, we used the coefficient of variation to assess the fluctuation of genes across the lung adenocarcinoma samples. The coefficient of variation can be calculated by Eq. 1.
Mean is the mean expression of the gene in all lung adenocarcinoma samples, and sd is the corresponding standard deviation. According to the distribution of all the coefficients of variation in the genes, we only screened genes whose absolute value of CV rank in the top 50% as genes that may be associated with lung adenocarcinoma. The remaining 50% of the genes are due to a small fluctuation in the vicinity of 0; therefore, these genes are supposed to have little impact on lung adenocarcinoma. To identify the specific genes for early or late phase, we used the limma method [12] to evaluate the differences compared with control group. The significance threshold of the P value was set to 0.05, and the cut off of the absolute value of the fold change was set to 1. The significantly differentially expressed genes in the early phase were marked as Δ0. The differentially expressed genes related to the late phase were marked as Δ1. The intersection between early and late phase specific genes was marked as Δ2. These genes can distinguish normal and cancer patients in the early phase and present dynamic changes with cancer development.

Co-expression correlation analysis
Since the interaction between genes changes from one stage to the next stage of lung adenocarcinoma, the two (1) CV = sd mean correlated genes tend to have a common biological effect in the non-disease state from the biological point of view. It appears as a mutual synergy or interaction. However, in the disease state, the function of the genes is abnormal, which may results in the change of correlation. Therefore, we examined the expression correlation in Δ0, Δ1 and Δ2 using the Pearson correlation coefficient [13] (greater than 0.5 is positively correlated, and below − 0.5 is negatively correlated).

Unsupervised clustering analysis
We utilize correlation analysis to construct the correlation coefficient matrix between genes and take advantage of hierarchical clustering [14] to achieve unsupervised clustering for samples and genes. From the results of unsupervised clustering, we observe and analyse the effect of distinguishing cancer patients from controls at the gene level. On the one hand, we can verify the Δ2 gene's ability to identify the specific stage of the lung adenocarcinoma sample; on the other hand, we can observe the gene expression patterns at different stages, such as genes in stage i that have a high expression in a group while having a low expression in stage j. The transformation of this expression pattern also indicates that the progression of the disease affects the function of gene regulation, which leads to some function level changes. The clustering results are visualized using heatmap thermal maps [15].

Specific and non-specific co-expression network analysis
Because the correlation between gene expression is different through disease status, the status specific networks should also reflect a significant difference in network characteristics from the perspective of the system network. We constructed specific networks of control, early and late phase based on the co-expression relationship between genes from Δ0 and Δ1. The non-specific network is constructed using genes from Δ2. If there are co-expression relationships between the two genes, there is an edge between genes. Because of the co-expression network constructed by gene co-expression in different disease states, it shows a significantly different topological nature, suggesting that the signal transmission efficiency of the system network is remarkably unlike in different vicious grades. Hence, we start analysing from six topological properties of network connectivity [16], namely, ASLP (average shortest path), Closeness Centrality, Cluster Coefficient, and Degree Distribution. If the edge of the network is missing, that is, the co-expression relationship between the genes disappears, the average shortest path of the network increases, and along with a decrease in the Degree Distribution, Clustering Coefficient, and Closeness Centrality, the network signal transmission performance is declining as well. The specific network of the control status reflects intrinsic relationships among genes. The specific network of the early phase reflects the changes in the initial stage of cancer. The specific network of the late phase reflects more variations occurring among genes. On the other hand, as the overlapped genes (Δ2) present deviations in both early and late phase, these genes can be regard as biomarkers associated with cancer progression. Therefore, we also construct a non-specific network using Δ2. Eventually, we use the degree of distribution of gene nodes in the network to evaluate the importance of genes; the higher the degree of the gene, the more adjacent genes are becoming affected by the abnormal gene. We convert the degree of all genes through the sigmoid function [17] to the weight of 0-1, and the weight of the gene not in the network defaults to the minimum.

Functional pathway enrichment
To further analyse the biological functions involved in the specific genes in the different disease states from the functional level, we apply the functional enrichment analysis on the Δ2 genes. By using Fisher's exact test [18], the significant pathways selected are the functions for regulation of these specific genes. Since these genes have dissimilar expression patterns in the different stages, at the same time, the genes that have a co-expression correlation tend to participate in the same biological ability; we deduce that these pathways can better indicate the efficacy of the genes in distinguishing stages. Beyond that, these pathways with abnormal functional levels in different stages can be used to explain the mechanisms of disease progression, and these pathways may contain potential drug targets or diagnostic markers.

Identification of significantly deviated pathways
Suppose we used the intersection gene to obtain N pathways through function enrichment. We first identify the differentially expressed genes in each gene pool and then use the inverse cumulative distribution function to convert the P value of ANOVA [19] to the Z value and multiply the weight of the gene. The differentially expressed genes in pathway P are calculated from the formula to calculate the deviation score A (P) of the pathway [20], as follows: In the calculation process, we first sort the differentially expressed genes from large to small; therefore, the larger the Z value is, the more significant is the gene expressed. Suppose that there are k differentially expressed genes in pathway P, then 2, 3……k genes are used to calculate the Z score mean value. When the Z score of the t genes is the maximum, the corresponding t genes have the greatest contribution to pathway P, and we calculate the deviation score A(P) of pathway P in the disease state.
To eliminate the influence of the size of the pathway itself, we calibrate the A(P) score with the correction method of the random perturbation principle. For the pathway P, the deviation score is A(P), and we recalculate a new A(p)' by k sets of genes in the random selected pathways. After 10,000 random cycles, we calculate the mean μ and standard deviation of the A(p)' background distribution. A calibrated A corrected is obtained using Eq. 3.
There are some problems to be considered when calculating the deviation score of pathway P. First, the number of genes differentially expressed in the pathway is large, but not all of them have a significant impact on the pathway, such as some genes belonging downstream of this pathway. The difference in expression of these genes is likely due to changes in upstream abnormal signals. In contrast, certain upstream genes, important enzymes, transporters, and other genes that have important regulatory effects may have a greater impact on the pathway. Second, the number of the genes in each pathway is significantly different; in order to eliminate the impact of the size of a certain pathway itself, we utilize the random punctuation process.

Early screening gene filtering using RFE
Significant pathways identified from random permutation present functional differences in the early and late phases of lung cancer development. This further suggests that the genes regulating these pathways have important roles in lung cancer. On the one hand, lung adenocarcinoma-related genes exhibit differences in the level of expression under the disease signal stimulation. On the other hand, the co-expression relation changes between the genes and thus affects the levels of downstream functional pathways. These genes involved in the regulation may be potential targets for lung cancer treatment or new clinical monitoring and diagnostic indicators. To accurately identify the optimal combination of gene features, we used the recursive feature elimination (RFE) algorithm [21] to filter genes. Finally, we screened out the risk-related genes of lung cancer to train the diagnostic model.

Establishing an early diagnostic model
Actually, as the SVM performs better in the binary classification, we combined stage I and stage II as the early benign group, while stage III and stage IV as the advanced malignant group. To distinguish patients belonging to benign and malignant groups using risk-related genes as features, we utilized the supervised classifier support vector machine (SVM) to train a diagnostic model [22]. Default parameters were used to initialize the model, including the RBF nonlinear kernel function, gamma of 0 and so on. A grid search algorithm was used to optimize all parameters [23]. Cross validation was used to calculate true positive rates and false positive rates, and an ROC curve was drawn to estimate the model's performance.

Survival analysis of independent validation data
To further validate the early risk genes, our screen can not only achieve the early diagnosis of lung cancer but also estimate the prognosis of patients to a certain extent, to provide the basis for the individual treatment strategy. We downloaded the lung adenocarcinoma samples from the TCGA database as independent validation data. We used a Cox regression analysis of the lung adenocarcinoma sample's overall survival data based on the risk genes. The P value was calculated to estimate the correlation between the survival time and risk gene variance.

Phase-specific gene recognition
We used lung adenocarcinoma samples of the early phase and late phase late to compare with healthy control groups, respectively. According to the result of the limma algorithm, 866 early phase-related differentially expressed genes were identified, including 136 up-regulated genes and 730 down-regulated genes; 913 late phase-related differentially expressed genes were identified, of which 419 were down-regulated genes and 494 were up-regulated genes. The distribution of the two groups of differentially expressed genes at the P value and logFC levels is shown in Fig. 1a.
In Fig. 1a, the down-regulated genes are labelled with green, and the up-regulated genes are labelled with red. It can be observed from the figure that the down-regulated genes are dominant in the early stage, whereas the up-regulated genes are dominant in the late phase as the cancer progresses. This suggests that an increasing number of genes are up-regulated during the malignant process of lung adenocarcinoma. The intersection genes of the early phase and late phase are shown in Fig. 1b.
We found 108 shared down-regulated genes in early and late lung cancer-related genes as well as 90 shared up-regulated genes. This indicates that the expression levels of these 198 genes in early and late lung cancer patients are both different. In the process of the malignant progression of lung adenocarcinoma, the 198 genes showed a dynamic change. On the one hand, this dynamic pattern of genes associated with the disease can be used as a clinical indicator to achieve an early diagnosis. On the other hand, the function of these genes is also likely to regulate the progress of cancer.

Co-expression correlation analysis
If the intrinsic correlation of a pair of genes is absent after entering a disease stage, this indicates one of these genes is associated with lung cancer. In contrast, if a pair of non-related genes, appears to be correlated in disease status, then it suggests that this two genes are likely initially divided into two parallel pathways. In the process of lung cancer progression, one of the pathways is dysfunctional; therefore, another pathway is activated as a compensatory pathway, thus reflecting the new gene co-expression relationship [24]. Other than these two groups of co-expression patterns, the remaining groups in any stage in the stable co-expression are those genes that always maintain the relationship along with disease progression. Therefore, it is essential to identify the coexpressed genes of these three different patterns, which is useful for explaining the specificity of the lung adenocarcinoma stage. In cancer studies, co-expression correlations between genes are more important because the correlation between genes varies dynamically with cancer progression. This dynamic change, on the one hand, provides a basis for the pathogenesis of cancer progression and, on the other hand, is also an important feature of dynamic monitoring of patient status. We use the Pearson algorithm to internally compute the correlation coefficient between any two genes for each stage-specific gene set. Table 1 calculates the number of genes in the phasespecific gene sets that have a remarkable correlation. Most of the relevant genes are positively correlated, and a few genes are negatively correlated. At the same time, we found that there are 3015 gene pairs in the three phases' intersections across the normal and cancer status. These gene pairs are stable in the expression of relevance within the three phases; 3015 gene pairs involve 164 genes in total, and therefore, we investigate the 164 genes in the three phases of the co-expression state, as shown in Fig. 1c.
We observe that the 164 genes of the intersection reflect apparent differences in the four stages. Most genes are still positively correlated, but the degree and type of correlation between any two genes are not consistent in the different stages. This result again suggests that the covariant correlation between the two genes changes with the progression of lung adenocarcinoma.

Unsupervised clustering analysis
Since the 164 specific genes shared by the early and late phases are expressed differently during the cancer process, it is an important feature to monitor the progression of the disease. Therefore, we take advantage of the 164 genes of the intersection to perform a cluster analysis on the three groups. By using the Pearson correlation coefficient, we construct a correlation matrix, and for the clustering method, we use hierarchical clustering. Unsupervised clustering analysis is performed on all the samples to examine the discriminant effects of these genes on different phase samples. The results are shown in Fig. 2a.
We can intuitively observe that almost all of the cancer samples are clustered together, and the control samples are clustered together. Therefore, we can conclude that at the 164 genes' molecular level, the control and cancer samples have remarkable differences. However, early-and late-phase cancer samples are mixed together and are not easy to distinguish. It is also found that the expression pattern of the gene changes significantly from the normal state to the tumour state. The genes that are highly expressed in the control group have a low expression Fig. 1 The distribution of differentially expressed genes in the early and late phase. a The horizontal axis is logFC, and the vertical axis is the P-value after the negative logarithm conversion. The left panel shows the distribution of early-phase lung cancer related genes, and the right panel shows the distribution of late-phase lung cancer related genes. Up-and down-regulated genes are marked in red and green, respectively. b Up/down-regulated genes in the early and late phases are marked with four colour markers. c Each colour block corresponds to the correlation coefficient of two genes: the red represents a positive correlation, and the blue represents a negative correlation. The heatmap matrix from left to right corresponds to the control and the early and late phases level in lung cancer samples, and vice versa. At the same time, we also found that a certain proportion of lung adenocarcinoma samples and other cancer samples are not clustered together but are mixed in the control sample. This suggests that some samples, although defined as being from lung cancer patients, have a molecular level that is close to the healthy control sample and that these samples are likely to have a better prognosis.

Specific and non-specific co-expression of network analysis
We continue to construct specific and non-specific networks in different phases. In networks, the gene is a node, and the correlation is an edge; if the two genes are positively correlated, then the edge is red, and if the correlation is negative, the edge becomes green. The Cytoscape software realizes the network construction, and the network topology is analysed by using the network analysis plug-in. The edge between nodes represents the correlation coefficient; the stronger the correlation is, the thicker the edge. It is clear that in each phase-specific network, some genes assemble into clusters. There is a significant co-expression correlation between the genes within each cluster, indicating that the genes within the cluster may be functionally consistent. We use the network analysis plug-into do the topology analysis for the four networks separately, and the results are shown in Fig. 2b.
The average shortest path measures the average state of the shortest path of a node to the other nodes in the network. Therefore, the shorter the average shortest path is, the more convergent the network and the higher the signal transmission efficiency. The degree distribution is a  Fig. 2 The expression and function property of the genes. a The horizontal axis represents the sample; the vertical axis represents the genes. We use three colours to mark the different levels of classification of lung adenocarcinoma: blue for the control group, green for the early phase, and red for the late phase. The red and green blocks represent the expression pattern of the gene; green represents low expression, and red represents high expression. b shows the analysis of four topological network properties, including the average shortest path, node degree distribution, closeness centrality, and clustering coefficient. Different groups are marked in different colours. c A-C corresponds to the pathway enrichment results of three phases. The horizontal axis is the pathway term, and the vertical axis is the P-value after the negative logarithmic transformation. We label the number of genes in the pathway by dark blue and light blue markers. The brighter the colour, the more genes that are enriched in the pathway; the darker the colour, the fewer genes that are enriched Fan et al. J Transl Med (2018) 16:205 measure of the number of adjacent nodes in the network. The higher the degree, the higher the number of adjacent nodes that can affect the gene and the higher the signal transmission efficiency. Closeness centrality reflects the proximity of a node to other nodes in the network; the closer to the centre, the stronger the network shrinkage is and the closer the distance is between genes. The clustering coefficient is a sub-module that represents the ability of adjacent nodes to form a complete graph in a graph, a high clustering coefficient, and a sub-module that may exist in the network. It can be seen from Fig. 2b that the changes in the early phase of lung adenocarcinoma in three specific networks are the most obvious, which is reflected in the average shortest path, the increase in the degree, the clustering coefficient and the closeness centrality. This series of network topological features suggests that in patients with early lung cancer signal stimulation, the patient's body produces a significant stress response to resist and compensate for abnormal molecular function. In the early stage of lung cancer, the specific network is obviously contracted, thus reducing the average shortest path, increasing the clustering coefficient and approaching the centre, and improving the performance of the network signal transmission as a whole. However, with the further progress of lung cancer, this stress response in the early phase disappears or is insufficient to compensate for functional abnormalities and reflected in a gradually decreased network performance. The average performance of a non-specific network is between the normal and disease states, and because of the greater number of nodes and edges in nonspecific networks, the degree distributions are relatively more discrete. Detailed network structures and genes information of specific or non-specific networks can be seen from Fig 3.

Functional pathway enrichment analysis
We apply a functional enrichment analysis to each stagespecific gene set. For the analytical method, we adopt the Fisher method, and the significance threshold is P < 0.05. We obtain the significance (P) of the path involved in each stage and the number of genes enriched in the corresponding pathway. The results are shown in Fig. 2c.
The three graphs in Fig. 2c correspond to early phase genes, late phase genes and their intersections. By observation, the significant functions involved in the early phase include primary immunodeficiency and RNA transport. The specific functions of the late phase concentrate in the T/B cell receptor signalling pathway, chemokine signalling pathway, NF-kappa B signalling pathway and primary immunodeficiency. Functions related to intersection genes concentrate on primary immunodeficiency. Functional enrichment analysis suggests that the mechanism of immune regulation in the process of lung adenocarcinoma significantly changes. Additionally, abnormal immune systems include innate immunity, T/B lymphocyte-regulated specific adaptive immunity, natural killer-regulated nonspecific immunity, and other infections and inflammation-related functions. This further suggests that abnormalities of the immune system are essential causes of the progression of lung adenocarcinoma.
Moreover, combined with the network topology in the early and late phase changes in lung cancer, we speculated that during the early phase of lung cancer, the collective immune system in the stress response process plays an important role. Abnormalities occurred in the functions related with the initiation of innate immunity and adaptive immunity in the early phase, and therefore, the system network performance has a brief increase. However, with the disease progression, the immune system is insufficient to compensate for abnormal function. Therefore, the functional level of the immune system is an important factor in determining the risk of lung adenocarcinoma and an important indicator of the early diagnosis of lung adenocarcinoma.

Functional pathway imbalance score
We use Eq. 3 to calculate the variance score for each enriched functional pathway. To investigate whether the abnormality of these functional levels is significantly present in different phases, we analyse the imbalance scores in each of the three groups of samples. Finally, we select 12 significant pathways in the ANOVA analysis, as shown in Table 2.
We can see that these 12 pathways have obvious differences in the three phases; the P value is less than 0.05. To more intuitively analyse the imbalance of each pathway in these three phases' samples, we use the scatter plot to visualize the dynamic process of 12 pathways, as shown in Fig. 4.
It can be seen that the functional levels in the normal phase are inconspicuous and fluctuate near 0; the functional levels from the early phase begin to produce a significant fluctuation. To further clarify the direction of imbalance and variation of each pathway in the three phases, we compare the average distribution of each pathway through the boxplot visualization, as shown in Fig. 5.
It can be intuitively seen that in most of the pathways, the scores of the three groups are gradually reduced. This suggests that these levels of function are inhibited during the progression of lung cancer, and the level of immune system function is reduced by cancer invasion. The functional level of the actin cytoskeleton pathway presents an upward trend in the early phase, and then, with the progression of cancer, the functional level further declines. This may be caused by the stress response of the body under the stimulation of the cancer. Therefore, functional abnormalities can be found to occur in the early phase of lung cancer, and thus, the use of genes regulating these functions as a diagnostic feature can achieve the early diagnosis of lung cancer.

Early-phase diagnostic marker screening using the RFE algorithm
We found that the functional level of variation occurs in the pathways enriched from early-and late-phase specific genes, compared with the normal state. To achieve the early diagnosis of lung cancer and identify the genetic markers that significantly change with the progression of lung cancer, we used the recursive feature elimination The specific and nonspecific coexpression networks. The specific and non-specific networks of different phases, from a to d, corresponding to control, early-, and late-phase related genes and the overlapping gene set. The closer the node colour is to the blue, the higher degree the node has in the network; the closer the node colour is to red, and the lower degree the node has algorithm on the 198 genes. Finally, 12 genes were screened as diagnostic markers. Figure 6a shows the RFE algorithm optimization process; when the number of feature genes is 12, the model has the highest precision. Therefore, we use these 12 genes as the features to train the diagnostic model.

Risk diagnostic model construction based on early screening genes
We take 12 genes as features, using an SVM (support vector machine) to construct the classification model. During the initialization process, we set all the model parameters as the default parameters and intensively test the initial accuracy of the model in the training set. The parameter optimization process utilizes the grid search algorithm and iterates to find the optimal combination of parameters. The final model classification prediction results are shown in Fig. 6b. Figure 6b illustrates the classification efficiency of the model by ROC curve evaluation. In the process of crossvalidation, we randomize samples each time, taking four to do the training. For a prediction graph, the horizontal axis is the false positive rate; the vertical axis is a true positive rate, and the final model achieves an average accuracy of 0.91. The five-time cross validation accuracy is similar to the accuracy of the model in the training set, indicating that the model does not exhibit significant over-fitting. Using our predictive model of training, we can achieve early prediction of lung adenocarcinoma and distinguish benign and malignant progression. Notably, it not only provides new insight into the pathogenesis of lung adenocarcinoma but also suggests a new diagnostic marker or therapeutic target.

Survival analysis for validation
We downloaded samples of lung adenocarcinoma from the TCGA database as independent data for survival analysis. We identified 12 marker genes that could predict the diagnosis of early lung cancer patients, suggesting that these genes were significantly altered by disease signal stimulation in the early stage of cancer. After the functional analysis, we found that these genes were involved in the inherent immune and adaptive immune system; therefore, we speculate that with the progression of lung cancer, the patient's immune system followed with abnormal functional levels. To investigate whether these genes affected the survival of the patients by interfering with the immune level of the patients, we set the samples with variations of these 12 genes to be considered as high-risk groups, and the samples with no significant variations were set as the low-risk group and combined with Cox regression. The survival analysis results are shown in Fig. 6c. The significance (P) value was 0.03, indicating that the two groups of samples present significant differences in the survival level.

Discussion
X-ray is so far the primary tool for the clinical diagnosis of lung cancer, but for early lung cancers, misdiagnosis is unavoidable because of its low sensitivity. However, another method to confirm the diagnosis is tissue puncture, which is also a diagnostic tool that may cause cancer cells to spread and metastasize. Lung adenocarcinoma in the early stage is under restriction in local lesions by using ultrasound or X-ray, and other detective methods may lead to a missed diagnosis [25]. However, in the early stages of cancer development, the level of molecular expression has changed. This change occurs because the microenvironment of cancer tissue lesions is usually accompanied by the local inflammatory response, which activates the body's stress response and immune response. In the initial stage, under the regulation of innate immunity and adaptive immunity, multiple genes are differentially expressed to compensate for function and resist cancer [26]. However, the low sensitivity and specificity impede the broad application of these biomarkers in clinics. Therefore, to develop new methods and novel diagnostic biomarkers is necessary for the early detection of lung adenocarcinoma.
In the late phase of lung cancer, malignancy is increased due to the impact of decompensation. Therefore, more genes tend to present diverse expression. At the same time, due to the different stages of pathogenesis and degrees of malignancy, the gene expression patterns are remarkably different [27]. Hence, precisely distinguishing early-and late-phase patients from gene expression patterns as well as the level of functional  abnormalities is a milestone for the early diagnosis of lung adenocarcinoma. In this study, we identified early-and late-phase specific genes and the shared genes that are always differentially expressed across the cancer process. The flow chart of the whole study was shown in Fig. 7. Combined with unsupervised clustering analysis, we found that there were significant differences in the expression patterns of 164 shared genes between early-phase lung cancer (stage I) and advanced lung cancer (stage II-IV). To further discover the specificity between the two phases, we separately analysed each phase-specific gene set. We used the Pearson correlation coefficient to calculate the similarity between any two genes in each phase-specific gene set [13]. By using the relevant gene relationship, we constructed the co-expression network and analysed the topology of the network. By comparing the specific characteristics of the co-expression networks, we found that the network structure changes when the degree of malignancy increases [28].
With the rise in the degree of malignancy of lung cancer, the average shortest path of the system network gradually increases, suggesting that the network is becoming looser, and the signal transmission efficiency is gradually reducing [29]. Similarly, the closeness centrality reduces as the progression of the disease also prompts the network performance to decline [30]. Changes in the network structure suggest that in the process of lung cancer progression, the interactions between the genes exhibit dynamic changes under the regulation of the immune The red curve is the initial model accuracy without optimization. The green curve is the accuracy of the model after feature selection and parameter optimization. The blue curve is the average accuracy calculated using a fivefold cross-validation method. c The survival analysis. The horizontal axis is the survival time (month). The vertical axis is the percentage of patients in different groups


extfieldsofstudy: Medicine Geography

text:
Comments on the genus Diplura C. L. Koch, 1850, with description of two new species (Araneae, Mygalomorphae, Dipluridae)

Abstract Two new species of Diplura C. L. Koch 1850 are described from Brazil: Diplura mapinguari sp. n., from the state of Rondônia in southeastern Amazonia, northern Brazil, and Diplura rodrigoi sp. n., known from southeastern and central west regions of Brazil. Diplura rodrigoi sp. n. is morphologically similar to D. lineata (Lucas, 1857), D. sanguinea (F. O. Pickard-Cambridge, 1896), and D. mapinguari sp. n. Comments on diagnostic characters of Diplura are included. The synonymy of D. maculata (Mello-Leitão, 1927) with D. catharinensis (Mello-Leitão, 1923) is corroborated. A classification of color pattern of the dorsum of the abdomen is given.

Introduction
Diplura C. L. Koch, 1850 is a Neotropical mygalomorph genus currently including 17 species distributed from Panama to Argentina World Spider Catalog 2018). Most species were described from southeastern and southern Brazil, with additional records from Panamá and numerous South American countries: Argentina, Bolivia, Colombia, Ecuador, Paraguay, and Venezuela . Recently, we have examined many specimens of several species from Peru. The genus is currently recognized by the following combination of characters : simple lyra, formed by a single series of strong and thickened setae, with tip spatulate and not curved (Maréchal and Marty 1998, Drolshagen and Bäckstam 2011, Pedroso and Baptista 2014, legs with thin and restricted scopula (Drolshagen andBäckstam 2011, Pedroso andBaptista 2014), tarsi with only a few cracks (Drolshagen andBäckstam 2011, Pedroso andBaptista 2014) and males with short and thickened palp tibia (Pedroso and Baptista 2014).
Herein we describe two new species of Diplura, one from the state of Rondônia, southeastern Amazonia, northern Brazil, and the other from southeastern and central west Brazil, both based on male and female specimens.

Materials and methods
The description of color pattern is based on specimens preserved in 75% ethanol. Information and photos of living animals were included, when available. Observations, photographs and measurements were made with a Leica DFC295 camera attached to a Leica M205C stereoscopic microscope. All photos were edited in the Photoshop CS5 software and plates were prepared with CorelDraw X7 software. Measurements are given in millimeters, unless otherwise noted. Scale bars represent 1 mm, unless otherwise noted. Body length was measured from the anterior margin of the chelicerae to the posterior border of the abdomen, without spinnerets. Carapace length was measured from anterior margin of the clypeus to the posterior border. Each article of the pedipalp and legs was measured in retrolateral view, from the basal condylus to the distal one. The receptaculum seminis (spermathecae) was examined through dissection of the genital region of the females, cleaned, and immersed in clove oil for clearing. Geographical coordinates were obtained from Geonames (2018). The distribution map was created using ESRI ARCGIS 10 software.
Abbreviations. Institutions (and curators): The specific name is taken from the folklore of Amazonian Indian tribes. "Mapinguari" is a magical creature, a huge long-haired animal, with long arms and large claws.
Diagnosis. Diplura mapinguari sp. n. shares the same oblique continuous light stripes on the dorsum of the abdomen with other species (D. lineata, D. sanguinea, and D. rodrigoi sp. n.), but differs by having the stripes unequally spaced. The three median stripes of D. mapinguari sp. n are wider, almost touching each other near the middle and merging with the ventral light background (Fig. 1). The male of D. mapinguari sp. n. has the longest embolus in Diplura, 2.5× longer than the bulb (Fig.  5). Females also have an elongated receptaculum seminis, with a thin stalk, and only three lobes: the internal lobe is lateral and larger; the other two are a pair at the apex (Fig. 10).
Female (paratype, IBSP 12336). Measurements: Body length 28.1; carapace length 11.3, carapace 9.2; abdomen length 12.3, abdomen width 8.7. Leg formula 4123, total length: I 37.3, II 34.7, III 35.2, IV 44.8. Females resemble males, except by the following characteristics. Carapace: Clypeus length around 1/2 the diameter of AME, with anterior margin bearing four thick setae, elongated, and turned forward. Eye tubercle with three elongated setae at anterior margin and five setae between the posterior eyes. Chelicerae with 12 promarginal teeth. Maxilae with 26 cuspules. Lyra with 13 elongated setae, slightly curved medially. Genitalia: Receptaculum seminis ( Fig. 10) paired, elongated, separated by a little less than its height, with thin and elongated stalk, bearing one large and elongated internal lobe and two distal lobes that are not as large as the basal one.
Color pattern. Carapace reddish brown, with thoracic sulci a bit darker. Eye region black. In live juveniles, the carapace presents an almost black color and the colors of all body parts are more vivid (Fig. 1). Chelicerae, labium, sternum and coxae reddish brown. Legs mostly reddish brown, with darker femurs. Abdomen black, dorsum with five light transverse stripes, of different sizes and shapes (Fig. 1). First and last stripes shorter and thinner (disappearing in some specimens). Three median stripes wider and shorter, widest at their middle region, where they almost touch one another. Compared to the last median stripe, the two first median stripes are wider, longer and are connected to the venter. The carapace is covered with abundant light brown setae.
Distribution. Known only from Porto Velho and Guajará-Mirim, Rondônia state, Amazon area, northern Brazil (Fig. 29). Etymology. This species is named after the biologist, Rodrigo de Cerqueira da Costa, who first directed Denis Pedroso and Alessandro Giupponi in their studies of zoology.
Diagnosis. Diplura rodrigoi sp. n. differs from other species displaying several oblique continuous light stripes on the dorsum of the abdomen (D. lineata, D. sanguinea, and Diplura mapinguari sp. n.) by its very characteristic stripe pattern. Its elongated oblique stripes are very thin near the median line of the dorsum but continuously widening and approaching each other towards the venter, where they merge in the light background (Fig. 25). The three median stripes are longer and not as wide and confluent at their middle portion as in Diplura mapinguari sp. n. In contrast to D. lineata and D. sanguinea, the light stripes are longer and clearly delimited, with well-defined borders. Males of D. rodrigoi sp. n. are the only Diplurinae with a prolateroventral swelling at the distal portion of metatarsus I, which is glabrous and lighter than the surrounding areas (Fig. 13). As in D. sanguinea, males of Diplura rodrigoi sp. n. bear a relatively wide sperm duct from the base of the embolus towards the apex (Fig. 17), in contrast to Diplura mapinguari sp. n. and D. lineata, where the sperm duct is strongly constricted and becomes almost filiform towards the apex. Diplura rodrigoi sp. n. males differ from D. sanguinea by having the megaseta of tibial spur of leg I a bit shorter than the basal portion of the spur (Fig. 12) and 10-11 setae in the lyra (Fig.  19), in contrast with megaseta longer than the spur and seven setae at the lyra in the latter species.
Female (Paratype, UFRJ 0102): Measurements: Body length 26.4; carapace length 9.6, carapace width 8.2; abdomen length 12.3; abdomen width 7.6. Leg formula 4123, length: I 30.4, II 24.8, III 23.7, IV 30.6. Females resemble males, except by the following characteristics. Carapace: Clypeus length similar to AME diameter, with anterior margin bearing five thick setae, elongated and turned forward. Eye tubercle with three thick and long setae and ten smaller setae, of variable size, at anterior margin. Area between posterior eyes with single long and thick setae. AME separated from each other by 80% of their diameter. Chelicera with 12-13 promarginal teeth. Plectrum with five thick and elongated setae. Maxillae with 15-17 cuspules. Lyra (Fig.  19) formed by 11-12 subequal setae, except for the first one, which is shorter and thinner. Setae slightly spatulated, with a thin and pointed apex. Genitalia: Receptaculum seminis (Fig. 20) paired, separated by a distance similar to its own height. Stalk relatively wide, its diameter a little smaller at the distal third, near the five or six apical large lobes, with all of them being similar in size.
Color pattern. Carapace reddish, with thoracic sulci a little darker. Eye region darkened. Chelicerae reddish brown. Labium, sternum and coxae orange brown. Sigillae slightly darker than sternum. Legs light brown. Dorsum of abdomen dark brown with five elongated pale brown transverse stripes, with their width diminishing posteriorly. Stripes continuous, with well-defined borders, thin and far from each other near the midline of the dorsum, but growing wide and closer toward the venter, where they blend into the pale background. When alive, they bear a vivid red carapace and dark brown legs and abdomen.
Natural history. The female paratype from Mendes (UFRJ 0102) was collected in a short tunnel, with a small silk lined opening, on the slope of an old road through second-growth Atlantic Forest. Other specimens have been collected under rocks and logs on forest floor. Most specimens were caught by pitfall traps in forested areas.
Records Distribution: Widespread and common in a wide area of southeastern and midwestern Brazil, from north Espírito Santo state and central-south Minas Gerais state, reaching west to the Distrito Federal and south, down to southeastern São Paulo state (Fig. 29). It is sympatric with D. lineata (Lucas, 1957), a similar species, in the municipalities of Casimiro de Abreu and Rio de Janeiro, both in Rio de Janeiro state, Brazil (see Pedroso et al. 2016).

Notes on the diagnosis and composition of Diplura
The diagnosis of Diplura in relation to other Diplurinae needs some comments and amendments. The presence of a lyra in the inner face of the palp maxilla is a key character that allows for separation between lyrate (Diplura C. L. Koch, 1850, Harmonicon  , 1879). The single line of stiff setae found in lyra of Diplura and Harmonicon is quite different from the plate-like lyra, with several series of setae, in Trechona (Pedroso et al., 2008). In Harmonicon, the lyra is formed by four to seven setae in adult specimens, with a flattened and curved tip (Drolshagen andBäckstam 2011, Pedroso andBaptista 2014). In Diplura, the number of setae in the lyra is quite variable. The described species of Diplura have stiff setae in the lyra ranging from four to seven in Diplura taunayi (Mello-Leitão, 1923) to 15 in Diplura studiosa (Mello-Leitão, 1920). However, we found specimens of Diplura from southern Brazil with up to 18 setae. On the other hand, specimens from southeastern and northeastern Brazil may have only one to three stiff setae in the lyra. One specimen we examined had only one seta in one maxilla, but the other maxilla was devoid of stiff setae. Those findings show that the number of setae in the lyra may increase or decrease in different species, and even the complete absence of the lyra may happen in some specimens. Therefore, as pointed out by Drolshagen and Bäckstam (2011) and Pedroso and Baptista (2014), the number of stiff setae in the lyra is not diagnostic for the genera.
The shape of the stiff setae of the lyra is much more variable than previously pointed out in the literature. In many described species of Harmonicon, including the typespecies Harmonicon rufescens F. P.-Cambridge, 1896, the setae in the lyra are strongly hook-shaped, with a strong distal curvature and an acute and not widened tip (F.  fig. 6), similar to some undescribed species we have examined from Brazil. The lyra of one of those undescribed species has a very wide and clearly spatulated tip. On the other hand, there are species of Diplura with curved stiff setae (ex. D. mapinguari sp. n., Fig. 9) and the tip not clearly spatulated (see D. lineata, Pedroso et al. 2016: fig. 22-24, and D. rodrigoi sp. n., Fig. 19). Therefore, we consider that the variation found on the shape of the stiff setae of the lyra in both Diplura and Harmonicon are too wide and clearly overlapping.
Besides the lyra, there are also other diagnostic characters in the maxilla. On its ventral face, there is a transversal suture just below the lyra. In Diplura, the area just below the basal portion of the maxillar suture is either glabrous or has just some small thin setae (Figs. 9, 19). On the other hand, there is a field of spiniform setae in that area in Harmonicon (Pedroso and Baptista 2014: fig. 6) and Trechona (Guadanucci et al. 2016: fig. 1g). Additionally, there is also a dense fringe of long and thin setae between the lyra and the basal portion of the suture in Trechona (Guadanucci et al. 2016: fig. 1g).
Specimens of Diplura have relatively short and thin scopula, sometimes almost absent in the metatarsus and tarsus of the legs. Scopulae in the other lyrate Diplurinae are formed by long numerous setae and cover most of the distal articles of the legs. In Diplura, the tarsi of the legs are not pseudosegmented and not very flexible, lacking large and numerous cracks; they have only a few ventral, thin cracks, that allow just a limited bending of the article (Drolshagen andBäckstam 2011, Pedroso andBaptista 2014).
An easy way to separate males of Diplura from other Diplurinae is by the short and incrassate pedipalp tibia (Fig. 21), which contrasts strongly with the elongated and thin tibia of other Diplurinae (Fig. 22). The ratio length/width of the pedipalp tibia reaches 3.6× in Diplura and is at least 5× in the other Diplurinae. This character was already pointed out by Pedroso and Baptista (2014) and Pedroso et al. (2016).
The receptaculum seminis of most Diplura species is formed by a short and relatively wide stalk topped by a cluster of many globular to ovoid lobes (fundus), of variable sizes (Fig. 20, Pedroso et al. 2016: fig. 13, 25). In D. mapinguari sp. n., the stalk is elongated and thinner and the lobes are separated in two groups (Fig. 10). In Harmonicon (Drolshagen and Bäckstam 2011: fig. 5) and Trechona Baptista 2004: fig. 1, Pedroso et al. 2008: fig. 4), there are usually two branches: the first one simple and at the inner side of the stalk, and the second one distal and topped by one large fundus, sometimes with lateral flaps. However, there are no clear multiple lobes as in Diplura.
Currently, 17 species of Diplura are considered as valid, distributed from Panama to Argentina World Spider, Catalog 2018). However, Diplura maculata (Mello-Leitão, 1937) has been already synonymized with Diplura catharinensis (Mello-Leitão, 1927) by Bücherl (1962: 261). Although this synonymy had been completely overlooked in the literature and catalogs, we agree with Bücherl after examining the types and additional material from the type-localities or nearby areas in Santa Catarina state, southern Brazil. As a result, Diplura now includes 18 species, taking into consideration the synonymy referred to above and the two new Brazilian species herein described.
Most Diplurinae species have a pattern of light stripes or dots over the dark dorsum of the abdomen. The exceptions are Harmonicon, and some species of Diplura, whose abdomen are uniformly dark, without a contrasting light pattern. The light pattern of most Diplurinae may be present only at the sides of the dorsum, but usually extend either to the middle area of the dorsum or to the venter, even merging with the background color in that area. The species of Diplura may be grouped below in six color pattern types following the different markings on the dorsum of abdomen, numbered from I to VI (Figs 23-28). These pattern types are an easy way to group species and also allows a clear distinction between these groups or isolated species.
Type IV (Fig. 26): Pattern formed by elongated transverse light blotches, sometimes more or less restricted to the sides of the dorsum: D. garbei (Mello-Leitão, 1923), D. petrunkevitchi (Caporiacco, 1955), D. sanguinea (F. O. Pickard-Cambridge, 1896). Type V (Fig. 27): Pattern formed by irregular large light spots, with most of them connected to one another, forming an irregular longitudinal stripe on the sides of the dorsum: D. lineata (Lucas, 1857). Type VI (Fig. 28): Pattern formed by three irregular very wide light spots, widest at their middle region, where they almost touch one another: D. mapinguari sp. n.
The only described species of Diplura not listed above is D. garleppi (Simon, 1892), from San Mateo, Bolivia. As its holotype has not been found in the Museum national d´Histoire naturelle, Paris (E. Leguin, pers. comm.), the short description does not give enough detail on the color pattern and we have not examined any specimen ascribable to D. garleppi, we were not able to include it in any of the proposed color pattern types.


Structural variants exhibit allelic heterogeneity and shape variation in complex traits

Despite extensive effort to reveal the genetic basis of complex phenotypic variation, studies typically explain only a fraction of trait heritability. It has been hypothesized that individually rare hidden structural variants (SVs) could account for a significant fraction of variation in complex traits. To investigate this hypothesis, we assembled 14 Drosophila melanogaster genomes and systematically identified more than 20,000 euchromatic SVs, of which ∼40% are invisible to high specificity short read genotyping approaches. SVs are common in Drosophila genes, with almost one third of diploid individuals harboring an SV in genes larger than 5kb, and nearly a quarter harboring multiple SVs in genes larger than 10kb. We show that SV alleles are rarer than amino acid polymorphisms, implying that they are more strongly deleterious. A number of functionally important genes harbor previously hidden structural variants that likely affect complex phenotypes (e.g., Cyp6g1, Drsl5, Cyp28d1&2, InR, and Gss1&2).Furthermore, SVs are overrepresented in quantitative trait locus candidate genes from eight Drosophila Synthetic Population Resource (DSPR) mapping experiments. We conclude that SVs are pervasive in genomes, are frequently present as heterogeneous allelic series, and can act as rare alleles of large effect.

Introduction
Understanding the molecular basis of heritable variation in complex traits is of central importance to evolution, animal and plant breeding, and medical genetics (Mauricio 2001, Goddard and Hayes 2009, Stranger et al. 2011. Over the last decade, short read (50-150 bp) genomic data appropriate for characterizing SNPs and small indels in non-repetitive genomic regions has accumulated at an exponential rate (Shendure andJi 2008, Bansal et al. 2010). This in turn has catalyzed hundreds of quantitative trait loci (QTL) mapping and genome wide association (GWAS) studies in model organisms, humans, and agriculturally important animals and plants (Varshney et al. 2009, Davey et al. 2011, Day-Williams and Zeggini 2011. Despite these efforts, for most traits, GWAS hits only explain a small fraction of known trait heritability (Manolio et al. 2009, Eichler et al. 2010. One hypothesis accounting for hidden genetic variation is that individually rare large mutations (>100 bp) that alter genome structure make significant contributions to complex trait variation (Frazer et al. 2009, Eichler et al. 2010). These structural variants (SVs) change the genome via duplication, deletion, transposition, and inversion of sequences. This hypothesis is attractive since rare causative variants are difficult to detect with GWAS (Spencer et al. 2009). Moreover, genotyping approaches based on short reads or microarrays miss a significant number of SVs Eichler 2016, Chakraborty et al. 2018). Finally, SVs are likely to be more deleterious and deleterious more often than SNPs (Emerson et al. 2008, Conrad et al. 2010, Cridland et al. 2013, Rogers et al. 2015.
Accurate and gap-free genomes provide a direct and reliable path to comprehensive identification of SVs (Alkan et al. 2011, Huddleston et al. 2017. To achieve this goal, we assembled reference-quality genomes for fourteen geographically diverse Drosophila melanogaster strains (Fig. 1a) using Single Molecule Real Time sequencing (Chakraborty et al. 2016). These assemblies are contiguous and complete BUSCO 22 99.9-100%)( Table 1, Fig. 1b, supplementary   Table 1), making them comparable to the D. melanogaster reference genome, arguably the best metazoan genome assembly. Thirteen of the fourteen strains are near isogenic founders of the Drosophila Synthetic Population Resources (DSPR) (King et al. 2012), a large set of advanced intercross recombinant inbred lines (RILs) designed to map quantitative trait loci (QTLs) . We also assembled the genome of Oregon-R, an outbred stock widely used as a "wild-type" strain both by Drosophila geneticists and by large scale community projects like modENCODE (The modENCODE Consortium et al. 2010, Graveley et al. 2011, Schwartz et al. 2012).

DNA extraction
Genomic DNA was extracted from females following the protocols described in Chakraborty et al. 2016 and the genomic DNA was sheared using 10 plunges of a 21gauge needle, followed by 10 of a 24-gauge needle (Jensen Global, Santa Barbara).
SMRTbell template library was prepared following the manufacturer's guidelines and sequenced using P6-C4 chemistry in Pacific Biosciences RSII platform at University of California Irvine Genomics High Throughput Facility. The total number of SMRTcell and base pairs sequenced, and read length metrics for each strain is given in supplementary Table 6.

Genome assembly
The genomes were assembled following the approach described in Chakraborty et al. (2016). For all calculations of sequence coverage, a genome size of 130Mbp is assumed (G =130×10 6 bp). For each strain, we generated a hybrid assembly with DBG2OLC (Ye et al. 2016) and longest 30X PacBio reads and a PacBio assembly with canu v1.3 (Koren et al. 2017) (supplementary Table 5). The paired end Illumina reads were obtained from King et al. (King et al. 2012). The hybrid assemblies were merged with the PacBio only assemblies with quickmerge v0.2 (Chakraborty et al. 2016, Solares et al. 2018) (l = 2Mb, ml = 20000, hco = 5.0, c = 1.5), with the hybrid assembly being used as the query. Because the PacBio assembly sizes were closer to the genome size of D. melanogaster, we added the contigs that were present only in the PacBio only assembly but not the hybrid assembly by performing a second round of quickmerge (Solares et al. 2018). For the second round of quickmerge (l = 5Mb, ml = 20000, hco = 5.0, c = 1.5), the PacBio assembly was used as the query and the merged assembly from the first merging round the reference assembly. The resulting merged assembly was processed with finisherSC to remove the redundant sequences and additional gap filling using raw reads (Lam et al. 2015). The assemblies were then polished twice with quiver (SMRTanalysis v2.3.0p5) and once with Pilon v1.16 (Walker et al. 2014). For Pilon, we used the same Illumina reads as used for hybrid assemblies.

Comparative Scaffolding
We scaffolded the contigs for each assembly based on the scaffolds from the reference assembly (Hoskins et al. 2015), following a previously described approach . Briefly, TEs and repeats in the assemblies were masked using RepeatMasker (v4.0.7) and aligned to the repeat-masked chromosome arms (X, 2L, 2R, 3L, 3R, and 4) of the D. melanogaster ISO1 assembly using MUMmer (Kurtz et al. 2004). After filtering of the alignments due to the repeats (delta-filter -1 ), contigs were assigned to specific chromosome arms on the basis of the mutually best alignment. The scaffolded contigs were joined by 100 Ns, a convention representing assembly gaps.
The unscaffolded sequences were named with a 'U' prefix.

BUSCO analysis
We ran BUSCO (v3.02) (Waterhouse et al. 2017) on the Pilon polished pre-scaffolding assemblies to evaluate the completeness of all the assemblies relative to the ISO1 release 6 ( r6.13 ) assembly. We used both the arthropoda and diptera datasets for the BUSCO evaluation. For the arhtropoda database, three orthologs (EOG090X0BNZ, EOG090X0M0J, EOG090X049L) weren't found in any of the 15 strains (ISO1, Oregon-R, and 13 DSPR founders). Further inspection of these orthologs revealed that they are present in ISO1 even though the BUSCO analysis misses them when applied to ISO1 (EOG090X0BNZ is CG3223, EOG090X0M0J is Pa1, and EOG090X049L is CG40178).
Consequently, we removed these three genes from consideration as uninformative.

Variant detection
For variant detection, we aligned each DSPR assembly individually to the ISO1 release 6 assembly (release 6.13) (Hoskins et al. 2015) using nucmer (nucmer -maxmatch -noextend) (Kurtz et al. 2004). We identified and classified the variants using SVMU 0.2beta (Structural variants from MUMmer) (n = 10) . SVMU classifies the structural differences between two assemblies as insertion, deletion, duplication, and inversion based on whether the DSPR assemblies have longer, shorter, more copy, or inverted sequence, respectively, with respect to the reference genome.
The variant calls for individual genomes were combined using bedtools merge (Quinlan 2014) and converted into a vcf file using a custom script (https://github.com/mahulchak/dspr-asm). TE insertions were identified by examining the overlap between RepeatMasker identified TEs and SVMU insertion calls using bedtools, requiring that at least 90% of RepeatMasker TE annotation overlap with svmu insertion annotation. 12.8% SV mutations, for which mutation annotation were complicated by secondary mutations, were flagged as 'complex' (CE=2 in the VCF file).
Additionally, 16.3% SVs that were located within 5Kb of a complex SV were often part of a complex event and were also assigned a tag (CE=1) to differentiate them from the unambiguously annotated SVs (CE=0).

Genotype validation
To determine the genotyping error rate, a set of randomly selected 50 simple (CE=0) SVs obtained from SVMU were manually inspected on UCSC genome browser representation of the multiple genome alignment of the 15 genomes (http://goo.gl/LLpoNH). Furthermore, to estimate the genotyping accuracy of the SVs occurring in the vicinity of the complex mutations, where mutation annotation is complicated by alignment ambiguities, we manually inspected 217 SVs occurring within 20Kb of 50 randomly selected complex (CE=2) SVs. Among these, 3/217 and 0/50 SVs were absent in the UCSC browser and therefore they are likely mis-annotated by our pipeline. The mis-annotated SVs (insertion in A1 and tandem array CNV in A7) are located in a complex, repetitive, structurally variable genomic region on chromsome 3L (3L:7669500-7679100) (supplementary Fig.5).

Comparing SV genotypes from de novo assemblies to short read only calls
TE genotypes for the founders (Cridland et al. 2013) were downloaded from flyrils.org and the insertion coordinates were lifted over to the current release (release 6) of the reference genome (Hoskins et al. 2015) using UCSC liftover tool (Kent et al. 2002). For detection of the duplicates, we have previously found that discordant read pair based method (Pecnv) (Rogers et al. 2014) was comparable to split read mapping (Ye et al. 2009) and more reliable than methods based on coverage alone (Abyzov et al. 2011), so we used Pecnv. Pecnv was run using the settings described before . Because svmu reports tandem duplicate CNVs as insertions (with appropriate CNV tags to separate from TE and other insertions) and Pecnv reports sequence range being duplicated, the SVMU CNV insertion coordinates were extended by 100 bp before comparison (bedtools intersect) between Pecnv output and svmu output was conducted. The non-TE indel genotypes were obtained from Pindel output (the "LI" and "D" events) using the commands described previously 15 . For determining population frequency of indel SVs (e.g. the reference FB element in InR), Pindel output based on the alignment bam files were used. We only estimate the false negative rate of short read only callers, but note that these methods also generate false positive SV calls.

Gene expression analysis
The preprocessed expression data for female heads ) and IIS/TOR expression data (Stanley et al. 2017) from whole bodies were downloaded from www.flyrils.org. Expression QTL analysis (supplementary Fig. 11) for Cyp28d1 and Gss1 using the head expression data were performed using the R package DSPRqtl following the instructions provided in the manual (DSPRscan,model = gene ~ 1,design = "ABcross"). When expression data for multiple isoforms were present, expression data only for the longest transcript that is expressed in the head was used. The genotype values at the eQTL were determined using the function DSPRpeaks included in the DSPRqtl package. No eQTL were found for InR so the genotype values for the InR expression data were obtained by assigning the founder genotypes to the RILs used in the IIS/TOR expression dataset, using the posterior probabilities of the forwardbackward decoding of the HMM for the panel B RILs available on www.flyrils.org. Drsl5 expression levels in A4 and A3 were obtained from a publicly available RNAseq dataset (Marriage et al. 2014).

Comparison of site frequency spectra
The histogram of allele frequencies (site frequency spectrum or SFS) was collated for four categories: synonymous SNPs, non-synonymous SNPs, duplicate CNVs, and TE insertions. The frequencies of SNPs were collected from the VCF file (King et al. 2012) using vcftools and bcftools (Danecek et al. 2011, Li 2011. The frequencies of SVs were collected from the column 4 of the combined SVMU output for the TE insertions and duplication CNVs from all DSPR strains (https://github.com/mahulchak/dspr-asm).
Complex mutations (CE=1 and CE=2) were excluded from the analysis. Let N be the sample size and xi be the number of sites in frequency class i, where 0 < i < N. The SFS was "folded", meaning we focused attention on the minor allele frequency (MAF), or yi = minimum(xi, N-xi). Pairwise comparisons between different SFS site categories were conducted using the χ 2 test on allele frequencies and site categories. For allele frequencies, two types of classifications were used: 1) every yi for 0 < i < N (N-1 df); and 2) considering singletons versus the other frequency categories, or yi for i = 1 versus 2

Candidate genes associated with mapped QTL
The candidate genes from DSPR QTL papers were selected based the following criteria: 1) The gene falls within the QTL peak; 2) additional functional data is cited by the authors of the respective study to highlight the gene; 3) the functional information cited by the authors did not use knowledge about structural variation affecting the candidate locus (supplementary Table 3). The additional data can either be expression data collected by the authors or existing functional data known about the genes. Only 44 candidate genes from 8 studies fulfilled these criteria but 3 among these fell outside the euchromatic boundaries used here (supplementary Table 1). Hence only 41 candidate genes were included in the SV enrichment analysis. Of the 41 candidate genes identified, 10 of them were at a single locus (GstE1-10). As a result, we carry out our analysis treating GstE1-10 as either a single gene or ten different genes (the qualitative outcome is unchanged). To test if candidate genes are longer than average genes, we considered all genes (supplementary table 4) as well as the dataset excluding the GstE1-10 genes (supplementary Table 4). The lengths of candidate genes were compared against the rest of the genome using a Mann-Whitney U test.

Candidate gene enrichment analysis
To determine if candidate genes are enriched for SVs relative to the rest of the genome, we analyzed the dataset both without merging and with merging the GstE1-10 into a single 13kb SV-burdened locus (supplementary Table 4). A Fisher's Exact Test was applied to the counts in categories of candidate gene vs. rest of the genome and SVburdened vs. SV-free genes. To account for the lengths of the candidate genes being longer than the rest of the genome, we performed a Monte Carlo resampling of the whole genome according to the histogram of gene sizes in the candidate gene lists (supplementary Table 4). We sampled from the genome by drawing from each gene length bin with a hypergeometric distribution, where n is the number of candidate genes in the candidate bin, K is the number of SV-burdened genes in the genome bin, and N-K is the number of SV-free genes in the genome bin (supplementary Table 4). We then tallied up the number of observed SVs. We repeated this 100,000 times to construct a Monte Carlo distribution of the SV burden expected of genes matching the size distribution observed in the actual candidate genes. This led to simulated size distributions that matched the observed size distributions (every Mann-Whitney U pvalue of Monte Carlo sample lengths compared against the observed candidate lengths > 0.1).

Calculating the SV burden in genes in diploid individuals
In order to calculate the distribution of SV burden expected in diploids, the haploid genotypes of each founder was paired with every other founder, for a total of 78 possible pairings. For each of these diploid pairings, the number of unique SV mutations for each gene in the genome was recorded. A mutation is said to affect a gene if it falls within the gene span, which is defined as affecting nucleotides between the start and end coordinates of the gene feature in the Drosophila melanogaster release 6.16 gff file (dos Santos et al. 2015). The number of SV mutations overlapping a gene in a given diploid combination is considered that gene's multiplicity for that combination. Any gene with a multiplicity ≥ 1 for a particular diploid comparison is considered SV-burdened for that diploid.

Results and Discussion
De novo assembly reveals a large number of previously hidden functionally important SVs: Our assemblies are extremely contiguous, with the majority of each chromosome arm represented by a single contig (Fig. 1b). We also close the two remaining gaps in the major chromosome arms of the euchromatic D. melanogaster reference genome (Chang and Larracuente 2018) in all our assemblies (supplementary Fig.1-3). We identified SVs by comparing each assembly to the reference ISO1 genome , focusing our attention on large (>100bp) euchromatic SVs (supplementary Table 2), and ignoring heterochromatin regions as they are gene poor (Smith et al. 2007) and require specialized approaches and extensive validation (Khost et al. 2017). Manual inspection of 267 randomly sampled SVs indicate that misannotations are rare (3/267), and occur in ambiguously aligned structurally complex genomic regions (supplementary Fig. 5; see Methods). We discovered 7,347 TE insertions, 1,178 duplication CNVs, 4347 indels, and 62 inversions in the 94.5 Mb of euchromatin spanning the five major chromosome arms across the DSPR founders ( Fig. 1c-d). Each founder strain exhibits 637 TE insertions, 134 duplications, 694 indels, and 7 inversions on average (Table 2).
A large fraction of the SVs (36% of non-reference TEs, 26% of deletions, 48% of insertions, 60% of duplication CNVs) present in the assemblies eluded detection using high specificity SV genotyping methods ) employing high coverage paired end Illumina reads (supplementary Fig. 6), Some of these novel events are likely to affect phenotypes. For example, extensive evidence links complex SV alleles of the cytochrome P450 gene Cyp6g1 to varying levels of DDT resistance (Daborn et al. 2002. Despite extensive study of this locus, we discovered three new SV alleles involving TE insertions that likely have different functional consequences (supplementary Fig. 7a-b). Similarly, we discovered a previously hidden tandem duplication of the antifungal, innate immunity gene Drsl5 (Yang et al. 2006) that exhibits >1000 fold higher expression relative to its single copy counterpart in line A4(supplementary Fig. 8a-b). Read pair orientation methods failed to detect this mutation because one allele bears a 5kb spacer sequence derived from the first exon and intron of Kst inserted between the gene copies (supplementary Fig. 8a).
Another duplicate allele of Drsl5 contains a Tirant LTR retrotransposon inserted into the same spacer sequence (supplementary Fig. 8a). We also easily detect the two SV mutations underlying the D. melanogaster recessive visible genes cinnabar (Warren et al. 1996) (cn) and speck (sp) present in the ISO1 reference genome (dos Santos et al. Fig. 9-10). In the case of sp a large insertion in the reference genome is mis-annotated as an intron. For cn a large exonic deletion is not identified as such (dos Santos et al. 2015). Both alleles are likely knock-outs.

SVs are deleterious:
Most TEs and duplicates are present in only one strain (Fig. 1e), indicating that purifying natural selection has prevented them from rising to higher frequencies. The folded site frequency spectrum (SFS) of the TEs and CNVs exhibit more rare variants than synonymous and non-synonymous SNPs ( Fig. 1e; p-value < 1e-10, χ2 test between frequency classes of SVs and non-synonymous SNPs), suggesting that SVs are on average under purifying selection (Emerson et al. 2008, Cridland et al. 2013, Rogers et al. 2015. Amongst SVs, TEs are more enriched for rare variants than duplicates, indicating that TE insertions are on average more deleterious than CNVs ( Fig. 1e; pvalue < 1e-10, χ2 test between frequency classes of TEs and CNVs). Under mutation selection balance models (Pritchard 2001, rare deleterious variants (minor allele frequency or MAF <1%) are predicted to contribute significantly to complex trait variation (Gibson 2012), yet are unlikely to be tagged by SNPs typically used in GWAS experiments (Manolio et al. 2009). Consequently, if individually rare SVs underlie complex trait variation, they will often go undetected in association studies (Spencer et al. 2009).

Functional structural variation at mapped QTL:
Segregation of multiple alleles at a causal genes (ie allelic heterogeneity) can mislead discovery of causative loci in GWAS experiments , though such genes can be readily identified in multi-parent panels (MPPs) via QTL mapping . However, mapping resolution is often poor, thwarting the identification of causative mutations bearing a variant in a candidate gene of obvious functional significance. Yet, putatively causative SVs are often hidden as they disproportionately escape detection by short read sequencing . This limitation can be solved in the DSPR and other MPPs, as de novo assemblies of the founders of the MPP allow genotypes at SVs to be imputed for the lines on which phenotypes are measured (King et al. 2012).
A nicotine resistance mapping study employing the DSPR identified differentially expressed cytochrome P450 genes Cyp28d1 and Cyp28d2 as candidate causative genes at a mapped QTL, but proposed no causative mutations (Marriage et al. 2014). A previous de novo assembly of one DSPR founder strain assembled a resistant allele possessing tandem copies of the Cyp28d1 gene separated by an Accord LTR retrotransposon fragment   (Fig. 2a; supplementary Fig. 11).
Our assemblies of the DSPR founder strains revealed a total of seven structurally distinct alleles in this region, including additional candidate resistant alleles harboring gene duplications (Fig. 2a-b). For example, the resistant strain A2 carries a tandem duplication of a 15Kb segment containing both Cyp28d genes. The expression level of Cyp28d1 in the adult female heads of RILs bearing the A2 genotype is highest among all founder genotypes measured (Fig. 2c). Consistent with this, DSPR RILs bearing the A2 genotype show the highest resistance to nicotine toxicity among the A genotype RILs (Marriage et al. 2014) (Fig. 2b). This implies that the extra copies of Cyp28d1 and/or Cyp28d2 account for the increased expression and concomitant resistance to nicotine. Similarly, the B4 allele comprises a tandem duplication of a 6 Kb segment, containing one extra copy of Cyp28d1 and a nearly complete copy of Cyp28d2 ( Fig. 2a; supplementary Fig. 11). RILs carrying the B4 genotype at the Cyp28d locus also show high resistance to nicotine, making the duplication a compelling candidate for the causative mutation. On the other hand, in two alleles, TE insertions disrupt Cyp28d gene structure and function. For instance, A1 has the same duplication as B4, but a 4.7Kb F element inserted in the 5th exon disrupts the protein coding sequence of the second Cyp28d1 copy, likely rendering the copy nonfunctional ( Fig. 2a; supplementary   Fig. 11). Consistent with the hypothesis that the duplication causes increased nicotine resistance, the A1 genotype is more susceptible to nicotine than B4 (Fig. 2b). All of these SV alleles are singletons, and thus represent a hidden allelic series composed of individually rare alleles.
SVs may also affect genes central to life history traits. Expression levels of the insulin signaling pathway genes show substantial variation in F1 hybrids between DSPR panel B RILs and the A4 founder (Stanley et al. 2017). Among these is Insulin Receptor (InR), which plays key roles in several life history traits related to lifespan and is likely a key molecular mediator of the tradeoff between reproductive success and longevity (Tatar et al. 2001, Toivonen and Partridge 2009, Paaby et al. 2014. Amino acid polymorphism in InR evolves under positive selection and some non-synonymous variants affect fecundity and stress response (Paaby et al. 2010, Paaby et al. 2014. Expression variation of InR also affects body size, lifespan, and fecundity (Brogiolo et al. 2001, Rauschenbach et al. 2015, suggesting that natural cis-regulatory variation might also be under selection. We discovered a 215 bp fragment of a DOC6 element within a 2nd intron enhancer (Wei et al. 2016) (Fig. 3b-c) of InR on the AB8 haplotype, and this allele exhibits reduced gene expression relative to reference genotypes (Fig. 3b). This mutation presumably disrupts the enhancer (supplementary Fig. 12), making it a plausible candidate for expression variation in InR. Another founder, A6, carries a 1,042 bp insertion of DMRT1A (LINE) in the 2nd intron and a 946 bp insertion of a fragment of PROTOP in the 3rd intron. Both affect known cis-regulatory elements (Wei et al. 2016) ( Fig. 3b-c). Except for A2 and A6, all strains, including ISO1, harbor an FB-NOF element (FB{}1698) inside the first intron of InR (Fig. 3a). Like many genes, the first intron of InR possess several transcription factor binding sites (TFBS), including those for factors Nejire and Caudal (Negre et al. 2011) (Fig. 3c). The FB-NOF element is inserted within this dense cluster of TFBS and active enhancer marks (Fig.3c).
Furthermore, the FB element is segregating at high frequency in the strains discussed here (13/15), a North American population (Mackay et al. 2012) (125/170), and a French population (Pool et al. 2012) (4/9), but is rare in populations derived from D.
melanogaster's ancestral range in Africa (Pool et al. 2012, Lack et al. 2015 (Cameroon: 0/10, Rwanda: 1/27, Zambia: 10/139) (Fig. 3d). This raises the possibility that the FB element is more common in temperate cosmopolitan populations, similar to a previously described adaptive amino acid variant in InR (Paaby et al. 2010). In total, InR harbors a remarkable amount of potentially functional structural diversity; Including these variants described, there are 9 TE insertions and two deletions throughout the gene, many of which impinge on candidate regulatory regions or transcribed portions of the gene (Fig.   3a,3c).
Public resources like modENCODE annotate molecular phenotypes (e.g., RNAseq, ChIPseq, DNAse1HSseq) against reference genomes which are often genetically different than the strains assayed (The modENCODE Consortium et al. 2010, Graveley et al. 2011, Negre et al. 2011, Schwartz et al. 2012. Canton-S (our DSPR founder A1) and Oregon-R are strains commonly used in phenotypic assays (The modENCODE Consortium et al. 2010, Graveley et al. 2011, Schwartz et al. 2012, and we observe SVs segregating between these two strains and the reference ( Table 2). Interpretation of functional genomics data such as RNA-seq can be misleading when gene copy number varies between strains. We explored the glutathione synthetase region (containing Gss1 and Gss2), which is just one example among hundreds in modENCODE that likely suffer from misleading annotations. A tandem duplication present in ISO1 has created two copies of Gss1 and Gss2, which are associated with toxin metabolism and linked to tolerance to arsenic (Ortiz et al. 2009) and ethanol induced oxidative stress (Logan-Garbisch et al. 2015). While this duplication segregates at high frequency in DSPR strains (9/13), it is absent in Oregon-R (Fig. 4a) and escapes detection with high specificity short-read methods. As a result, using transcript and ChIP data derived from Oregon-R (as used in modENCODE (Graveley et al. 2011, Schwartz et al. 2012) results in misleading annotations of the two copies in ISO1. Indeed, among the eight structurally distinct Gs alleles in our dataset, ISO1 is the sole representative of its allele (Fig. 4a). The two most common Gs alleles include one that contains only a single Gs gene (in four strains, including Oregon-R) and one carrying only a tandem duplication, creating the Gss1/Gss2 pair (in five strains, including Samarkand/AB8) (Fig.   4a). The remaining 6 alleles have SV genotypes represented by only a single individual in the sample. Collectively, this sample represents a haplotype network of structural variation involving 5 TE insertions, one duplication, one insertion comprising TE and simple repeats, and two non-TE indels. The single copy allele with a 5' insertion of a 14kb repetitive sequence comprising Nomad retrotransposon fragments exhibits the highest expression, followed by duplicate alleles, whereas single copy alleles and duplicate alleles with intronic TE insertions generally have the lowest expression levels (Fig. 4b).
Although hypotheses employing SVs to explain missing heritability have been proposed (Manolio et al. 2009, Sudmant et al. 2015, the systematic under-identification of SVs via short read-and microarray-based genotyping (Alkan et al. 2011) limits their explanatory power. Using our comprehensive SV map, we measured the prevalence of SVs at the candidate genes reported in 8 complex trait mapping experiments employing DSPR (supplementary Table 3). We consider only genes in mapped QTLs explicitly cited by the authors of the original QTL studies (supplementary Table 3; see Methods).
In total, we identified 31 candidate genes and a single, tandem array of 10 small genes from the same family (GstE1-10) (Kislukhin et al. 2013), which we consider as a single additional locus with structural variation. Half of these candidates (16/32) possess at least one SV in one founder strain, whereas only 23.4% (3,252/13,861) of all D. melanogaster genes harbor SVs (p = 0.001, Fisher's exact test). Although the candidates we tested (supplementary Table 3) are approximately twice as large as the genome-wide average (supplementary Table 4; p = 6.5×10-5), this enrichment of SVs at candidate genes is not merely a consequence of them being longer. SVs are enriched in 100,000 Monte Carlo samples matching the candidate gene length distribution (p = 0.021; Fig. 5a). These results persist when the GstE genes are considered individually instead of being merged (length p-value = 0.034 and enrichment p-value = 2.9×10-3; supplementary Table 4).
In order to illustrate how common SVs are in the genome we quantified the per gene SV burden per diploid D. melanogaster individual (Fig. 5b). Across the genome, SVs appear in 9.3% of genes in diploid individuals. Of those, more than a third (34.4%) involve multiple SV mutations (Fig. 5c). One or more SVs burden more than half of genes in and above the 20-35kb range (Fig 5b). Furthermore, individual genes bearing multiple SVs comprise more than a third of burdened genes between 20kb and 35kb in length and more than half of larger genes. These observations suggest that the contribution of rare SVs of large effect to complex traits could be pervasive.

Conclusion
Despite claims that a significant proportion of complex trait variation in humans, model organisms, and agriculturally important animals and plants are likely due to rare SVs of large effect (Eichler et al. 2010), systematic inquiry of this hypothesis has been impeded by genotyping approaches attuned to SNP detection (Alkan et al. 2011). As reference quality de novo assemblies of population samples for eukaryotic model systems become increasingly cost-effective, methodical evaluation of the contribution of SVs to the genetic architecture of complex traits becomes feasible. Our comprehensive map of SVs in Drosophila provides the means to systematically quantify the contribution of rare SVs to heritable complex trait variation (Fig 2a,3a,4a,5a). The value of comprehensive SV detection is underscored by the presence of SVs in ~50% of the candidate genes underlying mapped Drosophila QTL, and by the observation that a large fraction of Drosophila genes harbor multiple rare SV alleles. The genomes of humans and agriculturally important plants and animals harbor more SVs than Drosophila, and thus are likely more burdened with genic SVs.
The genetic heterogeneity hypothesis posits that a sizable fraction of human complex disease is associated with an allelic series consisting of individually rare causative mutations at several genes of large effect (McClellan and King 2010). Furthermore, models for complex traits under either stabilizing (Turelli 1984, Johnson andBarton 2005) or purifying selection (Pritchard 2001 with constant mutational input predict the existence of genes segregating several individually rare causative alleles that account for a sizable fraction of complete trait variation. We provide examples of SVs in genes of functional significance, and show that genes harboring SVs are overrepresented in a collection of QTL candidate genes. Hidden SVs are thus examples of collectively common but individually rare deleterious genetic variants predicted under the genetic heterogeneity hypothesis. Future de novo assemblies of other genomes, including humans, models, and agriculturally important species, would quantify the generality of observations from Drosophila.