BACKGROUND
cell lines are widely used for perturbation experiments that aim to understand disease mechanisms at a cellular level.
cells used in such experiments are only rarely an inherent biological model for the disease of interest.
most commonly, genetic or environmental perturbations are required to create cellular responses, which can then serve as an experimental disease model.
it is known that many cell lines carry major genetic variations, which would be lethal for humans at the stage of prenatal development.
the main advantage of highly proliferative cell lines relative to primary cells and induced pluripotent stem cells is a significantly greater capacity for experiments that require large amounts of clonal cells with identical genetic background, such as needed in state-of-the-art high-throughput screening, especially in proteomics or metabolomics.
the availability of clonal cells allows the possibility for comparative perturbation experiments aiming to compare phenotypic outputs derived from a set of single node perturbations  <cit> .
in addition, using cell lines avoids the ethical concerns arising out of human primary neuronal cell culture.
the human neuroblastoma cell line sh-sy5y first described in  <cit> , is a commonly used cell line in studies related to neuroblastoma and neurodegenerative disease.
the cell line is a sub-clone of the parent cell line sk-n-sh which was originally established from a bone marrow biopsy of a neuroblastoma patient  <cit> .
the three human diseases most frequently mentioned in literature for sh-sy5y were neuroblastoma, alzheimer’s disease , and parkinson’s disease  .
a complete genomic characterization of sh-sy5y would thus elucidate the applicability and possible limits for modelling these disease-specific processes in the context of this cell line.
undifferentiated sh-sy5y cells have been extensively used as an in vitro model for research in neuroscience  <cit> .
the cell line shows biochemical properties of immature catecholaminergic neurons  <cit> .
studies have found that undifferentiated sh-sy5y express only immature neuronal markers and lack mature neuronal markers  <cit> .
thus, undifferentiated sh-sy5y cell line might not represent an appropriate model for diseases such as pd, which primarily affect differentiated dopaminergic neurons  <cit> .
further, xie et al.
<cit>  reviewed  <dig> articles on sh-sy5y cell as in vitro model for pd research and found that differentiation of sh-sy5y under a certain treatment results in a more dopaminergic neuronal phenotype, which could be extremely useful for modelling selective dopaminergic cell death in pd.
however, they point out that an optimally differentiated sh-sy5y dopaminergic cell model requires further research.
a systems genomics analysis could thoroughly assess the genetic mutations in all neuronal markers and therefore place limits on how closely differentiation can model a dopaminergic cell.
furthermore, treatment with differentiation-inducing agents enable sh-sy5y cells to become morphologically similar to mature primary neurons  <cit> , and the different treatment agents  result in a variety of neuronal phenotypes   <cit> .
genomic mutations may, however, impose limitations on the possible phenotypes to which sh-sy5y can differentiate.
therefore, a complete genomic characterisation of the sh-sy5y cell line can inform the limitations on the possible phenotypes imposed by the genome.
systems genomics aims to integrate genomic variation, copy number and structural variation with high-throughput gene expression, metabolomics and proteomic data to explore the genetic architecture of complex traits and multi-factorial diseases.
for characterizing the sh-sy5y cell line we integrated information from whole genome sequencing, transcriptomics, and proteomics experiments .
firstly, whole-genome sequencing is used to determine the genetic background of the cell line.
secondly, transcriptomics and proteomics data is used to study correlation across biological levels.
finally, we integrate genomic variants using a network analysis approach to evaluate the suitability of the sh-sy5y cell line as an in vitro model to study various neurodegenerative diseases, including pd.figure  <dig>
systems genomics approach to assess sh-sy5y as a disease model.
RESULTS
whole genome sequencing
the whole genome of sh-sy5y cell line was re-sequenced using two different sequencing platforms, those from complete genomics  and illumina  hiseq <dig>  to generate a high-quality list of variants.
an earlier study has shown the utility of combining information from both platforms and although they show high concordance, each platform alone failed to detect a significant number of exonic variants spread over  <dig>  genes  <cit> .
for the sh-sy5y cell line, dna sequencing by cg  <cit>  and il  <cit>  produced a genome-wide coverage of 57× and 49× respectively.
more than 90% of the genome was covered by a minimum read depth of  <dig> .
further, more than 95% of the exome was called with minimum required confidence set by each of the proprietary genotype-calling procedures of the two respective sequencing platforms.
the two genome-sequencing platforms – cg and il – produced a total union of  <dig> , <dig> single nucleotide variations  with 84% concordance between both platforms.
filtering them based on criteria defined in reumers et al.
<cit>  resulted in  <dig> , <dig> snvs with 99% concordance between both platforms.
out of these high quality variants 96% were previously reported by dbsnp build  <dig>  <cit>  and  <dig> % by  <dig> genomes project  <cit> .
<dig>  snvs and small indels were earlier reported by the catalogue of somatic mutations in cancer   <cit>  version  <dig>  the total number of snvs and the percentage of snvs found in  <dig> genomes project are similar to controls  in another study  <cit>  taken from the human genome diversity project  <cit> .
among the high confidence set of mutations in sh-sy5y,  <dig> were confirmed as somatic in the cosmic database and were also rare in the general population .table  <dig>
high-confidence mutations in sh-sy5y that were rare in the population and also confirmed as somatic in the cosmic database
rare mutations are snvs and small indels that were found in less than 5% of the samples in  <dig> genomes project, exome sequencing project and complete genomics baseline genomes.
in the genes column, when distances are given, the mutations are found in intergenic regions and the first gene precedes the mutation whereas the second gene succeeds the mutation.
the number of private protein-altering  snvs and indels, i.e.
protein-altering variants not found in  <dig> genome project or the exome sequencing project  <cit>  were  <dig> and  <dig> respectively , whereas the corresponding numbers for controls were  <dig>  and  <dig> .
the set of control genomes refers to  <dig> genomes used in another study  <cit>  taken from the human genome diversity project  <cit> .
the reason for the discrepancy between the number of ppa snvs between sh-sy5y and the control genomes is mainly due to our strict filtering strategy as well as the different sequencing platforms that were used .
the overlap between ppa variants and the variants in the catalogue of somatic mutations in cancer  was  <dig> and  <dig>  for sh-sy5y and control genomes respectively.
the overlap between genes containing ppa snvs and indels and the sanger cancer gene census   <cit>  was  <dig>  similar to controls .
we also performed a gene ontology  enrichment analysis on the  <dig> genes containing filtered ppa snvs and indels and also on the  <dig>  genes using unfiltered ppa snvs and indels.
however, no go term was significantly enriched .
we judged the overall quality of the variant calls and also compared them between the two genome sequencing platforms - cg and il.
then, we integrated variant calls from both these platforms by using the filtering criteria provided by an earlier study  <cit>  that compared these two platforms .
validation of snvs and small indels
validation was performed using the illumina omni- <dig> quad microarray, which assays loci from hapmap phase 1– <dig>  <cit>  and  <dig> genomes project.
out of  <dig>  heterogeneous snvs that were queried by genotyping,  <dig> % were concordant, which increased to  <dig> % for the filtered variants .
functional prediction of snvs and small indels
the functional effect of snvs and small indels were predicted using annovar  <cit>  by annotation with labels for genomic regions  and coding effects .
the annotation was performed by taking a consensus across four databases  based on choosing the most damaging effect predicted in order to aim for sensitivity .
in general, the association of snvs to genomic regions for platform-specific and concordant snvs did not show any significant differences.
the association of snvs with exonic, intronic and intergenic regions was 1%, 26–36% and 46–50% respectively, similar to the results in lam et al.
<cit>  .
the functional prediction found  <dig> genes with rare non-synonymous snvs, indels or substitutions .
neuroblastoma-relevant genes
the genome sequencing of sh-sy5y found  <dig> genes with rare non-synonymous snvs, indels or substitutions that overlapped with the list of  <dig> genes containing somatic mutations in the complete genome sequence of  <dig> untreated primary neuroblastoma tumours  <cit>  .
using the hyper-geometric test, genes with rare non-synonymous snvs and indels were found to be significantly enriched among the genes with somatic mutations in primary neuroblastoma tumours .
out of these  <dig> genes, only mutations in  <dig> genes were predicted as damaging by sorting intolerant from tolerant   <cit> .
furthermore, only  <dig> genes with rare non-synonymous snvs, indels and substitutions overlapped with the  <dig> genes associated with neuroblastoma in the literature  and therefore were not significantly enriched .
however, genes with non-synonymous snvs, indels and substitutions, copy number variations, and structural variations  were significantly enriched among the  <dig> genes extracted for neuroblastoma from literature with an overlap of  <dig> genes .
we also compared genomic mutations in sh-sy5y to other genomic sequencing studies of primary neuroblastoma from patients.
among the seven genes mutated at a significant frequency in a study of  <dig> matched tumour -normal samples  <cit> , none of the genes carried a rare amino-acid changing mutation.
among the rare germ-line variants predisposing to neuroblastoma identified by the target study, palb <dig> contained two rare non-synonymous mutations in sh-sy5y.
among the  <dig>  coding somatic mutations found in an aggregate of  <dig> matched tumour/normal samples, we found  <dig> overlapping variants  in sh-sy5y inside  <dig> genes among which  <dig> were rare amino-acid changing mutations inside  <dig> genes .
the pediatric cancer genome project  conducted a study of  <dig> patients with metastatic neuroblastoma and found mutations in atrx and alk in 22% and 14% of the patients correspondingly  <cit>
.
however, sh-sy5y contained no rare amino-acid-changing mutations in these genes.
structural variations  and mobile element insertions
structural variation events consist of large deletions , duplications , inversions, translocations and complex variations - other combinations of such chromosomal rearrangements.
the majority of sv events in the sh-sy5y cell line found by cg were deletions .
around 2–3% of cancers show chromothripsis  <cit> , where tens to hundreds of genomic rearrangements occur in a cellular crisis event.
in a certain region, a high density of genomic rearrangements  combined with frequent oscillations between two copy number states  <cit>  and frequent occurrences of runs of homozygosity constitute the hallmark of chromothripsis.
cytogenetic methods confirm a high density of genomic breakpoints in single cells and indicate that these breakpoints are not a result of parallel rearrangements in different sub-clones  <cit> .
however, for sh-sy5y, dna sequencing found little to no evidence of chromothripsis.
most notably, no oscillations between copy number states of 100-kb binned read depth were observed except for some on chromosomes  <dig>   <dig>   <dig> and x .
in all of these four chromosomes, the regions with oscillating copy number were near the centromeres  or the telomere .
even though all of these regions have at least six junctions, these regions consist of highly repetitive dna  <cit>  and are therefore prone to mapping errors.
we filtered the svs based on their frequency in the cg baseline genome set  <cit>  and extracted high-confidence calls with a frequency of less than 10% in the baseline.
the cg baseline set comprises of  <dig> genomes of healthy, disease-free individuals, which could be used to filter technical artifacts and variations common in the population.
after filtering, the high-confidence deletions overlapped with  <dig> genes including the following cancer-related genes - pten , which is a key tumour suppressor gene  <cit> , ctnna <dig> , which is a cell contact inhibition gene whose mutation can promote cancer development and formation, mcc  and mtus <dig> .
mobile element insertions  refer to insertion of sequences that can change their position within the genome.
<dig>  mei events were detected in the cg data  and the majority of those were alu events  and l <dig> retro-transposons .
copy number variation
cg identifies discrete coverage levels corresponding to ploidy levels using the distribution of observed normalized coverage values and then uses a hidden markov model to assign copy number levels .
since cg essentially uses the relative coverage as a proxy for copy number, it cannot directly infer the absolute copy number levels.
however, m-fish  experiments have demonstrated that diploidy dominates the copy number level in the sh-sy5y genome  <cit> .
therefore, we calculated the absolute copy number levels assuming that the most commonly found copy number level is  <dig>
we compared the copy number levels found by cg with an earlier measurement using comparative genomic hybridization  arrays  <cit> .
the cgh arrays used had a lower limit one million base pairs to their resolution, a drawback not found in second-generation sequencing technologies such as cg and il.
to achieve even higher resolutions and reduce false positive and false negative copy number determinations, we normalized the coverage levels to those observed in a reference set of  <dig> in-house genomes also sequenced by cg.
a visual illustration of the copy number levels  confirms key features of the cell line described previously such as partial trisomy of chromosome  <dig>  gain of chromosome  <dig>  2p, 17q, and loss of 1p, 14q and 22q.
among the large-scale copy number features identified in the pcgp project, the gain of 17q  was common with sh-sy5y.
the remaining cnvs that were reported with a lower frequency in pcgp - loss of 1p  or loss of 11q  were not found in sh-sy5y.figure  <dig>
cnvs detected by whole genome sequencing and array-based cgh by do et al.
<cit> .
the results from whole genome sequencing were from complete genomics and are given in the left half of the chromosomes.
the results from the array-based cgh are on the right half of the chromosomes.
regions are highlighted for copy number gain  and loss .
the major events partial trisomy of chromosome  <dig> and  <dig>  complete trisomy of chromosome  <dig>  gain in 17q and loss in 22q were confirmed.
.
interestingly, the strongest clinical marker of neuroblastoma that also determines the aggressiveness of neuroblastoma  <cit>  - gain of the gene mycn  – actually had a relative coverage of only  <dig>  however, more recently,  <cit>  confirmed the finding that the mycn gene is not strongly amplified in sh-sy5y.
also, among the  <dig> patients sequenced in the pcgp project, only 23% carried an amplification of mycn, which was defined as > <dig> copies detect by quantitative-pcr.
transcriptome analysis
to characterize the sh-sy5y cell line gene expression, we sequenced poly-adenylated rna from undifferentiated sh-sy5y cells using a stranded rna-seq protocol.
we used the grch <dig> human genome version together with the ensembl release  <dig>  <cit>  gene annotation.
we generated about  <dig> million paired-end reads of lengths 50 nt.
the raw reads are available in the european nucleotide archive  database  <cit>  under the study accession prjeb <dig>  in total,  <dig> million paired-end reads and  <dig>  million singleton reads could be mapped after quality control  filtering and trimming to either the genome or the transcriptome annotation.
this yields in 95% mappable reads to either genome or cdna sequences.
for rna-seq, we called snvs and small indels using both samtools  <cit>  and the genome analysis toolkit   <cit> .
we retained the intersection of the two outputs and those caller-specific variants supported with a read depth greater than or equal to  <dig>  as a recent study has shown low concordance between multiple variant-calling pipelines  <cit> , we decided to increase the level of confidence by using two variant-calling pipelines.
consequently for rna-seq, out of  <dig>  snvs and small indels detected,  <dig> % were concordant with unfiltered variants from at least one dna sequencing platform, either cg or il .
the rna-seq also showed high sensitivity as  <dig> % of the  <dig>  filtered exonic variants found through dna sequencing in genes with fpkm >  <dig> were also detected by rna-seq.
fpkm refers to fragments per kilobase of exon per million fragments mapped, which is a measure of gene expression.
this sensitivity of rna-seq for variant detection was higher than that achieved by a recent approach for snv calling in rna-seq data  <cit>  but the absolute number of detected exonic variants  was slightly lower than theirs , which might be due to our strict quality filtering of the dna sequencing data.
as for neuroblastoma-relevant genes from the  <dig> neuroblastoma primary tumours, out of the  <dig> rare, amino-acid changing mutations ,  <dig> of them were detected by rna-seq with the same zygosity.
the remaining  <dig> were likely not detected because of low gene expression as all mutations detected by rna-seq were in genes with fpkm >  <dig>  and all mutations not found by rna-seq were in genes with fpkm <  <dig> .
genetic copy number vs. gene expression
we compared genes that were always expressed in sh-sy5y across  <dig> different conditions  from the geo  database  <cit>  to those that were never expressed, whenever the corresponding probe was present.
here, the threshold used to determine whether genes were always expressed in the geo dataset was that the expression values of the corresponding transcripts were always greater than the median for all the transcripts in each microarray experiment.
similarly, the threshold for genes never expressed in the geo dataset was that the expression values were always less than the median for all the transcripts in each microarray experiment.
the frequency of mutations per  <dig> kilo-base-pair region in genes that were never expressed was  <dig> , which decreased to  <dig>  for genes, which were always expressed.
similarly, the average copy number for genes always expressed was  <dig> , which was lowered to  <dig>  for genes never expressed .
the distribution of expression of genes with copy number greater than  <dig> was significantly greater  than expression of genes with copy number lower than or equal to  <dig> .
further, we compared the genetic copy number and the rna-seq data of the same sh-sy5y sample.
the copy number and the logarithm of fragments per kilo-base of transcripts per million mapped reads  showed a positive correlation .
proteome analysis
in order to confirm and correlate the biological expression of the identified genes, the analysis of genomic data was integrated with the analysis of proteomics data analysis of sh-sy5y cell lysates .
proteins isolated from whole sh-sy5y cell lysates were fractionated by sds-page, in-gel digested using trypsin, and the recovered peptides analysed by lc-esi ms/ms.
the spectra were analysed by maxquant against a combinatorial human protein database.
this database contained the human genome reference hg <dig> protein entries combined with the sequence variants found from the dna sequencing done here.
specifically, the coding sequences of ensembl genes were modified using all the homozygous exonic variants .
in total, out of  <dig>  acquired ms/ms spectra,  <dig>  were matched to a peptide that could be mapped to a total of  <dig>  protein identifiers .
from these,  <dig>  protein groups corresponding to  <dig>  proteins were identified in the human genome reference hg <dig> protein annotation.
further,  <dig> proteins were identified using sequence variants found from the dna sequencing done here.
see additional file  <dig> for the abundances of proteins detected.
to increase peptide spectrum matches, we extended the human genome reference hg <dig> protein database by inserting reference sequences containing the homozygous snvs and indels in the exonic regions identified in this study.
the peptides that were additionally mapped to the extended protein database were then used for validation of genetic variants .
this approach validated  <dig> snvs,  <dig> indels and  <dig> block substitution.
gene expression vs. protein expression
gene expression levels from rna-seq and the corresponding protein abundance were compared to detect correlation.
we plotted the logarithm of fpkm against the logarithm of intensity–based absolute quantification   <cit>  score of  <dig>  genes , which were detected in both rna-seq and proteomics, which equals 93% of  <dig>  proteins detected by proteomics.
as expected, the correlation was weakly positive .
network-based analysis of suitability of sh-sy5y cell line as in vitromodel
in order to evaluate the suitability of sh-sy5y as in vitro model we adopted a framework described in  <cit>  that integrated genetic sequence information, and topology analysis of either disease or process-specific networks.
we applied betweenness-centrality ratio , a metric that allowed us to quantitatively assess the impact of genes mutated in the cell line on the disease or process network.
bc-ratio is normalized to lie between  <dig>  and  <dig> .
it indicates how much the information flow in the network might be altered by changes in the cell line: the higher bc-ratio, the greater the impact of the mutations, the lower the genetic integrity  of the cell line in the context of a specific network.
.
the bc-ratio metric applied to four neurodegenerative diseases  ranked alzheimer’s  as highest bc-ratio =  <dig> , pbc-value =  <dig> ), followed by parkinson’s disease , huntington disease , and amyothrophic lateral sclerosis .
to investigate the appropriateness of the cell line for perturbation experiments in the context of pd, we also scored the different processes corresponding to the hallmarks of pd .
see figure  <dig> for a protein-protein interaction network resulting in the cell line suitability scores for neuroblastoma and als and figure  <dig> for the network for glycolysis and reactive oxygen species  metabolism in pd.
additional file  <dig> contains network visualizations for each disease and each module in pd.
the role of the visualization is to give an idea of how damaged genes are distributed in the network, to which genes they are connected and how they could alter processes of interest.table  <dig>
network statistics and cell line scoring of the neurodegenrative diseases and neuroblastoma
flow bc stands for flow betweenness centrality.
the pbc-value refers to the network randomization test described in the methods section.
flow betwenness is defined for connected networks.
if the entire network is not connected, flow betweenness is computed for connected components.
in the table the flow betweenness centrality of the largest connected component is given and marked by *.
nb – neuroblastoma; ad – alzheimer’s disease; pd – parkinson’s disease; hd – huntington’s disease; als – amyotrophic lateral sclerosis.
network statistics and cell line scoring of the parkinson’s disease modules
flow bc stands for flow betweenness centrality.the pbc-value refers to the network randomization test described in the methods section.
flow betwenness is defined for connected networks.
if the entire network is not connected, flow betweenness is computed for connected components.
in the table the flow betweenness centrality value of the largest connected component is given and marked by *.
protein-protein interaction network for two diseases – neuroblastoma  and for als .
red nodes refer to omim-derived genes mutated in the cell line.
orange nodes refer to the omim-derived genes which are intact in the cell line.
dark green nodes represent genes coming from the network expansion, which are mutated in the cell line.
light green nodes represent genes coming from the network expansion, which are intact in the cell line.
nodes are scaled to the magnitude of their betweenness centralities.
blue edges show connections between pairs of nodes in which at least one was damaged in the cell line.
for neuroblastoma, mutations in the central genes – nme <dig> and alk contribute to the high bc-ratio.
for the als network, mutations occur in genes with lower betweenness centrality which results in the lower bc-ratio.
protein-protein interaction network for two pd map modules - glycolysis  and ros metabolism .
red nodes refer to genes mutated in the cell line.
orange nodes refer to genes which are intact in the cell line.
nodes are scaled to the magnitude of their betweenness centralities.
blue edges show connections between pairs of nodes in which at least one was damaged in the cell line.
for glycolysis, mutations in the central genes lead to the high bc-ratio.
mutated genes lie on the periphery of the ros metabolism network resulting in a low bc-ratio.
discussion
the complete system-wide omics analysis of the undifferentiated sh-sy5y cell line provided many opportunities for studying concordance across biological levels .
firstly, protein-altering variants, copy number variations and structural variations show consistency with the expression level of genes in the rna-seq data.
for instance, variant calling from rna-seq data showed  <dig> % sensitivity for the  <dig>  filtered variants detected from dna-seq in genes with an fpkm >  <dig>  secondly, sequencing both the dna and the rna showed a high degree of concordant snvs and indels  and served as further validation of dna sequencing results and our filtering procedure for genomic variants.
finally, the knowledge of genomic variants enhanced protein identification in proteomics experiments where it increased the number of proteins detected from  <dig>  to  <dig> .
a comparison of the two whole genome sequencing platforms – complete genomics  and illumina  – confirmed a significant number of variants that were discordant, either due to regions, which were not called by one of the platforms or due to platform-specific errors.
considering that 1%  of these platform-specific variants residing in  <dig>  genes were associated with exonic regions, using both platforms would increase the detection of potentially functionally important snvs.
also only 75–80% of platform-specific variants were found in the known list of variants , indicating that a significant number of mutations are somatic mutations or sequencing errors that need further investigation.
further, the larger platform-specific biases for small indels and structural variations also strengthened the argument for using both platforms for greater sequencing accuracy.
neuroblastoma-relevant genetic variations
without the genome sequence of the healthy cells from the patient, the difficulty of identifying somatic mutations in sh-sy5y precludes any direct comparison to the somatic mutations in primary neuroblastoma tumours.
thus, we compared rare mutations in sh-sy5y with the somatic mutations found in  <dig> primary neuroblastoma tumours.
the union of genes from  <dig> different patient samples resulted in a list of  <dig> genes, which we considered as relevant to neuroblastoma.
the result that only  <dig> genes with rare amino-acid changing mutations predicted as damaging in sh-sy5y overlapped with the list of  <dig> genes with somatic amino-acid changing mutations in primary neuroblastoma was expected as there was considerable heterogeneity among the different tumour samples themselves.
specifically, the sequencing of  <dig> neuroblastoma primary tumours found few recurrent genes with amino-acid changing mutations - only  <dig> genes contained amino-acid changing mutations in more than  <dig> tumour sample and only three genes in more than two tumour samples.
in this context, an overlap of  <dig> genes with rare amino-acid changing mutations suggests that sh-sy5y includes elements of the genetic architecture of a variety of neuroblastoma tumours.
therefore, it might offer opportunities to investigate traits arising from amino-acid changing mutations found in different kinds of neuroblastoma tumours.
as for structural variations, recurrently occurring structural alterations in the neuroblastoma primary tumours in ptprd, odz <dig> and csmd <dig> were not found among high-confidence rare structural alterations in sh-sy5y.
finally, chromothripsis, which was earlier identified in 18% of primary neuroblastoma, was not detected in sh-sy5y.
therefore, sh-sy5y would not serve as an appropriate model for neuroblastoma tumours suffering from chromothripsis.
suitability of sh-sy5y as an in vitromodel for neurodegenerative diseases
the bc-ratio relies on the network betweennness centrality to quantify genetic changes in the cell line.
it is important to consider that the bc-ratio does not aim to interpret the functional effects of these mutations.
indeed, mutations could cause a loss, gain, change or no change in the function.
however, such interpretations require assumptions.
in the case of neuroblastoma, we know that the sh-sy5y cell line was derived from a neuroblastoma patient.
hence, it can be assumed that many genetic mutations, which were observed in the neuroblastoma network , are indeed typical neuroblastoma mutations.
specifically, copy number variations in the central neuroblastoma genes nme <dig> and alk , both of which are found to be mutated in neuroblastoma as well support this reasoning.
briefly, it can be assumed that sh-sy5y cells are a good model for neuroblastoma cells and a high bc-ratio  was expected.
however, this case where a neuroblastoma cell line is expected to model neuroblastoma is a special case.
this interpretation is similar to the approach of domcke et al.
<cit>  of choosing the ‘most suitable models’ for ovarian tumour among ovarian cancer cell lines based on genomic and mrna expression profiles.
typically, in contrast, experimentalists use cell culture models to study processes lacking such ideal identity matching.
routinely, perturbation experiments are used to cause controllable genetic defects in networks, which are overall assumed to be intact  <cit> .
however, potential genetic changes, which are ignored, can introduce errors into the interpretation of results.
the main scope of the bc-ratio is to provide distances between networks of interest in a cell line and their ideal reference networks.
a bc-ratio of  <dig>  indicates an unchanged network representing an ideal reference for perturbation experiments.
the pbc value evaluates if a bc-ratio >  <dig>  was due to chance or it resulted from strong correlation between the nodes with high betweenness centrality and genes mutated in the cell line.
note that a non-significant pbc-value does not exclude any functional effects of network changes and it is strongly recommended to collect additional information for evaluating potential functional effects.
while our aim was not to predict functional effects related to rna or protein abundance, we are convinced that the provided data will support case-specific interpretations of functional network integrity.
a central question of our investigation was whether context -specific networks are overall changed in sh-sy5y.
the absence of network changes was considered as ideal case for controllable perturbation experiments.
in the following paragraphs, we first discuss the integrity of disease-specific networks in sh-sy5y and secondly, focus on pd in more detail.
the cell line suitability  metric applied to four neurodegenerative diseases ranked them from those with the largest to the smallest changes in the network.
alzheimer’s disease  was ranked as the highest, followed by parkinson’s disease , huntington’s disease  and, finally, amyotrophic lateral sclerosis .
bc-ratios in combination with the pbc-values tell us that genes mutated in the cell line are more likely to be involved in the information flow control in ad  network than the same type of genes is in the networks of other diseases .
overall, except for the neuroblastoma, bc-ratios were rather low , and pbc-values were not significant for all but neuroblastoma and alzheimer’s disease.
these results indicate that the networks related to pd, hd and als are largely intact in the sh-sy5y cell line and the latter can be considered a suitable model for perturbation experiments targeting these diseases, while ad might require a special caution.
it is also important to bear in mind that none of these disease-specific networks were perfectly intact.
indeed, the bc-ratio for each network was computed as greater than  <dig> , which indicates that changes in these networks were detected.
a prediction of potential functional effects of these network changes, including a classification in gain or loss of function goes beyond the scope of the bc-ratio but the provided data on rna and protein abundances can support such interpretations.
genes from which the networks for the neurodegenerative diseases and neuroblastoma have been grown, were derived from the omim databse where they have been classified as having disease causative effect .
we also calculated the bc-ratios of the respective networks assuming that only these seed genes have had been damaged in the cell line.
in all cases the number of mutated genes decreased while the bc-ratios of the cell lines increased for all but neuroblastoma .
note that in the case of neuroblastoma, the main actors which contributed to the high bc-ratio computed on the real data, already were omim genes, damaged in the cell line.
therefore a lower bc-ratio for neuroblastoma was expected in the absence of any other damaged genes.
the bc-ratio of the huntington’s disease increased considerably.
this is mostly due to the gene htt followed by prnp, which have the highest betweenness centrality in the hd network.
alzheimer’s disease scores next with app gene making the main contribution to the magnitude of the bc-ratio, being the first on the list of all genes.
it is immediately followed by  <dig> more omim-derived genes while their betweenness centrality is more than two times smaller.
similar effect can be observed in the case of amyotrophic lateral sclerosis: sod <dig> has the highest betweenness centrality, followed by vcp and dctn <dig> with almost twice as much lower betweenness.
in this experiment, parkinson’s disease had the lowest bc-ratio although the list was dominated by  <dig> omim-derived genes.
the explanation might lie in the low betweenness centralization of the network .
network centralization  <cit>  reflects how much variation is there in the centrality scores  among the nodes.
the value of  <dig>  indicates that the scores are rather equally distributed among the nodes and there are no nodes with apparent brokering role.
therefore changing the node labels from intact to damaged in such network produces smaller effect than in the remaining diseases.table  <dig>
bc-ratios of the neurodegenerative diseases and neurobalstoma computed under assumption, that only omim-derived seed genes have been mutated
the increase in the bc-ratios of all neurodegenerative diseases with the causative omim genes being considered as damaged, suggests that the genes which indeed are mutated in the cell line have smaller positional advantages than omim genes in the corresponding networks and therefore might less be involved in the information flow control.
we believe that this observation supports our conclusion about the overall sh-sy5y cell line suitability for the experimental studies of the neurodegenerative diseases.
to further analyse the suitability of the cell line as an in vitro model for pd, we studied the integrity of disease related sub-networks including mitochondrial dysfunction, reactive oxygen species  accumulation and calcium homeostasis  <cit> .
fujita et al.
<cit>  recently provided a comprehensive pd map, which integrates metabolic reactions, gene regulation and signalling processes in this complex network.
the network consists of sub-networks , which correspond to pathways and processes.
in these networks, nodes represent genes, proteins or small molecules, and edges represent molecular interactions.
we selected seven modules corresponding to hallmarks of pd  <cit>  from the pd map and scored them according to the bc-ratio metric.
the idea was to segregate pd-related modules as more suitable for studying gain or loss of function of a module.
in descending order of bc-ratio, the modules were glycolysis, mitochondria, calcium signalling, apoptosis, dopamine metabolism, ubiquitin proteasome system and reactive oxygen species  metabolism .
the pbc values indicated no significant changes in these networks.
however, we found that the glycolysis module  is the most impacted , which suggests that information flow through this module could be significantly altered by the genetic background of sh-sy5y.
one explanation for glycolysis being highly impacted could come from the fact that sh-sy5y is a cancer cell line.
it has been shown that glycolytic rates are higher in cancer cells than healthy ones  <cit> , known as the warburg effect  <cit> , which has been suggested to confer a proliferative advantage to tumour cells.conversely, the ros metabolism module  was the least impacted , which would indicate that the genes, mutated in the cell line do not compromise the information flow in this network.
for the functional interpretation beyond the bc-ratio, it is important to consider that most of the genetic variations found in the pd modules do not carry mutations, but almost exclusively copy number gains.
indeed, this gain of copy numbers is expected because of the partial trisomy of chromosomes  <dig>   <dig>  and  <dig> and the full trisomy of chromosome  <dig>  however, variations that alter the protein sequence occur in very few genes.
among all the damaged genes found in the seven modules from the pd-map, only four contained non-synonymous snvs.
the sorting intolerant from tolerant   <cit>  tool predicted only  <dig>  as damaging .
the gene myo <dig>  found in the calcium signalling module and involved in intracellular vesicle transport, contains a heterozygous non-synonymous mutation with a sift score of  <dig> .
it also plays an important role in trafficking and activity-dependent recruitment of ampa receptors to synapses  <cit> .
the gene hadh, found in the mitochondria module and which catalyses several reactions in beta-oxidation, also contains a heterozygous non-synonymous mutation with a sift score of  <dig>  the gene gzmb, found in the apoptosis module and which codes for a serine protease that is used by activated cytotoxic t lymphocytes to induce cell apoptosis  <cit> , contains a rare non-synonymous snv with a sift score of  <dig> .
however, in each of these three genes, only one copy of the gene contains such rare non-synonymous mutations.
further, four genes  in the mitochondria module and one gene  in the apoptosis module had a copy number of one, leaving only a single working copy.
ndufa <dig> has been shown to have a high degree of nitration and be associated with the oxidative damage to mitochondrial complex i  <cit> .
tomm <dig> plays an important role in mitochondrial clearance controlled by pink <dig> –park <dig> pathways  <cit> .
atp5l <dig> is involved in hydrogen ion trans-membrane transport activity  <cit> .
xrcc <dig> codes for the ku <dig> protein, which repairs double bond breakage of dna through the non-homologous end-joining pathway.
therefore, one should consider the damage to these genes when designing loss of function perturbation experiments.
indeed, it is important to test if the mutations per se cause a loss of function and to keep in mind that a function which is already lost cannot be lost a second time.
except for the  <dig> genes discussed above, all the genes in the modules of the pd map were free of translated mutations in sh-sy5y.
genes related to dopamine metabolism in sh-sy5y
in summary, all the genes involved in the dopamine metabolism pathway  except cyp2d <dig>  adh1b and ugt1a <dig> were expressed with or without specific treatments.
further, cyp2d <dig>  which contributes to dopamine biosynthesis through an alternative cytochrome p450-mediated pathway shown to exist in rats  <cit> , has a copy number of one, which may contribute to the lower expression.
among the genes in the dopamine metabolism pathway that were expressed with or without specific treatments, only adh1c and dbh contained a rare protein-altering mutation and ddc and akr1b <dig> had a copy number of three.
only the mutation in adh1c was predicted as damaging by sift.
therefore, the large majority of genes involved in dopamine metabolism contains no genetic damage and can be expressed under certain treatments.figure  <dig>
dopamine biosynthesis and degradation and the genes coding key enzymes.
this metabolic pathway is a modified version of a figure from meiser et al.
<cit> .
in addition to the enzymes involved, the genes coding for those enzymes have been added using the humancyc  <cit> .
the mutations in sh-sy5y affecting enzyme-encoding genes have been annotated only if they were rare protein-altering snvs or indels, cnvs, or svs.
CONCLUSIONS
in summary, we provide here the first whole genome analysis of the sh-sy5y cell line, which is widely used as a model for various neurological diseases.
we characterise the different types of biological information of sh-sy5y – genomics, transcriptomics and proteomics - and compare the relationships between them in terms of expression levels and variants.
the data related to the genome sequencing, rna sequencing and proteomics can be found at the following url - http://systemsbiology.uni.lu/shsy5y.
additionally, we also compare the sh-sy5y genome using two widely used whole genome sequencing platforms and show that using the two platform-specific coverages and sensitivities to different indel sizes supplements each other.
finally, our analysis based on context-specific network integrity ranked the integrity of sh-sy5y for parkinson’s disease as higher than for alzheimer’s disease but lower than for amyotrophic lateral sclerosis and huntington’s disease.
in the context of different pd related sub-systems, the same approach ranked ros metabolism as the most intact followed by the ubiquitin proteasome system, dopamine metabolism, calcium signalling, mitochondria and glycolysis.
overall, most of the genes belonging to the major parkinson’s disease pathways and modules were intact in the sh-sy5y genome.
particularly, each analysed gene related to pd has at least one intact copy in sh-sy5y.
therefore, somatic mutations do not significantly alter pd-related pathways in sh-sy5y.
importantly, the systems genomics analysis was performed on the undifferentiated cell line although differentiated cell lines are considered better suited as in vitro models for pd.
our novel bc-ratio method for cell line scoring integrates a genomic characterisation of the cell line, core set of disease or process-related genes and a protein-protein interaction network.
this method is not restricted to a specific cell line or disease and can be broadly applied.
furthermore, if multiple sequenced cell lines are available, the bc-ratio scores can also be used as a guide to select the cell line best suited for studying a particular disease.

 cell lines are widely used for perturbation experiments that aim to understand disease mechanisms at a cellular level. cells used in such experiments are only rarely an inherent biological model for the disease of interest. most commonly, genetic or environmental perturbations are required to create cellular responses, which can then serve as an experimental disease model. it is known that many cell lines carry major genetic variations, which would be lethal for humans at the stage of prenatal development. the main advantage of highly proliferative cell lines relative to primary cells and induced pluripotent stem cells is a significantly greater capacity for experiments that require large amounts of clonal cells with identical genetic background, such as needed in state-of-the-art high-throughput screening, especially in proteomics or metabolomics. the availability of clonal cells allows the possibility for comparative perturbation experiments aiming to compare phenotypic outputs derived from a set of single node perturbations. in addition, using cell lines avoids the ethical concerns arising out of human primary neuronal cell culture. the human neuroblastoma cell line sh-sy5y first described in, is a commonly used cell line in studies related to neuroblastoma and neurodegenerative disease. the cell line is a sub-clone of the parent cell line sk-n-sh which was originally established from a bone marrow biopsy of a neuroblastoma patient. the three human diseases most frequently mentioned in literature for sh-sy5y were neuroblastoma, alzheimer’s disease, and parkinson’s disease. a complete genomic characterization of sh-sy5y would thus elucidate the applicability and possible limits for modelling these disease-specific processes in the context of this cell line. undifferentiated sh-sy5y cells have been extensively used as an in vitro model for research in neuroscience. the cell line shows biochemical properties of immature catecholaminergic neurons. studies have found that undifferentiated sh-sy5y express only immature neuronal markers and lack mature neuronal markers. thus, undifferentiated sh-sy5y cell line might not represent an appropriate model for diseases such as pd, which primarily affect differentiated dopaminergic neurons. further, xie et al. reviewed articles on sh-sy5y cell as in vitro model for pd research and found that differentiation of sh-sy5y under a certain treatment results in a more dopaminergic neuronal phenotype, which could be extremely useful for modelling selective dopaminergic cell death in pd. however, they point out that an optimally differentiated sh-sy5y dopaminergic cell model requires further research. a systems genomics analysis could thoroughly assess the genetic mutations in all neuronal markers and therefore place limits on how closely differentiation can model a dopaminergic cell. furthermore, treatment with differentiation-inducing agents enable sh-sy5y cells to become morphologically similar to mature primary neurons, and the different treatment agents result in a variety of neuronal phenotypes. genomic mutations may, however, impose limitations on the possible phenotypes to which sh-sy5y can differentiate. therefore, a complete genomic characterisation of the sh-sy5y cell line can inform the limitations on the possible phenotypes imposed by the genome. systems genomics aims to integrate genomic variation, copy number and structural variation with high-throughput gene expression, metabolomics and proteomic data to explore the genetic architecture of complex traits and multi-factorial diseases. for characterizing the sh-sy5y cell line we integrated information from whole genome sequencing, transcriptomics, and proteomics experiments. firstly, whole-genome sequencing is used to determine the genetic background of the cell line. secondly, transcriptomics and proteomics data is used to study correlation across biological levels. finally, we integrate genomic variants using a network analysis approach to evaluate the suitability of the sh-sy5y cell line as an in vitro model to study various neurodegenerative diseases, including pd.figure systems genomics approach to assess sh-sy5y as a disease model.  whole genome sequencing the whole genome of sh-sy5y cell line was re-sequenced using two different sequencing platforms, those from complete genomics and illumina hiseq to generate a high-quality list of variants. an earlier study has shown the utility of combining information from both platforms and although they show high concordance, each platform alone failed to detect a significant number of exonic variants spread over genes. for the sh-sy5y cell line, dna sequencing by cg and il produced a genome-wide coverage of 57× and 49× respectively. more than 90% of the genome was covered by a minimum read depth of. further, more than 95% of the exome was called with minimum required confidence set by each of the proprietary genotype-calling procedures of the two respective sequencing platforms. the two genome-sequencing platforms – cg and il – produced a total union of, single nucleotide variations with 84% concordance between both platforms. filtering them based on criteria defined in reumers et al. resulted in, snvs with 99% concordance between both platforms. out of these high quality variants 96% were previously reported by dbsnp build and % by genomes project. snvs and small indels were earlier reported by the catalogue of somatic mutations in cancer version the total number of snvs and the percentage of snvs found in genomes project are similar to controls in another study taken from the human genome diversity project. among the high confidence set of mutations in sh-sy5y, were confirmed as somatic in the cosmic database and were also rare in the general population.table high-confidence mutations in sh-sy5y that were rare in the population and also confirmed as somatic in the cosmic database rare mutations are snvs and small indels that were found in less than 5% of the samples in genomes project, exome sequencing project and complete genomics baseline genomes. in the genes column, when distances are given, the mutations are found in intergenic regions and the first gene precedes the mutation whereas the second gene succeeds the mutation. the number of private protein-altering snvs and indels, i.e. protein-altering variants not found in genome project or the exome sequencing project were and respectively, whereas the corresponding numbers for controls were and. the set of control genomes refers to genomes used in another study taken from the human genome diversity project. the reason for the discrepancy between the number of ppa snvs between sh-sy5y and the control genomes is mainly due to our strict filtering strategy as well as the different sequencing platforms that were used. the overlap between ppa variants and the variants in the catalogue of somatic mutations in cancer was and for sh-sy5y and control genomes respectively. the overlap between genes containing ppa snvs and indels and the sanger cancer gene census was similar to controls. we also performed a gene ontology enrichment analysis on the genes containing filtered ppa snvs and indels and also on the genes using unfiltered ppa snvs and indels. however, no go term was significantly enriched. we judged the overall quality of the variant calls and also compared them between the two genome sequencing platforms - cg and il. then, we integrated variant calls from both these platforms by using the filtering criteria provided by an earlier study that compared these two platforms. validation of snvs and small indels validation was performed using the illumina omni- quad microarray, which assays loci from hapmap phase 1– and genomes project. out of heterogeneous snvs that were queried by genotyping, % were concordant, which increased to % for the filtered variants. functional prediction of snvs and small indels the functional effect of snvs and small indels were predicted using annovar by annotation with labels for genomic regions and coding effects. the annotation was performed by taking a consensus across four databases based on choosing the most damaging effect predicted in order to aim for sensitivity. in general, the association of snvs to genomic regions for platform-specific and concordant snvs did not show any significant differences. the association of snvs with exonic, intronic and intergenic regions was 1%, 26–36% and 46–50% respectively, similar to the results in lam et al.. the functional prediction found genes with rare non-synonymous snvs, indels or substitutions. neuroblastoma-relevant genes the genome sequencing of sh-sy5y found genes with rare non-synonymous snvs, indels or substitutions that overlapped with the list of genes containing somatic mutations in the complete genome sequence of untreated primary neuroblastoma tumours. using the hyper-geometric test, genes with rare non-synonymous snvs and indels were found to be significantly enriched among the genes with somatic mutations in primary neuroblastoma tumours. out of these genes, only mutations in genes were predicted as damaging by sorting intolerant from tolerant. furthermore, only genes with rare non-synonymous snvs, indels and substitutions overlapped with the genes associated with neuroblastoma in the literature and therefore were not significantly enriched. however, genes with non-synonymous snvs, indels and substitutions, copy number variations, and structural variations were significantly enriched among the genes extracted for neuroblastoma from literature with an overlap of genes. we also compared genomic mutations in sh-sy5y to other genomic sequencing studies of primary neuroblastoma from patients. among the seven genes mutated at a significant frequency in a study of matched tumour -normal samples, none of the genes carried a rare amino-acid changing mutation. among the rare germ-line variants predisposing to neuroblastoma identified by the target study, palb contained two rare non-synonymous mutations in sh-sy5y. among the coding somatic mutations found in an aggregate of matched tumour/normal samples, we found overlapping variants in sh-sy5y inside genes among which were rare amino-acid changing mutations inside genes. the pediatric cancer genome project conducted a study of patients with metastatic neuroblastoma and found mutations in atrx and alk in 22% and 14% of the patients correspondingly. however, sh-sy5y contained no rare amino-acid-changing mutations in these genes. structural variations and mobile element insertions structural variation events consist of large deletions, duplications, inversions, translocations and complex variations - other combinations of such chromosomal rearrangements. the majority of sv events in the sh-sy5y cell line found by cg were deletions. around 2–3% of cancers show chromothripsis, where tens to hundreds of genomic rearrangements occur in a cellular crisis event. in a certain region, a high density of genomic rearrangements combined with frequent oscillations between two copy number states and frequent occurrences of runs of homozygosity constitute the hallmark of chromothripsis. cytogenetic methods confirm a high density of genomic breakpoints in single cells and indicate that these breakpoints are not a result of parallel rearrangements in different sub-clones. however, for sh-sy5y, dna sequencing found little to no evidence of chromothripsis. most notably, no oscillations between copy number states of 100-kb binned read depth were observed except for some on chromosomes and x. in all of these four chromosomes, the regions with oscillating copy number were near the centromeres or the telomere. even though all of these regions have at least six junctions, these regions consist of highly repetitive dna and are therefore prone to mapping errors. we filtered the svs based on their frequency in the cg baseline genome set and extracted high-confidence calls with a frequency of less than 10% in the baseline. the cg baseline set comprises of genomes of healthy, disease-free individuals, which could be used to filter technical artifacts and variations common in the population. after filtering, the high-confidence deletions overlapped with genes including the following cancer-related genes - pten, which is a key tumour suppressor gene, ctnna, which is a cell contact inhibition gene whose mutation can promote cancer development and formation, mcc and mtus. mobile element insertions refer to insertion of sequences that can change their position within the genome. mei events were detected in the cg data and the majority of those were alu events and l retro-transposons. copy number variation cg identifies discrete coverage levels corresponding to ploidy levels using the distribution of observed normalized coverage values and then uses a hidden markov model to assign copy number levels. since cg essentially uses the relative coverage as a proxy for copy number, it cannot directly infer the absolute copy number levels. however, m-fish experiments have demonstrated that diploidy dominates the copy number level in the sh-sy5y genome. therefore, we calculated the absolute copy number levels assuming that the most commonly found copy number level is we compared the copy number levels found by cg with an earlier measurement using comparative genomic hybridization arrays. the cgh arrays used had a lower limit one million base pairs to their resolution, a drawback not found in second-generation sequencing technologies such as cg and il. to achieve even higher resolutions and reduce false positive and false negative copy number determinations, we normalized the coverage levels to those observed in a reference set of in-house genomes also sequenced by cg. a visual illustration of the copy number levels confirms key features of the cell line described previously such as partial trisomy of chromosome gain of chromosome 2p, 17q, and loss of 1p, 14q and 22q. among the large-scale copy number features identified in the pcgp project, the gain of 17q was common with sh-sy5y. the remaining cnvs that were reported with a lower frequency in pcgp - loss of 1p or loss of 11q were not found in sh-sy5y.figure cnvs detected by whole genome sequencing and array-based cgh by do et al.. the results from whole genome sequencing were from complete genomics and are given in the left half of the chromosomes. the results from the array-based cgh are on the right half of the chromosomes. regions are highlighted for copy number gain and loss. the major events partial trisomy of chromosome and complete trisomy of chromosome gain in 17q and loss in 22q were confirmed.. interestingly, the strongest clinical marker of neuroblastoma that also determines the aggressiveness of neuroblastoma - gain of the gene mycn – actually had a relative coverage of only however, more recently, confirmed the finding that the mycn gene is not strongly amplified in sh-sy5y. also, among the patients sequenced in the pcgp project, only 23% carried an amplification of mycn, which was defined as > copies detect by quantitative-pcr. transcriptome analysis to characterize the sh-sy5y cell line gene expression, we sequenced poly-adenylated rna from undifferentiated sh-sy5y cells using a stranded rna-seq protocol. we used the grch human genome version together with the ensembl release gene annotation. we generated about million paired-end reads of lengths 50 nt. the raw reads are available in the european nucleotide archive database under the study accession prjeb in total, million paired-end reads and million singleton reads could be mapped after quality control filtering and trimming to either the genome or the transcriptome annotation. this yields in 95% mappable reads to either genome or cdna sequences. for rna-seq, we called snvs and small indels using both samtools and the genome analysis toolkit. we retained the intersection of the two outputs and those caller-specific variants supported with a read depth greater than or equal to as a recent study has shown low concordance between multiple variant-calling pipelines, we decided to increase the level of confidence by using two variant-calling pipelines. consequently for rna-seq, out of snvs and small indels detected, % were concordant with unfiltered variants from at least one dna sequencing platform, either cg or il. the rna-seq also showed high sensitivity as % of the filtered exonic variants found through dna sequencing in genes with fpkm >  were also detected by rna-seq. fpkm refers to fragments per kilobase of exon per million fragments mapped, which is a measure of gene expression. this sensitivity of rna-seq for variant detection was higher than that achieved by a recent approach for snv calling in rna-seq data but the absolute number of detected exonic variants was slightly lower than theirs, which might be due to our strict quality filtering of the dna sequencing data. as for neuroblastoma-relevant genes from the neuroblastoma primary tumours, out of the rare, amino-acid changing mutations, of them were detected by rna-seq with the same zygosity. the remaining were likely not detected because of low gene expression as all mutations detected by rna-seq were in genes with fpkm >  and all mutations not found by rna-seq were in genes with fpkm < . genetic copy number vs. gene expression we compared genes that were always expressed in sh-sy5y across different conditions from the geo database to those that were never expressed, whenever the corresponding probe was present. here, the threshold used to determine whether genes were always expressed in the geo dataset was that the expression values of the corresponding transcripts were always greater than the median for all the transcripts in each microarray experiment. similarly, the threshold for genes never expressed in the geo dataset was that the expression values were always less than the median for all the transcripts in each microarray experiment. the frequency of mutations per kilo-base-pair region in genes that were never expressed was, which decreased to for genes, which were always expressed. similarly, the average copy number for genes always expressed was, which was lowered to for genes never expressed. the distribution of expression of genes with copy number greater than was significantly greater than expression of genes with copy number lower than or equal to. further, we compared the genetic copy number and the rna-seq data of the same sh-sy5y sample. the copy number and the logarithm of fragments per kilo-base of transcripts per million mapped reads showed a positive correlation. proteome analysis in order to confirm and correlate the biological expression of the identified genes, the analysis of genomic data was integrated with the analysis of proteomics data analysis of sh-sy5y cell lysates. proteins isolated from whole sh-sy5y cell lysates were fractionated by sds-page, in-gel digested using trypsin, and the recovered peptides analysed by lc-esi ms/ms. the spectra were analysed by maxquant against a combinatorial human protein database. this database contained the human genome reference hg protein entries combined with the sequence variants found from the dna sequencing done here. specifically, the coding sequences of ensembl genes were modified using all the homozygous exonic variants. in total, out of acquired ms/ms spectra, were matched to a peptide that could be mapped to a total of protein identifiers. from these, protein groups corresponding to proteins were identified in the human genome reference hg protein annotation. further, proteins were identified using sequence variants found from the dna sequencing done here. see additional file for the abundances of proteins detected. to increase peptide spectrum matches, we extended the human genome reference hg protein database by inserting reference sequences containing the homozygous snvs and indels in the exonic regions identified in this study. the peptides that were additionally mapped to the extended protein database were then used for validation of genetic variants. this approach validated snvs, indels and block substitution. gene expression vs. protein expression gene expression levels from rna-seq and the corresponding protein abundance were compared to detect correlation. we plotted the logarithm of fpkm against the logarithm of intensity–based absolute quantification score of genes, which were detected in both rna-seq and proteomics, which equals 93% of proteins detected by proteomics. as expected, the correlation was weakly positive. network-based analysis of suitability of sh-sy5y cell line as in vitromodel in order to evaluate the suitability of sh-sy5y as in vitro model we adopted a framework described in that integrated genetic sequence information, and topology analysis of either disease or process-specific networks. we applied betweenness-centrality ratio, a metric that allowed us to quantitatively assess the impact of genes mutated in the cell line on the disease or process network. bc-ratio is normalized to lie between and. it indicates how much the information flow in the network might be altered by changes in the cell line: the higher bc-ratio, the greater the impact of the mutations, the lower the genetic integrity of the cell line in the context of a specific network.. the bc-ratio metric applied to four neurodegenerative diseases ranked alzheimer’s as highest bc-ratio = , pbc-value =  ), followed by parkinson’s disease, huntington disease, and amyothrophic lateral sclerosis. to investigate the appropriateness of the cell line for perturbation experiments in the context of pd, we also scored the different processes corresponding to the hallmarks of pd. see figure  for a protein-protein interaction network resulting in the cell line suitability scores for neuroblastoma and als and figure  for the network for glycolysis and reactive oxygen species metabolism in pd. additional file contains network visualizations for each disease and each module in pd. the role of the visualization is to give an idea of how damaged genes are distributed in the network, to which genes they are connected and how they could alter processes of interest.table network statistics and cell line scoring of the neurodegenrative diseases and neuroblastoma flow bc stands for flow betweenness centrality. the pbc-value refers to the network randomization test described in the methods section. flow betwenness is defined for connected networks. if the entire network is not connected, flow betweenness is computed for connected components. in the table the flow betweenness centrality of the largest connected component is given and marked by *. nb – neuroblastoma; ad – alzheimer’s disease; pd – parkinson’s disease; hd – huntington’s disease; als – amyotrophic lateral sclerosis. network statistics and cell line scoring of the parkinson’s disease modules flow bc stands for flow betweenness centrality.the pbc-value refers to the network randomization test described in the methods section. flow betwenness is defined for connected networks. if the entire network is not connected, flow betweenness is computed for connected components. in the table the flow betweenness centrality value of the largest connected component is given and marked by *. protein-protein interaction network for two diseases – neuroblastoma and for als. red nodes refer to omim-derived genes mutated in the cell line. orange nodes refer to the omim-derived genes which are intact in the cell line. dark green nodes represent genes coming from the network expansion, which are mutated in the cell line. light green nodes represent genes coming from the network expansion, which are intact in the cell line. nodes are scaled to the magnitude of their betweenness centralities. blue edges show connections between pairs of nodes in which at least one was damaged in the cell line. for neuroblastoma, mutations in the central genes – nme and alk contribute to the high bc-ratio. for the als network, mutations occur in genes with lower betweenness centrality which results in the lower bc-ratio. protein-protein interaction network for two pd map modules - glycolysis and ros metabolism. red nodes refer to genes mutated in the cell line. orange nodes refer to genes which are intact in the cell line. nodes are scaled to the magnitude of their betweenness centralities. blue edges show connections between pairs of nodes in which at least one was damaged in the cell line. for glycolysis, mutations in the central genes lead to the high bc-ratio. mutated genes lie on the periphery of the ros metabolism network resulting in a low bc-ratio. discussion the complete system-wide omics analysis of the undifferentiated sh-sy5y cell line provided many opportunities for studying concordance across biological levels. firstly, protein-altering variants, copy number variations and structural variations show consistency with the expression level of genes in the rna-seq data. for instance, variant calling from rna-seq data showed % sensitivity for the filtered variants detected from dna-seq in genes with an fpkm >  secondly, sequencing both the dna and the rna showed a high degree of concordant snvs and indels and served as further validation of dna sequencing results and our filtering procedure for genomic variants. finally, the knowledge of genomic variants enhanced protein identification in proteomics experiments where it increased the number of proteins detected from to. a comparison of the two whole genome sequencing platforms – complete genomics and illumina – confirmed a significant number of variants that were discordant, either due to regions, which were not called by one of the platforms or due to platform-specific errors. considering that 1% of these platform-specific variants residing in genes were associated with exonic regions, using both platforms would increase the detection of potentially functionally important snvs. also only 75–80% of platform-specific variants were found in the known list of variants, indicating that a significant number of mutations are somatic mutations or sequencing errors that need further investigation. further, the larger platform-specific biases for small indels and structural variations also strengthened the argument for using both platforms for greater sequencing accuracy. neuroblastoma-relevant genetic variations without the genome sequence of the healthy cells from the patient, the difficulty of identifying somatic mutations in sh-sy5y precludes any direct comparison to the somatic mutations in primary neuroblastoma tumours. thus, we compared rare mutations in sh-sy5y with the somatic mutations found in primary neuroblastoma tumours. the union of genes from different patient samples resulted in a list of genes, which we considered as relevant to neuroblastoma. the result that only genes with rare amino-acid changing mutations predicted as damaging in sh-sy5y overlapped with the list of genes with somatic amino-acid changing mutations in primary neuroblastoma was expected as there was considerable heterogeneity among the different tumour samples themselves. specifically, the sequencing of neuroblastoma primary tumours found few recurrent genes with amino-acid changing mutations - only genes contained amino-acid changing mutations in more than tumour sample and only three genes in more than two tumour samples. in this context, an overlap of genes with rare amino-acid changing mutations suggests that sh-sy5y includes elements of the genetic architecture of a variety of neuroblastoma tumours. therefore, it might offer opportunities to investigate traits arising from amino-acid changing mutations found in different kinds of neuroblastoma tumours. as for structural variations, recurrently occurring structural alterations in the neuroblastoma primary tumours in ptprd, odz and csmd were not found among high-confidence rare structural alterations in sh-sy5y. finally, chromothripsis, which was earlier identified in 18% of primary neuroblastoma, was not detected in sh-sy5y. therefore, sh-sy5y would not serve as an appropriate model for neuroblastoma tumours suffering from chromothripsis. suitability of sh-sy5y as an in vitromodel for neurodegenerative diseases the bc-ratio relies on the network betweennness centrality to quantify genetic changes in the cell line. it is important to consider that the bc-ratio does not aim to interpret the functional effects of these mutations. indeed, mutations could cause a loss, gain, change or no change in the function. however, such interpretations require assumptions. in the case of neuroblastoma, we know that the sh-sy5y cell line was derived from a neuroblastoma patient. hence, it can be assumed that many genetic mutations, which were observed in the neuroblastoma network, are indeed typical neuroblastoma mutations. specifically, copy number variations in the central neuroblastoma genes nme and alk, both of which are found to be mutated in neuroblastoma as well support this reasoning. briefly, it can be assumed that sh-sy5y cells are a good model for neuroblastoma cells and a high bc-ratio was expected. however, this case where a neuroblastoma cell line is expected to model neuroblastoma is a special case. this interpretation is similar to the approach of domcke et al. of choosing the ‘most suitable models’ for ovarian tumour among ovarian cancer cell lines based on genomic and mrna expression profiles. typically, in contrast, experimentalists use cell culture models to study processes lacking such ideal identity matching. routinely, perturbation experiments are used to cause controllable genetic defects in networks, which are overall assumed to be intact. however, potential genetic changes, which are ignored, can introduce errors into the interpretation of results. the main scope of the bc-ratio is to provide distances between networks of interest in a cell line and their ideal reference networks. a bc-ratio of indicates an unchanged network representing an ideal reference for perturbation experiments. the pbc value evaluates if a bc-ratio >  was due to chance or it resulted from strong correlation between the nodes with high betweenness centrality and genes mutated in the cell line. note that a non-significant pbc-value does not exclude any functional effects of network changes and it is strongly recommended to collect additional information for evaluating potential functional effects. while our aim was not to predict functional effects related to rna or protein abundance, we are convinced that the provided data will support case-specific interpretations of functional network integrity. a central question of our investigation was whether context -specific networks are overall changed in sh-sy5y. the absence of network changes was considered as ideal case for controllable perturbation experiments. in the following paragraphs, we first discuss the integrity of disease-specific networks in sh-sy5y and secondly, focus on pd in more detail. the cell line suitability metric applied to four neurodegenerative diseases ranked them from those with the largest to the smallest changes in the network. alzheimer’s disease was ranked as the highest, followed by parkinson’s disease, huntington’s disease and, finally, amyotrophic lateral sclerosis. bc-ratios in combination with the pbc-values tell us that genes mutated in the cell line are more likely to be involved in the information flow control in ad network than the same type of genes is in the networks of other diseases. overall, except for the neuroblastoma, bc-ratios were rather low, and pbc-values were not significant for all but neuroblastoma and alzheimer’s disease. these results indicate that the networks related to pd, hd and als are largely intact in the sh-sy5y cell line and the latter can be considered a suitable model for perturbation experiments targeting these diseases, while ad might require a special caution. it is also important to bear in mind that none of these disease-specific networks were perfectly intact. indeed, the bc-ratio for each network was computed as greater than, which indicates that changes in these networks were detected. a prediction of potential functional effects of these network changes, including a classification in gain or loss of function goes beyond the scope of the bc-ratio but the provided data on rna and protein abundances can support such interpretations. genes from which the networks for the neurodegenerative diseases and neuroblastoma have been grown, were derived from the omim databse where they have been classified as having disease causative effect. we also calculated the bc-ratios of the respective networks assuming that only these seed genes have had been damaged in the cell line. in all cases the number of mutated genes decreased while the bc-ratios of the cell lines increased for all but neuroblastoma. note that in the case of neuroblastoma, the main actors which contributed to the high bc-ratio computed on the real data, already were omim genes, damaged in the cell line. therefore a lower bc-ratio for neuroblastoma was expected in the absence of any other damaged genes. the bc-ratio of the huntington’s disease increased considerably. this is mostly due to the gene htt followed by prnp, which have the highest betweenness centrality in the hd network. alzheimer’s disease scores next with app gene making the main contribution to the magnitude of the bc-ratio, being the first on the list of all genes. it is immediately followed by more omim-derived genes while their betweenness centrality is more than two times smaller. similar effect can be observed in the case of amyotrophic lateral sclerosis: sod has the highest betweenness centrality, followed by vcp and dctn with almost twice as much lower betweenness. in this experiment, parkinson’s disease had the lowest bc-ratio although the list was dominated by omim-derived genes. the explanation might lie in the low betweenness centralization of the network. network centralization reflects how much variation is there in the centrality scores among the nodes. the value of indicates that the scores are rather equally distributed among the nodes and there are no nodes with apparent brokering role. therefore changing the node labels from intact to damaged in such network produces smaller effect than in the remaining diseases.table bc-ratios of the neurodegenerative diseases and neurobalstoma computed under assumption, that only omim-derived seed genes have been mutated the increase in the bc-ratios of all neurodegenerative diseases with the causative omim genes being considered as damaged, suggests that the genes which indeed are mutated in the cell line have smaller positional advantages than omim genes in the corresponding networks and therefore might less be involved in the information flow control. we believe that this observation supports our conclusion about the overall sh-sy5y cell line suitability for the experimental studies of the neurodegenerative diseases. to further analyse the suitability of the cell line as an in vitro model for pd, we studied the integrity of disease related sub-networks including mitochondrial dysfunction, reactive oxygen species accumulation and calcium homeostasis. fujita et al. recently provided a comprehensive pd map, which integrates metabolic reactions, gene regulation and signalling processes in this complex network. the network consists of sub-networks, which correspond to pathways and processes. in these networks, nodes represent genes, proteins or small molecules, and edges represent molecular interactions. we selected seven modules corresponding to hallmarks of pd from the pd map and scored them according to the bc-ratio metric. the idea was to segregate pd-related modules as more suitable for studying gain or loss of function of a module. in descending order of bc-ratio, the modules were glycolysis, mitochondria, calcium signalling, apoptosis, dopamine metabolism, ubiquitin proteasome system and reactive oxygen species metabolism. the pbc values indicated no significant changes in these networks. however, we found that the glycolysis module is the most impacted, which suggests that information flow through this module could be significantly altered by the genetic background of sh-sy5y. one explanation for glycolysis being highly impacted could come from the fact that sh-sy5y is a cancer cell line. it has been shown that glycolytic rates are higher in cancer cells than healthy ones, known as the warburg effect, which has been suggested to confer a proliferative advantage to tumour cells.conversely, the ros metabolism module was the least impacted, which would indicate that the genes, mutated in the cell line do not compromise the information flow in this network. for the functional interpretation beyond the bc-ratio, it is important to consider that most of the genetic variations found in the pd modules do not carry mutations, but almost exclusively copy number gains. indeed, this gain of copy numbers is expected because of the partial trisomy of chromosomes and and the full trisomy of chromosome however, variations that alter the protein sequence occur in very few genes. among all the damaged genes found in the seven modules from the pd-map, only four contained non-synonymous snvs. the sorting intolerant from tolerant tool predicted only as damaging. the gene myo found in the calcium signalling module and involved in intracellular vesicle transport, contains a heterozygous non-synonymous mutation with a sift score of. it also plays an important role in trafficking and activity-dependent recruitment of ampa receptors to synapses. the gene hadh, found in the mitochondria module and which catalyses several reactions in beta-oxidation, also contains a heterozygous non-synonymous mutation with a sift score of the gene gzmb, found in the apoptosis module and which codes for a serine protease that is used by activated cytotoxic t lymphocytes to induce cell apoptosis, contains a rare non-synonymous snv with a sift score of. however, in each of these three genes, only one copy of the gene contains such rare non-synonymous mutations. further, four genes in the mitochondria module and one gene in the apoptosis module had a copy number of one, leaving only a single working copy. ndufa has been shown to have a high degree of nitration and be associated with the oxidative damage to mitochondrial complex i. tomm plays an important role in mitochondrial clearance controlled by pink –park pathways. atp5l is involved in hydrogen ion trans-membrane transport activity. xrcc codes for the ku protein, which repairs double bond breakage of dna through the non-homologous end-joining pathway. therefore, one should consider the damage to these genes when designing loss of function perturbation experiments. indeed, it is important to test if the mutations per se cause a loss of function and to keep in mind that a function which is already lost cannot be lost a second time. except for the genes discussed above, all the genes in the modules of the pd map were free of translated mutations in sh-sy5y. genes related to dopamine metabolism in sh-sy5y in summary, all the genes involved in the dopamine metabolism pathway except cyp2d adh1b and ugt1a were expressed with or without specific treatments. further, cyp2d which contributes to dopamine biosynthesis through an alternative cytochrome p450-mediated pathway shown to exist in rats, has a copy number of one, which may contribute to the lower expression. among the genes in the dopamine metabolism pathway that were expressed with or without specific treatments, only adh1c and dbh contained a rare protein-altering mutation and ddc and akr1b had a copy number of three. only the mutation in adh1c was predicted as damaging by sift. therefore, the large majority of genes involved in dopamine metabolism contains no genetic damage and can be expressed under certain treatments.figure dopamine biosynthesis and degradation and the genes coding key enzymes. this metabolic pathway is a modified version of a figure from meiser et al.. in addition to the enzymes involved, the genes coding for those enzymes have been added using the humancyc. the mutations in sh-sy5y affecting enzyme-encoding genes have been annotated only if they were rare protein-altering snvs or indels, cnvs, or svs.  in summary, we provide here the first whole genome analysis of the sh-sy5y cell line, which is widely used as a model for various neurological diseases. we characterise the different types of biological information of sh-sy5y – genomics, transcriptomics and proteomics - and compare the relationships between them in terms of expression levels and variants. the data related to the genome sequencing, rna sequencing and proteomics can be found at the following url - http://systemsbiology.uni.lu/shsy5y. additionally, we also compare the sh-sy5y genome using two widely used whole genome sequencing platforms and show that using the two platform-specific coverages and sensitivities to different indel sizes supplements each other. finally, our analysis based on context-specific network integrity ranked the integrity of sh-sy5y for parkinson’s disease as higher than for alzheimer’s disease but lower than for amyotrophic lateral sclerosis and huntington’s disease. in the context of different pd related sub-systems, the same approach ranked ros metabolism as the most intact followed by the ubiquitin proteasome system, dopamine metabolism, calcium signalling, mitochondria and glycolysis. overall, most of the genes belonging to the major parkinson’s disease pathways and modules were intact in the sh-sy5y genome. particularly, each analysed gene related to pd has at least one intact copy in sh-sy5y. therefore, somatic mutations do not significantly alter pd-related pathways in sh-sy5y. importantly, the systems genomics analysis was performed on the undifferentiated cell line although differentiated cell lines are considered better suited as in vitro models for pd. our novel bc-ratio method for cell line scoring integrates a genomic characterisation of the cell line, core set of disease or process-related genes and a protein-protein interaction network. this method is not restricted to a specific cell line or disease and can be broadly applied. furthermore, if multiple sequenced cell lines are available, the bc-ratio scores can also be used as a guide to select the cell line best suited for studying a particular disease. 
most of the genes belonging to the major parkinson’s disease pathways and modules were intact in the sh-sy5y genome. the disease-specific network analysis approach ranked the genetic integrity of sh-sy5y as higher for pd than for alzheimer’s disease but lower than for huntington’s disease and amyotrophic lateral sclerosis for loss of function perturbation experiments. the human neuroblastoma cell line, sh-sy5y, is a commonly used cell line in studies related to neurotoxicity, oxidative stress, and neurodegenerative diseases. although this cell line is often used as a cellular model for parkinson’s disease, the relevance of this cellular model in the context of parkinson’s disease and other neurodegenerative diseases has not yet been systematically evaluated. we have used a systems genomics approach to characterize the sh-sy5y cell line using whole-genome sequencing to determine the genetic content of the cell line and used transcriptomics and proteomics data to determine molecular correlations. further, we integrated genomic variants using a network analysis approach to evaluate the suitability of the sh-sy5y cell line for perturbation experiments in the context of neurodegenerative diseases, including pd. sh-sy5ycell linewhole genome sequencingrna-seqproteomicscell line suitability evaluationparkinson’s diseaseissue-copyright-statement© the author 2014 specifically, each analysed gene related to pd has at least one intact copy in sh-sy5y.
 the human neuroblastoma cell line, sh-sy5y, is a commonly used cell line in studies related to neurotoxicity, oxidative stress, and neurodegenerative diseases. although this cell line is often used as a cellular model for parkinson’s disease, the relevance of this cellular model in the context of parkinson’s disease and other neurodegenerative diseases has not yet been systematically evaluated.  we have used a systems genomics approach to characterize the sh-sy5y cell line using whole-genome sequencing to determine the genetic content of the cell line and used transcriptomics and proteomics data to determine molecular correlations. further, we integrated genomic variants using a network analysis approach to evaluate the suitability of the sh-sy5y cell line for perturbation experiments in the context of neurodegenerative diseases, including pd.  the systems genomics approach showed consistency across different biological levels. most of the genes belonging to the major parkinson’s disease pathways and modules were intact in the sh-sy5y genome. specifically, each analysed gene related to pd has at least one intact copy in sh-sy5y. the disease-specific network analysis approach ranked the genetic integrity of sh-sy5y as higher for pd than for alzheimer’s disease but lower than for huntington’s disease and amyotrophic lateral sclerosis for loss of function perturbation experiments. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords sh-sy5ycell linewhole genome sequencingrna-seqproteomicscell line suitability evaluationparkinson’s diseaseissue-copyright-statement© the author 2014 