Integrating massive RNA-seq data to elucidate transcriptome dynamics in Drosophila melanogaster

SraMongo: downloading NCBI databases

We developed a Python package (SraMongo, https://sramongo.readthedocs.io/en/latest/) to query, download, parse and store records from the SRA. SraMongo uses the Entrez utility API to query Entrez databases and download their XML records. Each sample is an individual sequencing library assigned with a unique SRA accession number (SRX/ERX/DRX) that is associated with one or more sequencing runs (SRR/ERR/DRR), which constitutes technical replicates. Although only 24 K samples were analyzed, SraMongo is capable of handling larger numbers of datasets. Each XML record is parsed into a local document database (MongoDB Community server) and indexed by experiment accession (SRX/ERX/DRX). We used this local database to generate a processing queue, explore technical and biological metadata and aggregate metadata for use in downstream analyses.

MassiveQC: quality control of massive high-throughput sequencing data

MassiveQC obtained the download URL of each run (SRR) from the ENA Portal API (https://www.ebi.ac.uk/ena/portal/api/filereport) and used the command ‘ascp’ or ‘wget’ command to download the raw sequencing file (.fastq.gz). The downloaded files were examined by checking the sequencing reads. Single-end runs with at least 1000 reads and an average read length > 10 bp were considered successful. Meanwhile, pair-end runs with both R1 and R2 having at least 1000 reads and an average read length > 10 bp were considered successful. For runs where either R1 or R2 failed, the workflow kept the passing FASTQ files (R1 or R2) and re-annotated the library as single-ended. Runs that completely failed these criteria were flagged as having download issues and removed from further processing. Additionally, runs with malformed quality scores or sequenced using the ABI Solid color space were similarly flagged and removed.

Following ‘Check reads’, contamination levels were quantified using FastQ Screen [25] and adapter and poly(T/A) sequences were trimmed using Atropos. Reads shorter than 25 bp after trimming were removed. The runs were aligned to the reference genome using HISAT2 with the parameter‘—max_intronlen 300000—known-splicesite-file <splicesite-file>’. Additional alignment statistics were then generated using ‘samtools stats’ [26], ‘bamtools stats’ [27] and ‘CollectRNASeqMetrics’ (https://gatk.broadinstitute.org/hc/en-us). Library complexity was estimated using ‘MarkDuplicates’. Library strandedness was determined by calculating the percentage of reads mapping to each strand using ‘CollectRNASeqMetrics’. A run was assigned as forward stranded if 75% of the reads mapped to the same strand as the gene model, reverse stranded if 75% of the reads mapped to the opposite strand as the gene model or unstranded. Subsequently, MassiveQC generated a set of genomic counts (genic, intergenic, junctions) using ‘featureCounts’ [28]. After ‘Quality control’ and ‘Alignment’, we trained an iForest model (scikit-learn; v0.22.1) using the 21 features generated by the above process and detected outliers based on this model.

iForest model construction and evaluation in D. melanogaster RNA-Seq

The Library Strategy information obtained by SraMongo describes which sequencing method was used (e.g. RNA-Seq, EST, WGS, ChIP-Seq). Currently, the SRA allows for 36 different terms. Here we focused on the most common Library Strategy (RNA-Seq). In MassiveQC, we used an anomaly detection method (iForest) to identify and remove outlier RNA-Seq samples. Initially, the RNA-Seq samples were segregated into training and testing datasets, with the former being used to train the iForest model by utilizing 100 decision trees. We validated the RNA-Seq iForest model using RNA-Seq test data to generate a test dataset consisting of different proportions of EST, WGS or ChIP-Seq runs. To better understand why iForest classified a sample as an outlier, we used Shapley Values to help interpret the model [29]. Briefly, given the RNA-Seq iForest model, local Shapley values were calculated for each sample describing how each feature contributed to moving the sample toward or away from the mean.

Alignment workflow

The alignment workflow generated experiment level (SRX) deliverables, including genomic feature counts and coverage tracks (BigWigs). The workflow began by removing adapter sequences, poly(T/A) sequences and low-quality bases (−q 20) using Atropos [30]. Reads shorter than 25 bp after trimming are excluded. Runs were then aligned to the D. melanogaster Release 6 plus ISO1 MT (GCA_000001215.4) using ‘hisat2—max_intronlen 300000—known-splicesite-file <r6.11 GTF> —rna-strandness <forward or reverse>’. The library strandedness information was obtained from the pre-alignment workflow. Multi-mapping reads and low-quality alignments were removed using ‘samtools -q 20’, and all runs from a library (i.e. SRRs from an SRX) were merged using ‘samtools merge’. We also performed a post hoc analysis using correlation to flag libraries whose runs did not look like technical replicates. The workflow generated genomic coverage counts (genic, intergenic and junction) using ‘featureCounts’ with the options ‘-p -P -C -J’ for pair-end or ‘-J’ for single-end libraries. Stranded genome browser tracks were then generated using ‘bamCoverage—binSize 1—normalizeTo1x 129,000,000—ignoreForNormalization chrX—outFileFormat bedgraph—filterRNAstrand <forward or reverse>’ [31].

Normalization of biological metadata

The normalization of biological metadata was carried out for RNA-Seq samples, by creating a set of keywords for each sample using natural language processing. Briefly, for each sample, we tokenized all author biological metadata, including sample titles and used term frequency inter-document frequency (TFIDF), to identify and weight keywords using NLTK [32]. Next, we aimed to create a set of normalized biological attributes for the four main attribute classes (sex, developmental stage, tissue and cell type). We aggregated author-supplied attributes related to these categories and harmonized class labels by hand (e.g. developmental stage versus dev_stage). For each attribute class, we manually mapped values to the FlyBase controlled vocabulary and imputed missing values by searching the sample keywords.

Analysis of gene expression dynamics

All analyses in this study used FPKM (Fragments per Kilobase of transcript per Million mapped reads) to quantify gene expression. We applied the R package sva to remove batch effects among samples and estimate gene expression variations [33]. Genes with FPKM > 0.1 in at least one sample were retained for further analysis. The evolutionary age of genes was retrieved from GenOrigin [34] and GenTree [35]. The estimates of d_N/d_S values were provided by flyDIVaS [36]. The promoter of each gene was defined as the region 500 bp upstream to 500 bp downstream of the transcription start site. We measured sequence conservation of gene promoters using phastCons and PhyloP, with basewise conservation scores acquired from the UCSC Genome Browser [37]. Transcription factor (TF) binding information was downloaded from GTRD, a database of transcription factor binding sites generated from ChIP-seq data [38]. Bedtools intersect was used to identify the overlap between of TF binding sites and promoter regions [39]. The number of types of different TFs binding in a region was indicative of the TF diversity. The number and diversity of TFs represented the complexity of transcription regulation. GO enrichment analysis was performed using clusterProfiler [40] based on the background of all annotated genes.

MassiveQC detects outliers based on isolated forests using sample self-reporting

When processing a large-scale public sequencing data from SRA, it is essential to capture the data quality and technical experimental details [17]. The SRA stores metadata in addition to RNA-Seq reads. We found that SRA samples were sometimes mis-annotated, and more often, that important information was scattered across fields, databases and subsections of the underlying manuscripts. There was no central source for this information, which made it difficult for us to use a large amount of public data for analysis. Almost all quality control software considered only sequencing quality of the samples and did not address the above issues [18, 19]. Fortunately, many of the critical metadata affect the data and are thus self-reporting. Therefore, we developed a machine learning-based software, MassiveQC, that takes advantage of this self-reporting to (1) eliminate poor quality samples, samples collected for special purposes (e.g. 5′ or 3′ end analysis) and mis-annotated samples (e.g. DNA-seq labeled as RNA-seq) (Supplementary Figure S1) and (2) to separate a major type of sequencing samples of interest (RNA-seq) from a few other types of contaminated sequencing samples (e.g. ChIP-seq) using iForest [24].

To ensure the quality of the sequencing and to gain more information from their self-reporting, MassiveQC used a variety of methods to extract the features of sequencing samples (Figure 1). First, MassiveQC checked and filtered the reads in each FASTQ file after the download (Figure 1A) and integrated FastQ Screen [25] and Atropos [30] to detect contamination and trim adapters, respectively (Figure 1B). Meanwhile, the features of FASTQ were extracted from the results (‘FASTQ features’, Figure 1B and D). Then, the processed FASTQs with high quality were aligned to the reference genome using HISAT2 [41]. Picard and Bamtools were employed to extract alignment statistics, and count statistics were summarized from the featureCounts results (Figure 1C and D). Finally, MassiveQC obtained 21 features from each sample and used these features to build an RNA-Seq iForest model and to identify outliers among all samples (Figure 1E) (see Methods). Given that MassiveQC identifies outliers depending only on read features of samples provided by the users, it can be applied to any other species.

Compared to other tools designed for data quality control, such as FastQC, RNA-SeQC2 and RNA-QC-chain, the advantages of MassiveQC are 1) it can integrate features from multiple dimensions (read, alignment and expression) to make the results more accurate; 2) it generates cutoffs based on the self-reporting of all samples, which is not affected by artificial factors and 3) isolated forests are robust to multi-modal distributions and are applicable to publicly available heterogeneous data (e.g. different tissues and developmental stages) (Table 1).

To evaluate the performance of MassiveQC, we used the large RNA-Seq datasets with 952 runs in the SRA project SRP045429 as an example [42]. The median size of the runs is 346.4 Mb (Figure 2B) and the median running time is 11 minutes (Figure 2A). The above tasks were performed with 8 threads per run on a PBS cluster with 40 nodes whose CPU is Intel(R) Xeon(R) CPU E5-2690 v4. Each node has 256G running memory and 28 threads. Each sample was processed by 8 threads. We then collected 30 samples of other sequencing types as contaminating samples (Supplementary Table S1) to examine the accuracy of MassiveQC. We randomly added k samples (k = 5 to 30) from the 30 contaminating samples to the SRP045429 dataset and used MassiveQC to identify outliers. After 100 simulations, we calculated the mean probability that a contaminating sample was identified as an outlier. In all simulations, contaminating samples were identified as outliers (Figure 2C). In addition, we performed 100 simulations with n normal samples (n = 100 to 800) and 30 contaminated samples to evaluate the stability of the model in the samples of different sizes (Figure 2D). When samples of the target type accounted for more than 90% of the total samples, 99.8% of all contaminating samples would be identified as outliers.

A pipeline for collection and filtration of RNA-seq data from SRA using MassiveQC

Next, we built a pipeline based on MassiveQC for integrating massive RNA-seq data to generate a comprehensive transcriptome atlas of D. melanogaster, by curating the SRA [43], the well-known global repository of raw sequencing data (Figure 3). In this pipeline, we first queried the SRA for ‘D. melanogaster’ [organism]’ and downloaded all associated metadata from the NCBI BioProject, BioSample, Gene Expression Omnibus (GEO) and PubMed repositories using SraMongo which we developed to store and organize SRA metadata (Figure 3A). As of March 2020, there were 62,927 D. melanogaster samples (Supplementary Figure S2). We individually downloaded and processed all 73,539 runs using MassiveQC, except for 107 runs that could not be downloaded and 4923 ABI Solid runs that we excluded to streamline downstream processing (Figure 3B). These results are available on GEO (this study, GSE117217) and the composite tracks are available on FlyBase (‘Oliver lab SRA Aggregated RNA-Seq’) [44].

Then, we explored the marginal contributions of feature importance using Shapley values (Supplementary Figure S1C and D) [29] and found that the most important indicators of library mRNA quality were ‘% rRNA reads’ and ‘% alignment with transcripts’, and that metadata mis-annotation was indicated by ‘% mRNA’ and ‘Mid Gene Body Coverage’. To test the model, we mixed different proportions of samples from other library strategies [such as expressed sequence tags (EST) and whole genome sequencing (WGS)] (Supplementary Figure S1B) with RNA-Seq samples. Our data showed that WGS and ChIP-Seq libraries (DNA-Seq samples) were consistently identified as RNA-Seq outliers by iForest, confirming the ability of iForest to identify mis-annotations (Figure 3E). Using this method, we removed 875 RNA-Seq outliers (Supplementary Figure S1A). We also excluded libraries with low complexity (< 5 K expressed genes with fragments per kilobase of transcript per million mapped reads (FPKM) > 1) or low read depth (< 2 million mapped reads) (Supplementary Figure S3). As a result of these processes, we retained 12,049 high-quality RNA-seq samples for further analysis. After filtering, we used more detailed parameters to carry out the alignment workflow and generate experimental-level (SRX) deliverables including genomic feature counts and coverage tracks (BigWigs). In addition, we normalized metadata from SRA to facilitate downstream analysis. We identified and weighted keywords using NLTK [32] and manually aggregated them into 4 major attribute classes (sex, developmental stage, tissue and cell type) and mapped them to FlyBase-controlled vocabularies (Figure 3D). The metadata are also available at GEO (GSE117217).

Drosophila transcriptome atlas based on 12,049 high-quality public RNA-seq samples

A large collection of samples also facilitates the grouping of results by fly anatomy. We manually annotated the RNA samples by organ, developmental stage and sex using available metadata and explored the expression variation among samples to identify outliers. When parsed by organ anatomy (Figure 4A), we found 28 sources, with the two largest groups being whole body (n = 5225) and head (n = 2436, including the brain). We also parsed by five developmental stages, including egg, embryo, larva, pupa and adult. The adult stage has the largest sample size (n = 7835). Given that the technical biases among our data could confound biological outcomes [45], we applied the R package sva to remove batch effects between samples (Supplementary Figure S4). The large-scale dataset allows us to examine gene expression patterns in depth and systematically. We first examined the breadth of gene expression, defined as the frequency of samples in which a gene was expressed and revealed that 6455 (39.06%) genes were expressed (FPKM >1) in more than 50% samples (only protein-coding genes, lncRNAs and pseudogenes were considered, same below). For example, eukaryotic translation initiation factor 1A (eIF1A), required for mRNA translation [46], was expressed in almost all samples. Meanwhile, we found that 4559 (27.59%) genes were expressed in less than 10% samples (Supplementary Figure S5), suggesting their distinct gene expression spectrums. Another example was the olfactory receptors (ORs). The olfactory system of Drosophila contains 61 OR genes [47] and the saturation curve showed that the more head tissues, the higher the detection rates of OR genes, with Or43b detected in 93.44% of all head samples (2186) (Supplementary Figure S6).

It is not easy to visualize expression distances between 12 K samples. Therefore, we next subjected the samples to low-dimensional embedding and identified clusters [48]. The 2D Uniform Manifold Approximation and Projection (UMAP) plot showed that gene expressions of functionally relevant samples (organs, tissues or cells) were clustered together, for example, ovary and egg, as well as midgut and digestive system (Figure 4B), suggesting the common gene expression programs. Likewise, tissues were also grouped by developmental stage. Many adult samples were sexed and they clearly separated by sex in the UMAPs (Figure 4C and D), which was consistent with previous studies that sex dimorphism is fully manifested at the adult stage [3, 16, 49–51].

Dynamic pattern of gene expression

We further applied the two-dimensional uniform manifold approximation and projection (UMAP) implemented in Seurat to calculate the standardized variance of each expressed gene [48]. Among the 16,679 genes, roX2, a long non-coding RNA that is functionally redundant with roX1 and essential for the assembly of the dosage compensation complex [52, 53] exhibited the highest expression variance (Figure 5A). Consistently, the expression of roX2 was dynamic during development and was limited to a few tissues based on FlyBase expression data [50]. Furthermore, our observation was also consistent with a recent study that roX2 was the top male-biased specific gene at the single cell level [16], confirming our evaluation of gene expression dynamics. We also introduced another metric tau [54] to evaluate the expression specificity of genes and observed the strong consistency between gene expression dynamics and expression specificity (Spearman correlation, Rho = 0.77, P = 0) (Supplementary Figure S7). We then ranked all genes based on their expression variance and showed the expression pattern of the top dynamic genes (5%) and the bottom dynamic genes (also known as stably expressed genes) (5%). As expected, the most dynamically expressed genes were expressed in one or limited tissues (Figure 5B). For example, a set of genes was preferentially expressed in larval anal pad (box1) and another set of genes exhibited a higher expression level in the adult retina than in other tissues (box2). In contrast, the least dynamically expressed genes were found in a wide range of organs or tissues (Figure 5C).

For an intuitive understanding of their distinct expression dynamics, we displayed the expression profile of Lsm11, which was unevenly distributed in the ovary samples and expressed at the adult stage (Figure 5D). We also confirmed and extended these observations using single cell data from Fly Cell Atlas [16]. Once again, and congruently, Lsm11 was found to be expressed in cells of fat body and ovary (Figure 5E). Interestingly, when we further examined the expression pattern of Lsm11 in different cell types of ovary, we observed that it was preferentially expressed in young germ cell, 16-cell germline cyst in germarium region and oviduct (Supplementary Figure S8A), showing remarkable expression dynamics at single-cell resolution.

However, eukaryotic translation elongation factor 2 (eEF2) was widely expressed in nearly all samples (12,041 out of 12,049) (Figure 5F). Similarly, its expression profile at the single-cell level recapitulated the above observation (Figure 5G). We also showed the expression pattern of eEF2 in different cell types of ovary and found that it was expressed in the majority of cells (Supplementary Figure S8B). Considering that the ubiquitously expressed genes are associated with frequent protein-DNA-binding events [55], we further examined whether these genes were involved in complex expression regulation.

Regulation of genes with different expression dynamics

Gene expression in organism is modulated by spatial–temporally dynamic regulatory programs [56]. To further reveal whether the genes with dynamic expression have distinct evolutionary and regulatory patterns, we first investigated the association of gene types with expression dynamics at the basic level of encoded gene product. We observed that protein-coding genes exhibited the lowest expression variation, followed by lncRNAs and pseudogenes (based on FlyBase annotation) (Wilcoxon test) (Figure 6A), consistent with previous reports that lncRNAs and pseudogenes are expressed in a highly tissue- or condition-specific manner [57]. We also looked for correlations with the time when a gene first arose in the lineage. Interestingly, the origination time of genes (gene age) was negatively correlated with its expression dynamics (Spearman correlation, Rho = −0.43, P = 0), with younger genes displaying higher expression variation (Figure 6B, Supplementary Figure S9). Previous studies have shown that young genes tend to have a testis-specific expression and subsequently acquire broader expression (Supplementary Figure S10) [34, 58]. Our work showed that this gradual repurposing of genes is a general trend.

We next surveyed the forces governing the evolution of Drosophila genes at the protein level. First, we examined the evidence of selection on dynamically expressed genes using d_N/d_S value, which was defined as the ratio of the rates of nonsynonymous and synonymous codon substitutions, with high values (> 1) indicating positive selection and low values indicating purifying selection [59]. The expression dynamics were positively correlated with selection strength and relatively stably expressed (or widely expressed) genes were enriched for signatures of purifying selection (Supplementary Figure S11), suggesting that purifying selection might constrain the dynamics of gene expression.

While adult phenotypes are entirely dependent on the successful completion of increasingly complex gene expression patterns during development, it is thought that pathways evolve in reverse. In other words that as long as the final phenotype is achieved, it is possible to exchange mechanisms for initiating development [60]. If this is applicable to gene expression, then we might expect to see more variance in early developmental stages. Surprisingly, we discovered that late expressed genes showed higher expression variation than early-expressed genes (Figure 6C). This suggested that even if early development was evolutionarily fungible, there was tight regulation at micro-timescales. The overall trend was linear, but showed two peaks at adult 14-day and 24-day ages. This result might be explained by the mutation accumulation theory of aging, where the force of purifying selection decreases with age, as proposed by Sir Peter Medawar [61]. Moreover, the genes expressed in testis showed significantly larger expression variations than those in other organs (Figure 6D), the pattern might result from rapid gene evolution in testis and pervasively accessible chromatin allowing extensive transcription [56, 62]. We also observed that genes expressed in the ovary were much less dynamic than those in other tissues, which might because genes important for Drosophila ovarian function are highly conserved during evolution [63], suggesting that these genes might also be expressed in other tissues and be pleiotropic. Collectively, these observations (Figure 6B–D) provided a link between gene expression patterns in an extant species and the molecular evolution and life history theory of the lineage (Medawar’s hypothesis).

To evaluate whether our observations of the associations between expression dynamics and evolutionary constraints could be extended to regulatory programs, we examined the sequence conservation in promoter regions which governed the transcription initiation and gene activity [64]. Expectedly, the promoters of genes with higher expression dynamics were significantly, albeit slightly, less evolutionarily conserved (Supplementary Figure S12). Interestingly, we also found that the promoters of genes with lower expression dynamics harbored more TF-binding sites and higher diversity of TFs (Figure 6E, Supplementary Figure S13) (Spearman correlation, Rho = −0.53, P = 0). These results were reminiscent of high-occupancy target (HOT) regions, defined as the segments that were bound by high amount of TFs, and genes related to HOT regions were stably expressed [55, 65], with their expression variability resembling that of housekeeping genes [66]. Additionally, the transcripts derived from less dynamically expressed genes were bound by more RNA-binding proteins (RBPs) than those from highly dynamically expressed genes (Figure 6F, Supplementary Figure S14). These data indicated a stronger and more complex transcriptional and post-transcriptional regulation of these less dynamic genes. Given the diverse expression patterns, genes that need to be expressed widely require more TF-binding sites to ensure that they could be expressed in essentially all cell types. If the levels of proteins encoded by these genes were also in a narrow tolerance range, then they might be more subject to translational control by RBPs.

Function and phenotypic essentiality of genes with different expression dynamics

Evolutionary constraints and regulatory complexity are strong predictors of functional importance. To extend this inference, we hypothesized that the dynamically expressed genes might perform lineage-specific functions while the less dynamically expressed genes might be associated with the maintenance of basic cellular functions (acting as ‘housekeeping genes’). Considering the different functional paradigms between protein-coding genes and lncRNAs [57], we first focus on protein-coding genes. In support of our notion, the dynamically expressed genes were significantly enriched in developmental process, immune response and neuropeptide signaling pathway gene ontology terms (Figure 7A), while the less dynamically expressed genes were enriched in protein phosphorylation, cell junction organization and GTPase-mediated signal transduction terms (Figure 7B), indicating the different functional categories of these two types of genes.

To further assess whether the less dynamically expressed genes were required for housekeeping functions (i.e. are essential), we assigned genes to different classes of phenotypic severity based on a public large-scale transgenic short-guide RNA library that enabled efficient disruption of specific target genes [67]. We compared the degree of expression variation for genes in these different classes and revealed a clear association between expression dynamics and the severity of phenotypes when those genes were disrupted (Figure 7C). For example, genes required for lethality showed the lowest expression dynamics. Using another dataset of essential genes predicted by a machine learning algorithm [68], we found that essential genes were enriched for less dynamically expressed genes (Supplementary Figure S15). Collectively, our data reinforced and extended the idea that genes with less dynamic expression have important, conserved functions in Drosophila.

Considering our relatively poor understanding of specific lncRNA functions, it is difficult to compare expression dynamics with functional importance. To get around the poor functional annotation, we borrowed functional information from protein-coding genes based on a previous method [69], by constructing a co-expression network between protein-coding genes and lncRNAs. We used this functional relationship to infer lncRNA functions. With P < 0.01 and absolute Pearson correlation coefficient R > 0.65, a total of 21,613 co-expression pairs were identified between 2383 coding genes and 89 lncRNAs (Figure 7D). Protein-coding genes significantly correlated with dynamically expressed lncRNAs were involved in specific differentiation, such as the sperm-associated functions of spermatid development and differentiation (Figure 7E), while protein-coding genes correlated with stably expressed lncRNAs were involved in basic biological processes, such as cell–cell signaling and cell junction assembly (Figure 7F). These observations on protein-coding genes and lncRNAs indicated that expression dynamics and functional categories were intertwined in these gene pairs.

Conservation and divergence of gene content and expression profile between Drosophila and human

As a premier model species, Drosophila is routinely used to elucidate the progress and underlying molecular mechanisms of human development and diseases [4, 70–73], under two assumptions: 1) a large fraction of human genes, especially disease-associated genes, have counterparts in the Drosophila genome; 2) the genetic regulatory programs are largely conserved between the two. However, how these two assumptions interact has not been systematically characterized. For example, approximately ranging from 65% to 77% of human disease genes have orthologs in Drosophila [71, 74], but how these genes are regulated in the two species is underexplored. Our dataset provided important genome-wide information on the transcriptional regulation of orthologs.

To this end, we classified human genes into different categories based on their phenotypic severity by integrating a dataset of essential genes [75] with a dataset of disease-associated genes [76]. This analysis obtained 2460 vital genes, 1145 vital and disease-suppressing genes, 3150 disease-suppressing genes and 10,463 other genes, with gradually decreasing severity of loss-of-function phenotypes. We matched these with 20,400 ortholog pairs between human (9223 genes) and Drosophila (6937 genes) using Ensembl Compara [77]. Interestingly, we observed a clear and significant association between orthology in Drosophila and phenotypic severity in human loss-of-function (Figure 8A). Specifically, 68.5% of human vital genes were conserved in the Drosophila genome, followed by vital and disease-suppressing genes (65.33%), disease-suppressing genes (49.3%) and other genes (42.34%), suggesting that functionally important genes have been conserved by evolutionary forces. The utility of Drosophila is partly due to the lower number of paralogs, but we also discovered that 40.57% of human vital genes had a one-to-one (1:1) orthology in Drosophila, significantly higher than genes in other three phenotypic groups. Overall, these analyses illustrated that the vital or disease-suppressing genes in human ‘obey’ strict 1:1 orthology relationships, suggesting that most of them could be investigated in Drosophila avatars.

Much of the knowledge of the genetic programs underlying human biology was first discovered in model organisms like rhesus macaque, mouse and Drosophila. Although the expression conservation between human and mammalian models has been elucidated [78], such analyses on human genes and Drosophila orthologs are understudied. Therefore, we evaluated the expression conservation and dynamics of 1:1 orthologs in 9 matched tissues (Figure 8B–D, Supplementary Figure S16) and found that the expression correlations ranged from 0.08 to 0.67 (Spearman correlation coefficient), with the lowest value in muscle (Figure 8B) and the highest value in ovary (Figure 8C), suggesting the high conservation of gene expression in ovary over hundreds of millions of years of divergence. Interestingly, the general linear relationship did not depend on gene category, but when looking closely, we found that the vital genes exhibited the highest expression conservation than other types of genes in 6 of the tissues, with the highest value of 0.65 in ovary, while the other genes exhibited the lowest conservation in 7 tissues, with the lowest value of 0.08 in muscle. We also found that the expression similarity of genes expressed in testis was considerably lower than other tissues, especially for disease-suppressing and other human genes (Figure 8D). These observations strongly suggested that the clinical phenotypes of human vital and vital and disease-suppressing genes were more likely to be recapitulated by the Drosophila model. For example, the high expression conservation and low expression variability of human Drosophila ortholog pairs in ovaries made them an ideal organ for translating disease etiology research between human and Drosophila.

Meng-Wei Shi is currently a PhD candidate at Hubei Key Laboratory of Agricultural Bioinformatics, College of Life Science and Technology and College of Biomedicine and Health, Huazhong Agricultural University, China.

Dan-Yang Wang is currently a PhD candidate at Hubei Key Laboratory of Agricultural Bioinformatics, College of Life Science and Technology and College of Biomedicine and Health, Huazhong Agricultural University, China.

Justin M. Fear is a postdoc researcher at Section of Developmental Genomics, National Institute of Diabetes and Kidney and Digestive Diseases, National Institutes of Health, Bethesda, USA.

Lu Chen is a postdoc researcher at Hubei Key Laboratory of Agricultural Bioinformatics, College of Life Science and Technology and College of Biomedicine and Health, Huazhong Agricultural University, China.

Yi-Xuan Tu is currently a PhD candidate at Hubei Key Laboratory of Agricultural Bioinformatics, College of Life Science and Technology and College of Biomedicine and Health, Huazhong Agricultural University, China.

Hong-Shan Liu is a master’s graduate at Hubei Key Laboratory of Agricultural Bioinformatics, College of Life Science and Technology and College of Biomedicine and Health, Huazhong Agricultural University, China.

Yuan Zhang is a master’s graduate at Hubei Key Laboratory of Agricultural Bioinformatics, College of Life Science and Technology and College of Biomedicine and Health, Huazhong Agricultural University, China.

Shuai-Jie Zhang is a master’s graduate at Hubei Key Laboratory of Agricultural Bioinformatics, College of Life Science and Technology and College of Biomedicine and Health, Huazhong Agricultural University, China.

Shan-Shan Yu is a master’s graduate at Hubei Key Laboratory of Agricultural Bioinformatics, College of Life Science and Technology and College of Biomedicine and Health, Huazhong Agricultural University, China.

Brian Oliver is a professor of Section of Developmental Genomics, National Institute of Diabetes and Kidney and Digestive Diseases, National Institutes of Health, Bethesda, USA.

Zhen-Xia Chen is a professor of Hubei Hongshan Laboratory, Hubei Key Laboratory of Agricultural Bioinformatics, College of Life Science and Technology, College of Biomedicine and Health, Interdisciplinary Sciences Institute, Huazhong Agricultural University, Wuhan 430070, China; Shenzhen Institute of Nutrition and Health, Huazhong Agricultural University, Shenzhen 518000, China; Shenzhen Branch, Guangdong Laboratory for Lingnan Modern Agriculture, Genome Analysis Laboratory of the Ministry of Agriculture, Agricultural Genomics Institute at Shenzhen, Chinese Academy of Agricultural Sciences, Shenzhen 518,000, China. Her research group has been working in evolutionary and developmental genomics.