SPsimSeq is a semi-parametric simulation method to generate bulk and single-cell RNA-sequencing data. It is designed to simulate gene expression data with maximal retention of the characteristics of real data. It is reasonably flexible to accommodate a wide range of experimental scenarios, including different sample sizes, biological signals (differential expression) and confounding batch effects

Assefa, Alemu Takele

Thas, Olivier

Vandesompele, Jo

English

Ghent University Academic Bibliography

Gene expressionSPsimSeq: semi-parametric simulation of bulk andsingle-cell RNA-sequencing dataAlemu Takele Assefa 1,*, Jo Vandesompele2,3,4 and Olivier Thas1,3,5,61Data Analysis and Mathematical Modeling, 2Biomolecular Medicine, 3Cancer Research Institute Ghent, 4Center for Medical Genetics,Ghent University, Ghent, Belgium, 5National Institute for Applied Statistics Research Australia (NIASRA), University of Wollongong,Wollongong, Australia and 6Data Science Institute, I-BioStat, Hasselt University, Hasselt, Belgium*To whom correspondence should be addressed.Associate Editor: Inanc BirolReceived on June 19, 2019; revised on January 2, 2020; editorial decision on February 8, 2020; accepted on February 11, 2020AbstractSummary: SPsimSeq is a semi-parametric simulation method to generate bulk and single-cell RNA-sequencingdata. It is designed to simulate gene expression data with maximal retention of the characteristics of real data. It isreasonably flexible to accommodate a wide range of experimental scenarios, including different sample sizes, bio-logical signals (differential expression) and confounding batch effects.Availability and implementation: The R package and associated documentation is available from https://github.com/CenterForStatistics-UGent/SPsimSeq.Contact: alemutakele.assefa@ugent.beSupplementary information: Supplementary data are available at Bioinformatics online.1 IntroductionThe number of computational tools for the analysis of bulk andsingle-cell RNA-sequencing (RNA-seq) data is growing rapidly(Zappia et al., 2018). Several methods have been introduced for asingle task, e.g. testing for differential expression (DE). These toolstypically pass through an evaluation process, often focusing on falsediscovery rate control and sensitivity. Although such an evaluationoften relies on simulated data with a built-in truth, to realistically as-sess the performance of these data analysis tools, the simulated datamust faithfully recapitulate the data characteristics of real data(Soneson and Robinson, 2018; Weber et al., 2019).Various methods have been proposed for simulating either bulk orsingle-cell RNA-seq data. The starting point is typically a distribution-al assumption of the gene expression levels, e.g. the (zero-inflated)negative binomial distribution (Zappia et al., 2017). Although theseparametric simulation methods are flexible and allow simulating vari-ous scenarios by generating synthetic data with good fit to the realdata (Soneson and Robinson, 2018), such strong distributionalassumptions do not hold in general. Due to the intrinsic biologicalvariability and technical noise, single-cell RNA-seq data sometimesshow multimodal distributions (Bacher and Kendziorski, 2016).There are also fully non-parametric approaches that employ subsam-pling from real data (Benidt and Nettleton, 2015). Although non-parametric simulators generate realistic synthetic data, they havelimited flexibility and require a large source dataset to subsamplefrom Assefa et al. (2018) and Benidt and Nettleton (2015).Here, we present a new simulation procedure for simulatingbulk and single-cell RNA-seq data. It is designed to maximallyretain the characteristics of real RNA-seq data with reasonable flexi-bility to simulate a wide range of scenarios. In a first step, the loga-rithmic counts per millions of reads (log-CPM) values from a givenreal dataset are used for semi-parametrically estimating gene-wisedistributions and the between-genes correlation structure. In par-ticular, the estimation of the probability distributions uses the fastlog-linear model-based density estimation approach developed byEfron and Tibshirani (1996) and Lindsey (1974). The method makesuse of the Gaussian-copulas (Cario and Nelson, 1997) to retain thebetween-genes correlation structure, as implemented by Hawinkelet al. (2019) for parametric microbiome data simulation. Arbitrarilylarge datasets, with realistically varying library sizes, can be sampledfrom these distributions while maintaining the correlation structurebetween the genes. Our method has an additional step to explicitlyaccount for the high abundance of zero counts, typical for single-cellRNA-seq data. This step models the probability of zero counts as afunction of the mean expression of the gene and the library size(read depth) of the cell (both in log scale). Zero counts are thenadded to the simulated data such that the observed relationship(zero probability to mean expression and library size) is maintained.In addition, our method simulates DE by separately estimating thedistributions of the gene expression from the different populations(e.g. treatment groups) in the source data, and subsequently sam-pling a new dataset from each group.Our simulation procedure enables benchmarking of statisticaland bioinformatics tools with realistic simulated datasets. In Section3 and Supplementary Material, we demonstrate that the simulateddata from our method retains the characteristics of the source datain terms of variability, distribution of mean expression levels,VC The Author(s) 2020. Published by Oxford University Press. 1This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/),which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contactjournals.permissions@oup.comBioinformatics, 2020, 1–3doi: 10.1093/bioinformatics/btaa105Advance Access Publication Date: 17 February 2020Applications NoteDownloaded from https://academic.oup.com/bioinformatics/advance-article-abstract/doi/10.1093/bioinformatics/btaa105/5739438 by guest on 06 May 2020fraction of zero counts and the dependence between genes (Fig. 1and Supplementary Material). The details of the SPsimSeq proce-dures, implementations and benchmarking results can be found inthe Supplementary Material. Data simulated with our procedure arecompared with the original real source data and with data simulatedwith the parametric Splat procedure (Zappia et al., 2017), whichuses a gamma-Poisson hierarchical model (splatter R Bioconductorpackage, v1.6.1; Zappia et al., 2017).2 DatasetFor the demonstration and benchmarking of the SPsimSeq method,we used three publicly available datasets: (i) neuroblastoma bulkRNA-seq data retrieved from Zhang et al. (2015; GEO accessionGSE49711), (ii) neuroblastoma NGP cells single-cell RNA-seq dataretrieved from Verboom et al. (2019; GEO accession: GSE119984)and (iii) peripheral blood mononuclear cell single-cell RNA-seq dataretrieved from (www.10xGenomics.com). The details of the datasetscan be found in the Supplementary Material.3 ResultsUsing the three source RNA-seq datasets (one bulk and two single-cell), we benchmarked the novel SPsimSeq simulation method. Inparticular, we compared the simulated data (using SPsimSeq andSplat) with the real data with respect to various gene and sample(cell) level characteristics as suggested by Soneson and Robinson(2018) and Zappia et al. (2017). To simulate bulk RNA-seqdata using Splat, we disabled its feature for adding dropouts(dropout.type ¼ “none”), which is specifically designed forsingle-cell RNA-seq data simulation.The results generally show that our simulation procedure suffi-ciently captures the properties of the real data both for bulk andsingle-cell RNA-seq datasets (Fig. 1 and Supplementary Figs S1–S14). The variability (in terms of variance and coefficients of vari-ation), distribution of mean expression level of genes, the fraction ofzero counts (per gene and sample/cells), the relationship between themean and variability gene expressions, the relationship between themean expression and fraction of zero counts and the dependence be-tween genes in SPsimSeq simulated data resemble that of the realdatasets. When compared with Splat, SPsimSeq generates more real-istic data with respect to all the considered metrics. In theSupplementary Material, we present the detailed benchmarkingresults, including the application of SPsimSeq for simulating single-cell RNA-seq data with read-counts and UMI-counts (unique mo-lecular identifier).FundingThis work was supported by the UGent Special Research Fund ConcertedResearch Actions [GOA grant number BOF16-GOA-023].Author’s contributionsA.T.A.: Conceptualization, Data curation, Formal analysis, Methodology,Software, Visualization, Writing—original draft and Writing—review andediting. J.V.: Conceptualization, Funding acquisition, Resources, Supervisionand Writing—review and editing. O.T.: Conceptualization, Funding acquisi-tion, Methodology, Supervision and Writing—review and editing.Conflict of Interest: none declared.ReferencesAssefa,A.T. et al. (2018) Differential gene expression analysis tools exhibitsubstandard performance for long non-coding RNA-sequencing data.Genome Biol., 19, 96.Bacher,R. and Kendziorski,C. (2016) Design and computational analysis ofsingle-cell RNA-sequencing experiments. Genome Biol., 17, 63.Benidt,S. and Nettleton,D. (2015) SimSeq: a nonparametric approach to simu-lation of RNA-sequence datasets. Bioinformatics, 31, 2131–2140.Cario,M.C. and Nelson,B.L. (1997). Modeling and generating random vectorswith arbitrary marginal distributions and correlation matrix. TechnicalFig. 1. (A) Distribution of the coefficients of variations (CV) from the real and simulated (SPsimSeq and Splat) bulk RNA-seq datasets. (B) The relationship between the meangene expression levels (in log-CPM) and three other gene-level characteristics (fraction of zeroes, variance and CV) from the real and simulated single-cell RNA-seq datasets(read-counts). The right panels show the grid differences between the simulated and real data with respect to the corresponding gene-level characteristic. In particular, themean log2-CPM is partitioned into 1000 grids and the average characteristic per grid is used to calculate the differences. (C) The cumulative distribution of the fraction of zerocounts per gene from the real and simulated single-cell RNA-seq datasets (read-counts)2 A.T.Assefa et al.Downloaded from https://academic.oup.com/bioinformatics/advance-article-abstract/doi/10.1093/bioinformatics/btaa105/5739438 by guest on 06 May 2020report. Department of Industrial Engineering and Management Sciences,Northwestern University, Evanston, Illinois.Efron,B. and Tibshirani,R. (1996) Using specially designed exponential fami-lies for density estimation. Ann. Stat., 24, 2431–2461.Hawinkel,S. et al. (2019) A broken promise: microbiome differential abun-dance methods do not control the false discovery rate. Brief. Bioinform., 20,210–221.Lindsey,J. (1974) Construction and comparison of statistical models. J. R.Stat. Soc. B, 36, 418–425.Soneson,C. and Robinson,M.D. (2018) Towards unified quality verification ofsynthetic count data with countsimQC. Bioinformatics, 34, 691–692.Verboom,K. et al. (2019) SMARTer single cell total RNA sequencing. NucleicAcids Res., 47, e93.Weber,L.M. et al. (2019) Essential guidelines for computational methodbenchmarking. Genome Biol., 20, 125.Zappia,L. et al. (2017) Splatter: simulation of single-cell RNA sequencingdata. Genome Biol., 18, 174.Zappia,L. et al. (2018) Exploring the single-cell RNA-seq analysislandscape with the scRNA-tools database. PLoS Comput. Biol., 14,e1006245.Zhang,W. et al. (2015) Comparison of RNA-seq and microarray-based mod-els for clinical endpoint prediction. Genome Biol., 16, 133.SPsimSeq 3Downloaded from https://academic.oup.com/bioinformatics/advance-article-abstract/doi/10.1093/bioinformatics/btaa105/5739438 by guest on 06 May 2020

SPsimSeq : semi-parametric simulation of bulk and single-cell RNA-sequencing data

A Summary: SPsimSeq is a semi-parametric simulation method to generate bulk and single-cell RNA-sequencing data. It is designed to simulate gene expression data with maximal retention of the characteristics of real data. It is reasonably flexible to accommodate a wide range of experimental scenarios, including different sample sizes, biological signals (differential expression) and confounding batch effects.This work was supported by the UGent Special Research Fund Concerted Research Actions [GOA grant number BOF16-GOA-023].Assefa, AT (corresponding author), Univ Ghent, Data Anal & Math Modeling, Ghent, Belgium.
alemutakele.assefa@ugent.b

THAS, Olivier

Document Server@UHasselt (Universiteit Hasselt)

SPsimSeq: semi-parametric simulation of bulk and single-cell RNA-sequencing data

Document Server@UHasselt

https://biblio.ugent.be/publication/8660778/file/8660875.pdf

SPsimSeq : semi-parametric simulation of bulk and single-cell RNA-sequencing data

Abstract

Similar works

Full text

Available Versions

Ghent University Academic Bibliography

Document Server@UHasselt (Universiteit Hasselt)

Document Server@UHasselt