Whole-exome sequencing data of suicide victims who had suffered from major depressive disorder

Suicide is one of the leading causes of mortality worldwide; it causes the death of more than one million patients each year. Suicide is a complex, multifactorial phenotype with environmental and genetic factors contributing to the risk of the forthcoming suicide. These factors first generally lead to mental disorders, such as depression, schizophrenia and bipolar disorder, which then become the direct cause of suicide. Here we present a high quality dataset (including processed BAM and VCF files) gained from the high-throughput whole-exome Illumina sequencing of 23 suicide victims – all of whom had suffered from major depressive disorder - and 21 control patients to a depth of at least 40-fold coverage in both cohorts. We identified ~130,000 variants per sample and altogether 442,270 unique variants in the cohort of 44 samples. To our best knowledge, this is the first whole-exome sequencing dataset from suicide victims. We expect that this dataset provides useful information for genomic studies of suicide and depression, and also for the analysis of the Hungarian population.


Background & Summary
Suicide is the 10th leading cause of mortality in the world 1,2 . It is a complex, multifactorial behavioural condition in which proximal and distant risk factors, such as mental disorders, genetic and epigenetic factors, social, cultural and economic environment, gender, age, personality, eating disorders, as well as drug or alcohol abuse 3-7 play important roles. The presence of early child abuse has also been described as a contributing factor to suicidal behaviour 8 . Studies have repeatedly shown that suicide is most common among patients with mood disorders, such as major depressive disorder and bipolar disorder 9 .
Whole-exome sequencing (WES) is a cost-effective and powerful tool for the analysis of complex and rare genetic diseases 10 . WES technique allows a base-pair comparison of exomes and consequently the examination of rare genetic variants, which may play a role in suicide. Alterations in copy number variations (CNVs) and short-tandem repeats (STRs), as well as the viral infections can be also analysed by WES.
We performed high-depth, paired-end, short-read WES on 23 suicide victims suffered from major depressive disorder (MDD; this information is available in EGA sample metadata) and 21 control patients who died from other diseases (EGA sample metadata). Harvesting of brain tissues was approved by the local ethics committee 11 . Informed consents from the next kin have not been obtained from the victims because the local regulations did not require consent for autopsy. All data were handled anonymously. Both of the cohorts were of Hungarian ethnicity. The Hungarian population is known to have had a very high annual suicide rate in the previous century 12,13 . The DNA samples derived from post-mortem brains were sequenced using the Agilent SureSelectXT Human All Exon V5 + UTRs kit and Illumina HiSeq 2000 sequencing platform resulting in an average of 117 million paired-end reads.
We obtained very high coverage per base position in both cohorts; in the suicide cohort 95.4-97.5% (quartiles) of target regions had higher than 20-fold and 82.8-90.5% of target regions had higher than 40fold coverage. In the Hungarian control cohort these values were 96.3-97.5% for 20-fold and 81.6-91.6% for 40-fold coverage, respectively. We called on average~130,000 variants per sample. Altogether we observed 442,270 unique variants in the control and suicide cohorts. Furthermore, we determined a total of 320,976 single nucleotide variants (SNVs) and 47,313 insertions/deletions (indels).
The library preparation kit used in this study, as well as the obtained very high coverage of sequencing reads that we obtained enabled us to capture information not only about the exons but also about the 5′and 3′-UTRs, the promoters to a certain length, and even about off-target sequences, such as introns, intergenic regions and infecting viruses. Our dataset had 286,754 high coverage regions (exceeding 20fold coverage) in all samples that can be suitable for downstream analysis of CNVs.
Here we describe the sample collection methods, the library preparation and sequencing method, the currently available data records, and technical validations for our dataset. A schematic overview of this study, including the experimental procedure and the bioinformatics workflow, is also presented (Fig. 1).
To our knowledge, this dataset is the first individual-level, open-access WES data release targeted to suicide. The dataset described here has been thoroughly analysed in our related manuscript 6 , where it has been shown that the TGF-β signalling pathway and three voltage-gated calcium ion channels may play an important role in the pathogenesis of suicidal behaviour. Our analysis proposed that suicide may be an umbrella-disease like disorder.
Here we provide the processed BAM files, as well as the processed Variant Call format (VCF) files for each of the samples containing the analysis ready variants of GATK HaplotypeCaller pipeline that allows for the testing of different variant filtration strategies that we had used in our original manuscript or for hypothesis testing of candidate genes (covered by SureSelect V5 All Exon plus UTR kit) without tedious bioinformatic processing.
Our dataset provided here in this manuscript can be re-used to investigate suicide and MDD. These data can be a valuable resource for investigating genetic variants, genes and signalling pathways to identify novel factors related to these disorders, and may provide novel information for the investigation of the Hungarian population or in general for studies of genetic polymorphisms of human population.

Methods
The methods described here are an expanded version of those described in our recently published, related publication 6 .

Medical history of the suicide victims
This study is based on whole-genome sequence data from occipital, cerebellar and somatomotor cortex regions of the brains of 23 suicide victims (15 males and 8 females) and 21 control patients (14 males and 7 females) who died suddenly from other causes (this information is available in EGA sample metadata). Samples were obtained by autopsy at the Department of Forensic Medicine of the Semmelweis University Medical School. For all of the suicide victims, a psychiatric diagnosis of major depression was on record. The diagnoses had been made by experienced psychiatrists on the basis of Diagnostic and Statistical Manual of Mental Disorders Fourth Edition (DSM-IV). Patients with a history of schizophrenia, epilepsy and other disorders were excluded from the study. The medical, psychiatric, and drug history of the suicide victims was obtained by psychological autopsy, which included interviews with the attending physician and with family members, as well as data obtained from medical records. The participants included in this study had no history of drug or alcohol abuse, and had not used antidepressant www.nature.com/sdata/ SCIENTIFIC DATA | 6:190010 | https://doi.org/10.1038/sdata.2019.10 medication for at least two months prior to death. This information was confirmed by toxicology tests on blood samples. Suicide victims died by hanging (n = 16), drug overdose (n = 6), or jumping from height (n = 1). Causes of death in control subjects are presented in EGA sample metadata. Harvesting of tissues was approved by the local ethics committee.

Tissue collection, dissection, and storage
Brains were obtained 1 to 10 h after death (this information is available in EGA sample metadata). After removal from the skull, the brains were cut in six major pieces (four cortical lobes, basal gangliadiencephalon, and lower brain stemcerebellum), rapidly frozen on dry ice, and stored at −70°C until dissection (2 days to 2 months). At the time of dissection, the brain samples were sliced into 1 mm-thick coronal sections at a temperature of 0 to 10°C. The cortical areas were cut out of the sections by a fine microdissecting (Graefe's) knife. In all cases, cortical samples were always taken from the right hemisphere. The samples were stored in airtight containers or plastic tubes at −70°C until further use.

DNA purification
Genomic DNA (gDNA) samples were isolated from the occipital, cerebellar or somatomotor cortex regions (Table 1 (available online only)) using the DNeasy Blood and Tissue Kit (Qiagen) according to the manufacturer's recommendations for the spin-column protocol, using 30 mg staring tissue material from each sample. In short, tissue samples were cut into small pieces, and then lysis buffer (provided by the kit) and proteinase K were added. Lysis reactions were carried out at 56°C until complete lysis was obtained. DNeasy Mini spin columns (kit's component) were used for the isolation of gDNA from the lysate. Elution was carried out twice to a final volume of 100 μl per elution.

Quality and quantity check of DNA
The quality of the DNA was checked prior to downstream analysis by NanoDrop (Thermo Fischer Scientific) measurement. High-quality samples with an OD 260/280 ratio ranging from 1.8 to 2.0 were subjected to the next step. Qubit 2.0 fluorometer (Invitrogene) was used for genomic DNA quantification before library preparation (Table 1 (available online only)).

Whole-exome sequencing
Paired-end, indexed libraries for Illumina sequencing were prepared from the DNAs of the 23 suicide victims and the 21 controls using post-mortem brain cortex tissues as a source. Whole-exome sequencing was carried out as previously described 14 with slight modifications. In short, 200 ng from each of the qualified gDNA samples in 100 μl volume were sheared into fragments of approximately 200-450 bps in Covaris microTUBEs (Covaris microTUBE AFA Fiber Pre-Slit Snap-Cap 6 × 16 mm) with Covaris s220 high-performance ultrasonicator. The quality of the DNA was checked prior to downstream analysis, using the Agilent Bioanalyser 2100 (Agilent Technologies) and the High Sensitivity DNA Analysis Kit (Agilent Technologies; Fig. 2). The end-repair, A-tailing (3 0 adenylation) and ligation of the paired-end adaptors steps were executed with the Agilent SureSelect Human All Exon V5 + UTRs kit (Agilent Technologies) as was described in the following manual: SureSelect XT Target Enrichment System for Illumina Paired-End Sequencing Library Illumina HiSeq and MiSeq Multiplexed Sequencing Platforms Protocol version 1.6 (SureSelect XT ). Amplification of the adaptor-ligated libraries was performed using the Herculase II Fusion DNA polymerase (Agilent Technologies) following instructions as set forth in the above mentioned manual. Purification after the enzymatic reactions was carried out by using Agencourt AMPure XP beads (Beckman Coulter). Quality and quantity of the libraries were determined by 2100 Bioanalyser instrument and DNA 1000 Assay (Agilent Technologies; Fig. 3). DNA library samples were hybridized with target-specific Capture Library using the SureSelect XT Reagent Kit, HSQ for the HiSeq platform (Agilent Technologies) based on the SureSelect XT protocol. Incubation of the hybridization mixture was carried out at 65°C for 24 h. After hybridization, the targeted molecules were captured on streptavidin-coated magnetic beads. The SureSelect enriched libraries were amplified by PCR with 6-bp indexing primers (Table 2 (available online only)) following the SureSelect XT manual's postcapture PCR recommendations using 12 cycles. The amplified libraries were purified using AMPure XP beads. The quality and quantity of the indexed libraries were checked by Agilent Bioanalyser and High Sensitivity DNA Assay (Fig. 4, Table 3). The quantities were also measured by Qubit 2.0 (Table 3). Two or three indexed libraries were pooled (Tables 2 (available online only

Sequence QC
FastQC (version v0.11.2) was used to generate quality reports (Data Citation 2) using the aligned and base quality recalibrated BAM files.

Variant calling
The variants were called by the GATK HaplotypeCaller (version: 3.3-0-g37228af) algorithm using the GVCF mode cohort analysis according to the Best Practice Guidelines [15][16][17] . Raw Variants were recalibrated by GATK VariantRecalibrator according to the best practice recommendations using the HapMap, Omni, dbSNP, 1000 G, Mills datasets for training included in the GATK GRCh37 bundle resulting analysis ready variants (Data Citation 2).

Data records
The raw FASTQ files, the GRCh37 aligned BAM files, and the variants in VCF 4.1 format associated with the samples analysed in this study are available on request at EGA (Data Citation 1). We also deposited the QC metric files of Picard MarkDuplicate, and QC reports of FastQC of the analysis ready Bam files at Figshare (Data Citation 2) and in Supplementary Table 1.

Technical validation
Quantitation of the purified DNA samples The isolated DNA samples were quantified by Qubit (Life Technologies) fluorometer using Qubit dsDNA HS (High Sensitivity) Assay Kit, which is highly selective for double-stranded DNA over RNA and is designed to be accurate for initial sample concentration ranging from 10 pg/μl to 100 ng/μl. DNA samples were diluted to 4 ng/μl with 1X Low TE Buffer. NanoDrop (Thermo Fischer Scientific) measurements were also performed to assess quantity and quality of DNA, 260:280 and 260:230 ratios greater to 1.8 were accepted.

Quality control of the sheared DNA samples
The quality of the sheared DNA samples (200 ng of each) were checked prior to downstream analysis, using the Agilent Bioanalyser 2100 (Agilent Technologies), and High Sensitivity DNA chip and reagent kit. The electropherogram showed a DNA fragment size peak (for each of the samples) at around 150 bps (Fig. 2).

Quality check of the amplified samples
Agilent Bioanalyser 2100 (Agilent Technologies) and DNA 1000 Assay were used for the quality and quantity control of the libraries after PCR. The sample fragments sizes were between 225 and 275 bps (Fig. 3).

Assessment of the quantity and quality of the indexed library DNAs
The quality of the amplified, indexed libraries were determined before multiplexing using the Agilent Bioanalyser 2100 (Agilent Technologies) and High Sensitivity DNA Assay (Table 3). The peak of the DNA fragment size was between 200-450 bps for each of the samples with an average of 355 bps (Fig. 4, Table 3).

Technical replicates
Two independent technical replicate libraries were prepared and sequenced from the same patients (three suicide victims and two controls; Table 4).

Quality control of raw reads
Illumina BCL files were converted to FASTQ files by the standard Illumina protocol to remove lowquality reads and adaptors.

Quality control of mapped reads
We used Picard tools MarkDuplicate algorithm to identify PCR or optical duplicates. According to MarkDuplicates statistics, we had very low percentages of duplicated reads (on average 0.624%) in our cohort. Genome Analysis Toolkit (GATK version 3.3) BaseRecalibrator base quality score recalibration tool was used to generate the final quality recalibrated BAM files for downstream analysis and variant calling according to best practices [15][16][17] .
To check the quality of our mapped reads we used FastQC to generate quality reports of base quality score recalibrated BAM files. According to the FastQC reports we had no adapter contamination or overrepresented sequences indicating contamination in any of the analysed sequences. Furthermore, the majority of sequences were of high quality (average Phred quality score >34) and even the last two base pairs of the sequences had an average Phred quality score >30. On average 99.9% of reads were mapped, 99% were properly paired, and we had less than 0.1% singletons. The average read counts in the cohort were 117 013 694.9 while the average insert sizes were in the range of 157. .0 base pairs (Table 5 (available online only)). According to the main quality indicators of our dataset all data files fall within acceptable parameters as shown in Table 5 (available online only).

Variant calling
The identification of different genomic variants including SNVs, indels, multiple nucleotide variants (MNVs), Genome Analysis Toolkit (GATK version 3.3) HaplotypeCaller algorithm was used with the GVCF mode cohort analysis according best practice [15][16][17] . HaplotypeCaller was run with the "-stand_emit_conf 10 -stand_call_conf 30" parameters. Variants were recalibrated by GATK VariantRecalibrator according to the best practice recommendations using the HapMap, Omni, dbSNP, 1000 G, Mills datasets for training included in the GATK GRCh37 bundle. We used 100.0, 99.9, 99.0, 90.0 tranches for variant recalibration and the 99.0 truth sensitivity level was used at ApplyRecalibration algorithm to identify QC PASS-ed variants.

Usage Notes
Ethical approval for the study was obtained from the following institutional review boards: The brain samples were collected by the Human Brain Tissue Bank (HBTB; Semmelweis University, Budapest). The activity of the HBTB has been authorized by the Committee of Science and Research Ethic of the Ministry of Health Hungary (ETT TUKEB: 189/KO/02.6008/2002/ETT) and the Semmelweis University Regional Committee of Science and Research Ethic (No. 32/1992/TUKEB) to remove human brain tissue samples, collect, store and use them for research.
The HBTB is a member of the BrainNet Europe II. consortium. The activity of the HBTB is in accordance with the ethic and safety rules and the scientific requirements of the consortium.  Conducting genetic testing on tissue samples and then sending them abroad has been authorized by the Semmelweis University Regional Committee of Science and Research Ethics (No. 34/2002/TUKEB).
The Institutional Review Board (IRB) at Stanford has made the following determination about the activity of the study, titled "Exome sequencing analysis of suicidal behavior using high-throughput DNA sequencing", based on Office for Human Research Protections (OHRP) and Food and Drug Administration (FDA) regulations and guidance: this project does not require submission to the Stanford IRB, because the study does not involve human subjects, it is not about living individuals (NOT-H3 1/1; November 26, 2014).

The following explains how to access the dataset provided in this manuscript
The data is made available through the European Genome-phenome Archive (EGA). Search for the Data Access Committee (DAC) ID of our study (EGAD00001004184) at the homepage of EGA (https://egaarchive.org/). This will show detailed information on our dataset, data providers, DAC and any other related documentation. You can also find information about who to contact about access to this dataset. Access to data will be allowed to qualified researchers for appropriate health related studies.
A DAC group is responsible for the access of the datasets deposited in EGA. Request for data access will be referred directly to our Data Access Committee: https://ega-archive.org/ datasets/EGAD00001004184. Data access decisions can be passed to the EGA in two ways: 1) By emailing helpdesk@ega-archive.org with the email address of each applicant and confirmation of the dataset/s to provide access. The EGA will then create an EGA account with the relevant access permissions. 2) By using the EGA DAC Admin tools (available to DAC's dealing with more than 5 data access applications/month).
If you need to request access to this data set, please contact: DAC for Hungarian human exome team (Department of Medical Biology, University of Szeged, Hungary) Contact person: Dr. Zsolt Boldogkői Email: boldogkoi.zsolt@med.u-szeged.hu Applicants will be asked to complete the Data Access Agreement (DAA) (including a brief summary of the proposal, proposed usage of the dataset, the storage of data, so the DAC can determine if the planned usage falls within the consents) and to agree to the terms and conditions of the DAA. The DAA must be signed by the applicant and the relevant Head of Department, or equivalent. If applications include a named collaborator then the collaborator's Institution must sign and submit a separate DAA. A template DAA can be found on the EGA website: https://www.ebi.ac.uk/ega/sites/ebi.ac.uk.ega/files/documents/ Example%20DAA.doc.