REVA as A Well-curated Database for Human Expressionmodulating Variants

Yu Wang ,Fang-Yuan Shi ,Yu Liang ,Ge Gao

1State Key Laboratory of Protein and Plant Gene Research,School of Life Sciences,Biomedical Pioneering Innovation Center(BIOPIC)&Beijing Advanced Innovation Center for Genomics(ICG),Center for Bioinformatics(CBI),Peking University,Beijing 100871,China

2Human Aging Research Institute,School of Life Sciences,Nanchang University,Nanchang 330031,China

Abstract More than 90% of disease-and trait-associated human variants are noncoding.By systematically screening multiple large-scale studies,we compiled REVA,a manually curated database for over 11.8 million experimentally tested noncoding variants with expression-modulating potentials.We provided 2424 functional annotations that could be used to pinpoint the plausible regulatory mechanism of these variants.We further benchmarked multiple state-of-the-art computational tools and found that their limited sensitivity remains a serious challenge for effective large-scale analysis.REVA provides high-quality experimentally tested expression-modulating variants with extensive functional annotations,which will be useful for users in the noncoding variant community.REVA is freely available at http://reva.gao-lab.org.

KEYWORDS Noncoding variant;Expression-modulating variant;Massively parallel reporter assay;Database;Benchmark

Introduction

Noncoding regions occupy the majority of the human genome[1].It has been demonstrated that noncoding variants can affect the regulation of genes [2],and more than 90%of disease-and trait-associated variants are noncoding variants [3].Noncoding variants that could affect gene expression can be considered as expression-modulating variants[4].Several experimental assays have been developed to characterize expression-modulating variants.Genome editing technologies such as transcription activator-like effector nucleases(TALENs),zinc finger nucleases(ZFNs),and clustered regularly interspaced short palindromic repeats with Cas9 nuclease (CRISPR/Cas9) provide high-quality validated data but are generally low throughput [5-7].Recently developed massively parallel reporter assays(MPRAs)can identify transcriptional regulatory elements in an efficient way,allowing systematic screening of tens of thousands of genetic variants for pinpointing the causal variants of complex traits [4,8,9].All expression-modulating variants stored in MaveDB[10]are validated by the MPRA experiments.MPRAs have generated over 10 million human expression-modulating variants [11];however,only around 30 thousand of them have been collected by MaveDB without any functional annotation,which hinders the further utilization of these data.

Although experimental assays for characterizing noncoding expression-modulating variants have generated a huge amount of data,it is still inadequate for covering all noncoding variants identified in human genomes.Therefore,multiple computational tools have been developed for identifying expression-modulating variants(Table 1).Transcription factors(TFs)could regulate genes through binding to sequence motifs [12],and noncoding variants could affect gene regulation by changing motifs [13].FunSeq2 integrated a module for detecting motif-breaking and -gain events through the change of position weight matrix (PWM) and other functional annotations to prioritize cancer driver mutations [14].Methods based on machine learning have been used wildly in biological researches [15].CADD [16] used support vector machine (SVM) to classify variants into functional and nonfunctional variants,and GWAVA [17] used random forest to predict disease-related variants.Both CADD and GWAVA were based on supervised learning methods,while Eigen [18] implemented unsupervised learning methods to classifyvariants.All these tools highly depend on existing annotations at corresponding loci.In 2015,Alipanahi et al.[19]developed DeepBind based on convolutional neural networks (CNNs) to predict the binding affinity between TFs and DNA or RNA binding proteins and RNA.DeepSEA[20]applied similar methods to predict the effect of noncoding variants on binding affinity and then classified variants through logistic regression into functional or nonfunctional groups.All tools mentioned above identified expression-modulating variants through indirect inference,because they were not trained on expression-modulating variants or expression-related data.EnsembleExpr [21]used MPRA data to train an ensemble-based model for characterizing expression-modulating variants directly.ExPecto [22]ab initiopredicted the variants’ effects on gene expression from 40-kb promoter-proximal sequences and then pinpointed expression-modulating variants.However,there is no comprehensive evaluation of these computational tools based on high-quality expressionmodulating variants;therefore,it is difficult for users to choose appropriate tools for their tasks.

Here,we present a repository for expression-modulating variants (REVA).The current release of REVA consists of over 11.8 million experimentally validated expressionmodulating variants in the human genome,curated with extensive functional annotations.We further benchmark seven popular computational tools in identifying expression-modulating variants[14,16-18,20-22]based on high-quality data in REVA.All data and benchmarking results are publicly available at http://reva.gao-lab.org.

Construction and content

Data collection and integration

To ensure unified and high-quality data,all records in REVA were collected and curated using a standard procedure (Figure 1).We used a list of keywords,“MRPA”,“STARR-seq”,“CRE-seq” with “mutation”,“variant”,and“variation”,to retrieve publications from PubMed (https://pubmed.ncbi.nlm.nih.gov/) and then manually checked the abstracts and full texts of the matching publications to obtain literatures that experimentally validated the effects of noncoding expression-modulating variants.

Figure 1 Overview of the structure of REVA

For filtered literatures,we extracted related information of the variants from the main texts as well as supplementary materials of publications and converted them to the same format(Table 2).Variants that failed to be mapped to bothGRCh37 and GRCh38 were removed.Variants only mapped to the coding region were also removed.For missing information,we used “.” as a placeholder.In addition,the detailed protocols and raw data of the experiments were also extracted.

Table 1 Properties of involved computational tools

Table 2 Variant information extracted during the data collection process

For variants with the same chromosome,genome location,reference allele,alternative allele,and experimental cell line from different publications were subjected to a meta-analysis to integrate data.The harmonic meanPvalue (HMP) method [25] was used in the meta-analysis,and the cutoff for the metaPvalue was set to 0.001 to generate the meta-label.The variants involved in the metaanalysis but without a rawPvalue were also removed.

The label of variants was given based on the cutoff for the adjustedPvalue or metaPvalue,and then variants were classified into positive variants and negative variants based on label or meta-label.If the variant’s label was 1,the variant was a positive variant and considered to have effects on gene expression;otherwise,it was a negative variant without effect on gene expression.

Database construction

All manually curated variant data,as well as metainformation,were stored in MongoDB (https://www.mongodb.com/)at two levels:accession and variant data(Figure S1).Each accession entry consisted of an accession number,created time,last updated time for the accession,information about the assay used in the publication (method type,original reference genome version,link to raw data,and summary of the assay),and the reference.Variant data included all related information about the variant,and each variant data entry was linked to one accession.For the data involved in the meta-analysis,the variant data contained the results of the meta-analysis and were linked to all related variants and accessions.

We also integrated DisGeNET v7.0 (https://www.disgenet.org/) [26] variant-disease associations,GWAS catalog (https://www.ebi.ac.uk/gwas/),ClinVar (https://www.ncbi.nlm.nih.gov/clinvar/),COSMIC(https://cancer.sanger.ac.uk/cosmic),and three-dimensional interacting genes and chromatin state from 3DSNP [27] to our database for providing more variant information.

Variant annotation

In efforts to pinpoint plausible regulatory mechanisms for these variants,we used 2403 trained CNNs to annotate the functional effects of sequence variations[28]based on 1249 TF binding profiles,766 histone modification profiles,280 DNA accessibility profiles,and 108 DNA methylation profiles from the recent Encyclopedia of DNA Elements(ENCODE) data.

NVIDIA Tesla P100 Graphics Processing Units with the implementation on the deep learning framework TensorFlow (https://www.tensorflow.org/) and Python(https://www.python.org/) were used for training models.We adopted stochastic gradient descent (SGD) as the optimizer,and the initial learning rate was 0.01.

The final output layer of the CNN model was a fully connected layer with a sigmoid function used to scale the output between 0 to 1.The input layer was a one-dimensional convolution layer with the thresholded rectified linear unit(ReLU)as the activation function.Next,the maxpooling layer was performed to reduce the complexity of the data.Then,the dropout layer was used to mitigate the overfitting problem.The next two layers were a fully connected layer with thresholded ReLU as the activation function and a dropout layer.

For TF binding,histone modification,and DNA accessibility models,the positive data for training CNNs were the 200-bp sequences centered on the peak in ENCODE profiles.Then we removed positive sequences from the human reference genome and split the rest into 200-bp bins.Random sampled 200-bp bins with the same number of positive data were used as negative data.For DNA methylation,the 200-bp sequences centered on the target base with the methylation rate more than 0.5 or less than 0.5 in wholegenome bisulfite sequencing(WGBS)data were considered as positive data and negative data,respectively.One-hot encoding was conducted to transform each sequence to a 200 ×4 binary matrix for model training.

A five-fold cross-validation strategy was used to train models.During each iteration of model training,15%of the input data were randomly selected as the independent testing dataset to evaluate model performance.The remaining data were split with 70%to train models and 15%as the validation dataset to optimize parameters.Model performance was evaluated with the area under the receiver operating characteristic curve(AUROC)and area under the precision-recall curve (AUPRC) to test the sensitivity and specificity,and models with the best performance were selected for variant annotation.An average AUROC and an average AUPRC of 2403 models were reported.

To character the binding affinity changes of the variant,we used 2403 trained CNNs to predict on 200-bp sequences centered with the reference allele and alternative allele,respectively.For each chromatin profile,the log2fold change (as the method shown in DeepSEA) [20] was calculated as the variant effect on chromatin profile.Specifically,

where Prwas the prediction of sequence with reference allele,and Pawas the prediction of sequence with alternative allele.

Furthermore,we incorporated 13 DNA physicochemical properties and 8 evolutionary features into the annotation pipeline.The 13 physicochemical properties were calculated as described by Li et al.[29],and 8 conservation scores were downloaded from UCSC Genome Browser(http://genome.ucsc.edu/).

Benchmarking

To prepare the benchmarking dataset for evaluating the performance of state-of-the-art computational tools in calling expression-modulating variants based on the curated data in REVA,we first excluded loci tested in mice (n=15,152).There are overlapping variants between the training datasets of state-of-the-art tools and the REVA benchmarking dataset.If the benchmarking dataset contains these variants,the performance of related tools will be overestimated.To avoid the influence of these variants and make a fair comparison,we further removed variants (n=47,518)that were either found in the GWAVA [17] and EnsembleExpr[21]training datasets or used to compute the empirical background distributions by DeepSEA [20].For the remaining 5,809,991 loci (37,816 positive and 5,772,175 negative),we ran CADD (v1.4,https://cadd.gs.washington.edu/download) [16],DeepSEA(http://deepsea.princeton.edu/),EnsembleExpr(https://github.com/giffordlab/EnsembleExpr/),and ExPecto (https://hb.flatironinstitute.org/expecto/)[22],and used precomputed score sets of Eigen (v1.1,http://www.columbia.edu/～ii2135/download.html) [18],FunSeq2 (v2.1.6,http://funseq2.gersteinlab.org/downloads) [14],and GWAVA (https://www.sanger.ac.uk/science/tools/gwava),to obtain the corresponding predicted score for evaluation.The thresholds used in the evaluation were those recommended by the corresponding papers or official websites (Table S1).

All variants in the benchmarking dataset were variants with expression-modulating potential.One of the biological mechanisms by which disease-related or phenotype-related variants function is having effects on gene expression regulation [30].Pinpointing disease-related or phenotyperelated variants is more useful for biomedical researches.Therefore,we further selected the GWAS,ClinVar,and HGMD subsets of the benchmarking dataset to test these tools’ power.

Results

Characterization and distribution of expression-modulating variants

All curated expression-modulating variants were validated by experiments,and we applied standard data collection and integration procedure to ensure the high-quality data with the unified format.By the end of November 2019,REVA consisted of 11,862,367 entries covering 5,948,789 experimentally tested noncoding loci across 18 cell cultures from 14 publications[4,8,11,31-41].We first excluded loci tested in mice (n=15,152) and with more than one alternative allele (n=26,276).Among the remaining 5,907,361 loci (34,700 positive and 5,872,661 negative),most were located in intergenic(positive:49.96%,negative:53.83%)and intronic(positive:35.96%,negative:39.62%)regions(Figure 2A;Table S2).We found that both positive and negative variants were unevenly distributed on chromosomes,and no variants were located on the Y chromosome (Figure 2B;Table S3).Specifically,fewer positive variants were located on chromosomes 1,3,5,13-15,and 21,and the X chromosome,and more positive variants were located on chromosomes 6,8,10-12,and 16-20.Fewer negative variants were located on chromosomes 9,13-15,21,and 22,and the X chromosome,and more negative variants were located on chromosomes 1-8,10-12,and 16-20.Biochemical activities were detected for 93.53%positive and 90.80% negative cases in at least one cell culture (Figure S2A;Tables S4 and S5).Of note,more positive than negative variants were found in TF binding regions,highlighting the contribution of TF binding changes to expression modulation(Figure S2B).

Figure 2 Annotation of the variants in REVA

Figure 3 Performance of involved tools on the benchmarking dataset

Extensive functional annotation of expression-modulating variants

We used 2403 trained CNNs to annotate the functional effects of expression-modulating variants[28].Most of the trained CNNs were accurate,with an average AUROC of 0.908 and an average AUPRC of 0.904.Among the 5,789,688 variants annotated,both positive and negative variants were found to lead to significant changes in binding affinity for 22 and 12 TFs on average,respectively,which also suggested that expression-modulating variants may affect gene expression regulation through changing the binding affinity of TFs.Moreover,8.72% positive and 3.56% negative variants were located at evolutionary conserved loci (phastCons100way score ＞0.6).

Benchmarking of state-of-the-art computational tools

To evaluate the power of state-of-the-art computational tools in calling expression-modulating variants,we further benchmarked multiple state-of-the-art computational tools based on the curated data in REVA.With the benchmarking dataset containing 5,809,991 loci (37,816 positive and 5,772,175 negative),we found that 1289 could not be predicted by DeepSEA(since their evolutionary features were not available),and 560,577 were not included in the

precomputed score set of Eigen,FunSeq2,and GWAVA,so we further excluded these 561,866 cases from follow-up analysis.Meanwhile,as EnsembleExpr could not finish the whole benchmarking dataset in a reasonable amount of time,we assessed its performance based on the average metrics over 5 randomly sampled sub-datasets with 368 positive and 56,026 negative cases on average.

Overall,the best-performing tool was DeepSEA,with the highest AUROC and F1 score (Figure 3A and B;Table S6).All tools performed well in terms of specificity but poorly in terms of sensitivity.EnsembleExpr had the highest sensitivity but the lowest specificity,whereas ExPecto showed the best specificity and worst sensitivity(Table S6).

There were 52,672 variants in the benchmarking dataset that overlapped with the GWAS catalog(v1.0.2),and 658 of them were positive variants.All tools performed better on variants overlapping with the GWAS catalog,and DeepSEA still had the best performance (Figure 3C).Meanwhile,ExPecto and FunSeq2 showed better performance on variants at evolutionarily conserved loci,while DeepSEA displayed moderate performance(Figure 3D).

The coverage and quality of training data may contribute significantly to the performance of machine learning-based models[44].To test whether variants from different cell lines would affect the performance of these tools,we further evaluated these tools separately on seven cell lines(Table S7).On GM18507,GWAVA-Unmatch performed best;on HEK293T and NA12878&NA19239,Eigen-PC had the highest F1 score;DeepSEA had the best performance on HepG2,K562,and K562_GATA1;and CADD performed best on SH-SY5H (Figure 3E),which suggested that the diversity of the original training data contributes to the performance differences of these tools.Of note,thus far,only ExPecto outputted cell type-specific scores for various tissues.

To provide a further explanation of the potential mechanisms of disease-related variants,we evaluated the benchmarking dataset on disease-related variants.There were 1400 variants in the benchmarking dataset that overlapped with HGMD (2019.3 professional),and 69 of them were positive variants(Table S8).Moreover,8 of 69 variants were verified to regulate gene expression by independent experiments;40 of 69 variants were associated with diseases such as colorectal cancer,nervous system diseases,and autoimmune diseases.To test computational tools’power on disease-related variants,we compared their performance on these variants.All tools performed better on variants overlapping with HGMD,and DeepSEA still had the best overall performance(Figure 3F),same on variants with class “DM?” and “DFP”.Eigen-PC showed the best performance on variants with class “DP”.Interestingly,ExPecto performed best on variants with class “FP” but worst on variants with other classes.We also evaluated variants overlapping with ClinVar (2019.10.08),and DeepSEA had the best overall performance,and Eigen showed the best performance on “Drug response” related variants(Figure S3;Table S9).

Web interface

REVA(http://reva.gao-lab.org)provides an interactive web interface for users to explore all data entries and analysis results(Figure 4,Figure S4).Users can start a quick search by chromosome position,rs ID,gene name,ensemble gene ID,or disease name.“Advanced search” provides a customized search and batch search for users.The query result is presented as a table,which includes basic information,expression information (such as the label,effect size,and adjustedPvalue),and the related genomic region.Users can directly click the link of position and rs ID to access UCSC Genome Browser and dbSNP (https://www.ncbi.nlm.nih.gov/snp/) for more information.Users can also click the “details” link for more information.The detail page contains eight modules:“Basic Information”,“Cell Line and Expression”,“Three-dimensional Interacting Gene”,“Chromatin State”,“Disease and Phenotype”,“Meta Sources” (only available for variants involved in meta-analysis),“Accession”,and “Annotation”.In the“Annotation” module,chromatin profile features are rendered as a heatmap by cell line and a boxplot by category,and DNA physicochemical properties and evolutionary features are presented as responsive tables.Users can download the annotation for further analysis.Moreover,we also provide benchmarking results of state-of-the-art computational tools.Users can download all variants in REVA and the benchmarking dataset through the “Download” page.

Explore plausible regulatory mechanisms of expressionmodulating variants

Autoimmune diseases are caused by the abnormal immune response to attack and damage functional tissues due to complex interactions between environmental and genetic factors [45].GWAS and fine-mapping studies have identified thousands of noncoding variants associated with autoimmune diseases [46].Since the mechanisms of autoimmune disease are complicated,pinpointing causal variants and exploring their possible functional mechanisms remain a challenge [47].

Ankylosing spondylitis is a kind of chronic autoimmune disease,but the pathogenesis remains unclear[48].On the advanced search page of REVA(Figure S5),we filtered the label to positive and searched with “ankylosing spondylitis”.The search result contained 8 entries,and among them,the variant rs4456788 (near theICOSLGlocus) had the largest effect size tested in HepG2 cell line and was considered to repress expression.It was also tested in K562 cell line and resulted in the same conclusion.Through the annotation module of the detail page,we found that in HepG2 cell line,the alternative allele of rs4456788 could decrease the binding affinity of TFs MAZ and FUS.MAZ has been proven to have bidirectional transcriptional regulation [49],and FUS has a transcriptional activation function[50].It could be the possible regulatory mechanism of rs445678,and this might be helpful for further researches on the mechanism of pathogenesis of ankylosing spondylitis.

Figure 4 Illustration of the Web interface of REVA

Discussion

REVA is a database specifically designed for storing experimentally validated expression-modulating data.It currently consists of 11,862,367 entries covering 5,948,789 experimentally tested noncoding loci across 18 cell cultures.Both experimentally validated expression-modulating variants and meta-information about assays were curated.Comparing with the existing database,REVA is the largest database designed for curating experimentally validated expression-modulating noncoding variants specially.Besides,we provide 2424 functional annotations,including TF binding,epigenetic modifications,DNA accessibility,DNA physicochemical properties,and evolutionary features.

Most of the variants in REVA were located in intergenic and intronic regions and were unevenly distributed on chromosomes.Several factors may contribute to the uneven distribution.First,it has been well demonstrated that the functional elements are unevenly distributed across chromosomes [51,52].Consistently,we found that the numbers of both positive and negative variants were highly correlated with the gene numbers across all chromosomes(Pearson’sr=0.80,P=2.6 × 10-6for positive variants;Pearson’sr=0.82,P=7.1 × 10-7for negative variants).Moreover,technical challenges counted too.In particular,the Y chromosome had long been taken as a “genetic wasteland” [53] and excluded from genomic analyses for quite some time due to its genetic and structure complexities[54].Although this idea has been shifted with more researches on chromosome Y,the underrepresentation of chromosome Yon commonly used arrays still exists[55].We also noticed that certain experimental designs may lead to reporting bias [8,33,38-41].However,after removing data generated from studies designed for assessing particular regions [38]or elements[8,33,39-41],we found that the uneven distribution remains.

Furthermore,we provide a high-quality benchmarking dataset for evaluating state-of-the-art computational tools designed for identifying expression-modulating variants as well as benchmarking results of multiple published computational tools as a reference for users to select the best tools for their particular tasks.Overall,all seven tools have high specificity but low sensitivity.DeepSEA has the best performance on the whole benchmarking dataset in terms of AUROC and F1 score,and all tools have better performance on disease-related or phenotype-related variants,suggesting that the diversity of the original training data of these tools contributes to different performance across different benchmark subsets.We noticed that not all tools involved in the benchmark were designed for identifying expression-modulating variants originally,and a “negative” expression-modulating noncoding variant might also be associated with disease via nontranscription mechanisms like epigenetic marks [56] or chromatin structuration [57].

It should be noted that not all variants collected in our database were tested by identical experimental protocols.Non-saturation mutagenesis-based studies examine several elements at a time,and each fragment usually contains one variant,with the effect size calculated by counting reads directly [8] or employing a linear model [32].Meanwhile,saturation mutagenesis-based studies focus on a few elements;each fragment contains two or more variants,and the effect size is calculated through linear regression [41].Protocol details for each variant were documented during curation to help users interpret records effectively(Figure S6).

We believe that this database will be useful for not only computational but also bench biologists in genomics,bioinformatics,and genetics communities,and we will keep the resources updated with new data and annotations that emerge in the coming years.

Data availability

REVA is freely accessible at http://reva.gao-lab.org.

Code availability

Source code for all analyses and benchmarking is available on GitHub at https://github.com/gao-lab/REVA-Data_-Source_Code.

CRediT author statement

Yu Wang:Methodology,Software,Data curation,Formal analysis,Visualization,Writing -original draft,Writing -review&editing.Fang-Yuan Shi:Methodology,Software,Data curation,Formal analysis.Yu Liang:Data curation.Ge Gao:Conceptualization,Project administration,Supervision,Funding acquisition,Resources,Writing-review&editing.All authors have read and approved the final manuscript.

Competing interests

The authors declare no competing interests.

Acknowledgments

This work was supported by funds from the National Key R&D Program of China(Grant No.2016YFC0901603),the National High Technology Research and Development Program of China (Grant No.2015AA020108),and the State Key Laboratory of Protein and Plant Gene Research and the Beijing Advanced Innovation Center for Genomics(ICG)at Peking University,China.The research of Ge Gao was supported in part by the National Program for Support of Top-notch Young Professionals.Part of the analysis was performed on the Computing Platform of the Center for Life Sciences of Peking University and was supported by the High-performance Computing Platform of Peking University.

Supplementary material

Supplementary data to this article can be found online at https://doi.org/10.1016/j.gpb.2021.06.001.

ORCID

0000-0003-2799-0369 (Yu Wang)

0000-0003-4185-8129 (Fang-Yuan Shi)

0000-0003-0695-2304 (Yu Liang)

0000-0001-6470-8815 (Ge Gao)