SupplementCollins et al.
Supplement for “Identifying bipolar disorder susceptibility loci in a densely affected pedigree” by Collins et al.
Table of Contents
Supplemental Methods
Clinical......
Figure S1: Spanish bipolar pedigree......
Genomic technologies......
Table S1: Summary of individuals used for each technology......
Data analysis......
Supplementary Results
Regions shared IBD......
Table S2: Candidate genomic regions shared by BIP cases......
Identification of genetic variation in genomic regions shared by cases......
Genome-wide variant search......
Figure S2: Segregation of MBD1 variant in Spanish bipolar family......
Figure S3: Segregation of NR4A1 variant in Spanish Bipolar family......
Table S3: SNVs from high throughput sequencing......
Table S4: Validated CNVs not in 1000 genomes data
References......
Supplemental Methods
Clinical
Figure S1: Spanish bipolar pedigree
Samples were ascertained from Spain. The pedigree has a high prevalence of mood disorders, particularly bipolar disorder (type 1) and recurrent major depressive disorder(Figure S1).Best estimate DSM-IV diagnoses were made by three independent raters and based on an interview with the Schedule for Affective Disorders and Schizophrenia–Lifetime Version (SADS-L)1as well as on information from other family members.All participating persons gave written informed consent before participating in the study and the study was approved by the relevant ethical committees.
Genomic technologies
All genomic positions are given with reference to UCSC hg19. A summary of the technologies applied is given in Table S1.
Table S1:Summary of individuals used for each technology.
Technology / Purpose / Number of subjectsMicrosatellites / Genome-wide linkage / 14BIP; cases unaffected (19meioses)
Illumina Omni / Shared segments (IBD), CNVs / 11 (BIP cases)
Exome sequencing / Exonic SNPs / 5 (BIP cases)
Genome sequencing / SNPs, CNVs / 3 (BIP cases)
Selection of individuals. Selection of the 11 individuals for SNP array genotyping was based on individuals affected with BP1 or BP2 from across the pedigree, intentionally including more distant relatedness, to reduce regions shared IBD by chance. The individuals selected for exome (N=5) and whole genomesequencing (N=3) were chosen based on sex (female) and BIP I diagnostic status to reduce heterogeneity and including more distant relatedness.
SNP genotyping arrays.SNP genotyping was performed on DNA extracted from whole blood from11 individuals affected with either Bipolar I or Bipolar II from this pedigree.Illumina HumanOmni Quad chips (San Diego, California) with 1,140,419 SNPs were used.Genotypes were called using GenomeStudio, and we used PLINK2for quality control.We removed SNPs with no variation among the eleven individuals and removed SNPs if missing in all subjects.All subjects passed QC, and 647,415 SNPs passed QC.
High-throughput sequencing.For exome sequencing, capture was performed using Agilent SureSelect (all exome, version 2, Santa Clara, CA).Sequencing was performed on Illumina’s Genome Analyzer IIx system with paired-end 76base pair reads.Whole genome sequencing was also performed on Illumina’s HiSeq2000 instrument using TruSeq version 3 chemistry with paired-end 100bp reads. All procedures were per the manufacturers’ protocols.
Data analysis
Genome-wide linkage analysis.Microsatellite genotyping was described in Schumacher et al.3Data were available for three founders and eleven non-founders, resulting in19 effective meioses that could be incorporated into the analysis. Analysis was performed using a non-parametric approach and calculating non-parametric linkage (NPL) scores and P-values as implemented in Genehunter.4 Individuals with BIP1, BIP2, or schizoaffective disorder-bipolar type were considered affected (definition “ASDII” in the original report).
Shared segment analysis.SNP genotypes from Illumina HumanOmni arrays for 11 BIP cases were used to identify DNA segments shared in all or most cases. We used four different methods to identify genomic regions shared by cases: Germline (pairwise regions shared over >0.5cM)5, Beagle (pairwise sharing)6, PEDIBD7, andthe method ofThomas et al.8 (assumes large stretches of IBS are likely IBD andidentifies the largest segments shared by cases). Any regions identified by two or more methods were considered candidate regions.
Sequencing QC, alignment, variant calling, and filtering. QC began with evaluation of the rawfastq files using FASTQC ( followed by alignment to UCSC hg19 using BWA.9We followed Genome Analysis Toolkit (GATK) next-generation sequencing data workflows for additional QCand variant calling.10 In brief, we removed duplicated reads from the mapped reads and performed local realignment around indels and base quality recalibration.
Calling of alternate alleles for single nucleotide variants (SNV) was performed with the GATK Unified Genotyper module which callssingle nucleotide and insertion/deletion variants based on multiple samples. It uses a Bayesian genotype likelihood model to estimate simultaneously the most likely genotypes and allele frequency in a population of N samples, emitting an accurate posterior probability of there being a segregating variant allele at each locus as well as for the genotype of each sample.
For SNVs passing GATK QC, additional filtering was performed using Annovar.11 Annovar’s variant filtration portion labeled SNVs according to locations (genic or non-genic) and predicted functional consequences (synonymous, non-synonymous, or splice site). It additionally filtered against dbSNP130 and 1000 Genomes Project data (11/2010 release). Manual culling was also performed based on viewing of reads in IGV12 and searching dbSNP132 to remove non-novelSNVsand ones with evidence of false calls.
Copy number variation (CNV).To identify CNVs present in all or most affected individuals, CNV calls were made from the Illumina arrays on 11 individuals using PennCNV.13For QC, we excluded all CNV calls with: >50% overlap with genomic gaps (e.g., centromeres) or known segmental duplications; PennCNV confidence score <10; orsupported by <10 probes. These steps were intended to remove CNVs of lower confidence. The resulting CNV call set rangedhad 185-293 CNVsper subject.
To verify the presence of the CNVs detected by PennCNV, CNVs were also evaluated using whole genome sequence data. We merged the .bam files for the three individuals from whole genome sequencing, then called CNVs from this files usingforestSV.14
Risk profile scoring.Common variant risk profile or polygenic score analysis was conducted using PLINK as previously described.15Briefly, we used the Psychiatric Genomics Consortium BIP mega-analysis results16 ( as the discovery set, and computed per subject risk profiles scores in target samples: BIP cases from the GAIN study17, GAIN controls17, and 11 BIP cases from the Spanish pedigree. These per subject scores were weighted sums of the number of BIP risk alleles.Analysis was performed using 10,196 SNPs common to all subjects.Comparisons between the three groups were performed using PROC GLM in SAS.
Supplementary Results
Regions shared IBD
Genome-wide linkage analysis.Based on 435 microsatellite genotyped in 13 individuals,3 this pedigree showed a maximum multipoint NPL score of 6.1 (p = 9.8x10-4) which does not withstand multiple testing correction.We selected regions of potential interest based on p0.01 for the NPL results (Table S2).
IBD in BIP affecteds. The genome-wide linkage analyses included only a subset of the pedigree. We therefore applied a different approach in which we used SNP genotypes of only BIP cases to identify regions shared IBD. We used 647,415SNP genotypes from Illumina arrays on 11 BIP case. We used four different approaches to sharing to identify regions predicted by multiple methods. As each method had advantages and limitations, we ran four methods, and included regions identified by two or moremethods.
Beagle and Germline evaluated pairwise sharing (i.e., a maximum of 55 pairs in the 11 BIP cases) to identify genomic regions with increased sharing. The maximum pairs sharing for either test was 35 pairs (equivalent to a maximum of 9 of 11 individuals). We also searched for genomic regions with runs of identity-by-state in all 11 BIP cases and identified chr15:82.2-84.6Mb) as having the longest IBS run of markers (204 SNPs over 2.4Mb). This same region also was highlighted by PEDIBD. Chromosome 11 also has a region of potential sharing in all 11 individuals. While this region includes the centromere, the high degree of sharing is based on ~1500 markers, and it is possible this is not a false positive.
Candidate regions.To determine candidate regions, we included the linkage regions with p<0.01 and the top findings from≥2 of the following methods: PEDIBD>45 pairs sharing; runs of IBS>100 consecutive SNPs IBS in all individuals; Germline>25 pairs sharing; and Beagle>15 pairs sharing. The six genomic regions (93.7Mb)used as candidate regions are shown in Table S2.
Table S2: Candidate genomic regions shared by BIP cases
Method / Subjects sharing / Chr / bp1 / bp2 / NoteBeagle
Runs of IBS / 6-7
11 / 6 / 25,576,377 / 29,118,813 / Extended MHC
Genome-wide linkage / N/A / 6 / 77,461,313 / 104,674,073
Genome-wide linkage / N/A / 10 / 103,196,234 / 122,742,912
Beagle
PEDIBD / 6
11 / 11 / 48,050,995 / 56,611,798 / Overlaps centromere
Germline
Beagle / 7-8
6-7 / 13 / 80,445,558 / 105,078,780
Runs of IBS
PEDIBD / 11
10-11 / 15 / 81,636,267 / 85,282,071
Genome-wide linkage / N/A / 22 / 17,859,280 / 24,488,718 / Includes 22q11.21VCF region
Identification of genetic variation in genomic regions shared by cases
SNVs.To identify possible functional variants within the candidate regions in Table S2, we used data from exome sequencing of five BIP cases, and whole genome sequencing from three of the same cases. There were 26 novel, potentially functional SNVsdetected in all three individuals withwhole genome sequencing after filtering for the regions in Table S2. All of these SNVs were within three olfactory receptors (OR4C3, OR4C45, and OR9G9).We repeated this process for the exome sequencing data, and identified 21 SNVs (17 were also identified by whole genome sequencing), and all were again in olfactory genes (OR4C3, OR9G9 and OR8U8/OR8U1).
We examined the read alignments for these SNVs using IGV. Given that we are aligning short sequences (76bp for exome and 100bp for whole genome sequencing), it is possible for reads to only have a few basepairs different between one gene and another from the same gene family or an unannotated pseudogene of similar sequence.When this happens, existing alignment algorithms can place the reads in the wrong place.Because a gene which truly has many SNVs will have some of the alternate alleles on the maternal chromosome and some on the paternal, one can look at the reads for closely positioned SNVs to determine phasing. Misalignments can be indicated by clusters of SNVs and when both reads which carry multiple alternate allelesand overlapping reads with no called alternate alleles.Because multiple variants were present in each gene, the existence of other highly similar members of this large gene family, and alignment patterns, we concluded these SNVs were likely due to misalignment. However some degree of tolerance of SNPs within olfactory receptors leading to this clustering cannot be excluded.
CNVs. We used the Illumina arrays to look for CNVs in the genomic regions in Table S2. No deletions >1kb were identified within these candidate region. One smaller deletion (839bp) was identified as a homozygous deletion in all 11BIP cases by PennCNV, and was also called using the whole genome sequencing data. This deletion is not reported in the 1000 Genomes Project CNV deletion calls, overlaps with calls in the Database of Genomic Variation (DGV).18This small deletion is in a non-genic region, and is unlikely to have functional consequences.
We show below all CNVs present in ≥9 BIP cases that were not identified in HapMap were confirmed using whole genome sequencing data. None of these CNVs include genic coding regions.
Genome-wide variant search
Given that prioritization of variants in Table S2 did not identify compelling genetic variation of plausible impact on risk for BIP, we widened our search to the genomic regions assessed by the technologies used in this study.
SNVs.We used whole genome sequencing data to identify functional, novel variants in all three BIP cases. We removed non-genic or synonymous SNVs along withthose in the 1000 Genomes Project (11/2010 release)or in dbSNP v130. After filtering for variants found in all three BIP cases and removing the olfactory receptor SNVs in OR4C3,OR4C45, andOR11H1, 34 variants remained(Table S3). We annotated these variants using the exome sequencing results and identified suspect SNV calls by looking at read alignment using IGV.
We prioritized eight variants for verification using Sanger sequencing (shown in green in Table S3).SNVsin ADAMTS7 and PRSS1were not confirmed by Sanger sequencing.Of the six SNVs confirmed by Sanger sequencing, the MLL2SNV was present in only four of 11BIP cases sequenced, and was excluded from further consideration. The remaining five variants were sequenced in 42 individuals in the pedigree for whom DNA was available. Three SNVs showed random patterns of segregation (NRXN1,ERMP1, and MYO9A). Two SNVs(MBD1 and NR4A1)showed patterns of potential interest (Figure S2-S3).Due to lack of informative matings in this pedigree, linkage scores were not calculable.
Figure S2: Segregation of MBD1 variant in Spanish bipolar family.
Figure S3: Segregation of NR4A1 variant in Spanish Bipolar family.
Page 1
SupplementCollins et al.
Table S3: SNVs from high throughput sequencing
WGS / ExomeGene / location / 042 / 005 / 105 / 042 / 005 / 105 / 012 / 107 / REF / ALT / VAR1 / Observations
FSHR / chr2:49295422 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/0 / 0/0 / C / T / nonsyn / Exome shows ref/ref individual
NRXN1 / chr2:50724576 / 0/1 / 0/1 / 0/1 / G / A / nonsyn
MUC4 / chr3:195509795 / 0/1 / 0/1 / 0/1 / C / T / nonsyn / Large deletion in dbSNP; other SNPs nearby
MUC4 / chr3:195512186 / 1/1 / 1/1 / 1/1 / T / C / nonsyn / Indel in dbSNP132
DDIT4L / chr4:101108877 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/0 / T / C / nonsyn / >1% dbSNP
PRIM2 / chr6:57398226 / 0/1 / 0/1 / 0/1 / T / A / Nonsyn / >1%dbSNP; >2 haplotypes
PRSS1 / chr7:142460335 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / A / G / nonsyn
ERMP1 / chr9:5832782 / 0/1 / 0/1 / 0/1 / 0/1 / 0/0 / 1/1 / ./. / 1/1 / G / C / nonsyn / 0/0 call too low coverage, may not be 0/0
FAM35A / chr10:88950273 / 1/1 / 0/1 / 0/1 / 1/1 / 0/1 / 0/1 / 0/0 / 0/0 / G / A / nonsyn / Exome shows ref/ref individual
MUC6 / chr11:1016665 / 0/1 / 0/1 / 0/1 / C / T / nonsyn / >1% dnSNP132
MUC6 / chr11:1017338 / 0/1 / 0/1 / 0/1 / C / A / nonsyn / >1% dnSNP132
MUC6 / chr11:1017691 / 0/1 / 0/1 / 0/1 / C / G / nonsyn / >1% dbSNP132
MUC6 / chr11:1017858 / 0/1 / 0/1 / 0/1 / A / G / nonsyn / >1% dnSNP132
OVCH2 / chr11:7717219 / 1/1 / 1/1 / 1/1 / A / T / splicing / Alignment questionable
MLL2 / chr12:49434717 / 0/1 / 0/1 / 0/1 / C / T / nonsyn
TROAP / chr12:49717507 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/0 / C / A / nonsyn / Exome shows ref/ref individual
CELA1 / chr12:51740409 / 1/1 / 1/1 / 1/1 / 1/1 / 0/1 / 1/1 / 1/1 / 0/0 / T / G / Nonsyn / >1%dbSNP132
CELA1 / chr12:51740410 / 1/1 / 1/1 / 1/1 / 1/1 / 0/1 / 1/1 / 1/1 / 0/0 / A / G / nonsyn / >1%dbSNP132
NR4A1 / chr12:52452495 / 0/1 / 0/1 / 0/1 / 0/1 / ./. / 0/0 / 1/1 / 0/1 / C / T / nonsyn / 0/0 call too low coverage, may not be 0/0
KRT82 / chr12:52799923 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/0 / G / A / stopgain SNV / Exome shows ref/ref individual
KRT78 / chr12:53238386 / 0/1 / 1/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/0 / G / A / nonsyn / Exome shows ref/ref individual
MYO9A / chr15:72192205 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / C / G / nonsyn
ADAMTS7 / chr15:79058730 / 0/1 / 1/1 / 0/1 / 0/1 / 1/1 / 0/1 / 0/1 / 0/1 / A / G / nonsyn
KCNJ12 / chr17:21318782 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / G / A / nonsyn / >1% dbSNP132;>2 haplotypes
KCNJ12 / chr17:21319069 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / G / A / nonsyn / >1% dbSNP132
KCNJ12 / chr17:21319087 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / 0/1 / G / A / nonsyn / >1% dbSNP132
TBC1D3C / chr17:34587732 / 1/1 / 1/1 / 1/1 / C / G / nonsyn / mapping quality poor
TBC1D3C / chr17:34587756 / 1/1 / 1/1 / 1/1 / T / C / nonsyn / mapping quality poor
TBC1D3C / chr17:34587758 / 1/1 / 1/1 / 1/1 / G / A / nonsyn / mapping quality poor
MBD1 / chr18:47793974 / 0/1 / 0/1 / 0/1 / G / C / nonsyn
ADAMTSL5 / chr19:1506135 / 0/1 / 0/1 / 0/1 / 0/0 / 0/1 / 0/0 / 0/1 / 0/0 / C / T / nonsyn / Exome shows ref/ref individual
FCGBP / chr19:40392631 / 1/1 / 0/1 / 0/1 / 1/1 / 0/1 / 0/1 / 1/1 / 0/1 / G / A / nonsyn / >1% dbSNP132
ZNF814 / chr19:58385748 / 0/1 / 0/1 / 0/1 / G / A / nonsyn / Uneven allelic balance
C21orf62 / chr21:34166190 / 1/1 / 1/1 / 1/1 / 1/1 / 1/1 / 1/1 / 1/1 / 1/1 / A / T / nonsyn / dbSNP132
WGS=whole genome sequencing. The numbers in the second row are pedigree ID numbers. The genotype is represented as a 0 for reference allele and 1 for an alternate allele.Nonsyn represents a nonsynonymous variant and stopgain is a variant which creates a stop codon.
Page 1
SupplementCollins et al.
CNVs.We filtered all CNV calls to retain novel CNVs (i.e., not in 1000 Genomes Project) that were either deletions or duplications (not both)in ≥8 of 11 BIP cases, >1kb in size, covered ≥1 exon, and verified by whole genome sequencing.No CNVs remained after filtering.
Restarting with the initial CNV call set, several CNVs were present in ≥9 of 11 BIP cases, verified using whole genome sequencing, not presentin 1000 Genomes data, but all were reported in DGV.Additionally, the only one that overlaps any coding genes was present in the deleted and the duplicated form in the family (Table S4).
Table S4: Validated CNVsnot in 1000 genomes data°
Chr / Copy number* / Probes / Size (kb) / # of BIP cases / Geneschr1:158867802-158871453 / 0 (11 indiv.) / 15 / 3.65 / 11 / None
chr1:187716675-187722124 / 1 (8 indiv.);
0 (3 indiv.) / 14 / 5.45 / 11 / None
chr3:32102055-32106725 / 1 (2 indiv.)
0 (7 indiv.) / 12 / 4.67 / 9 / ENSEMBL unprocessed pseudogene
chr8:47406167-47407676 / 1(2 indiv.)
0 (7 indiv) / 19 / 1.51 / 9 / None
chr9:131412549-131452672 / 0(7 indiv)
3 (4 indiv.) / 13 / 40.12 / 9 / WDR34; SET
chr10:107950711-107951550 / 0 (11 indiv.) / 14 / 0.84 / 11 / Intron of ENSEMBL lincRNA
chr16:76540062-76543447 / 0 (10 indiv.) / 15 / 3.39 / 10 / Intro of CNTNAP4
*Copy number is given, with the number of individuals with that copy number in parenthesis
° All CNVs in this table overlap a region reported in DGV
References
1.Endicott J, Spitzer RL. A diagnostic interview: the schedule for affective disorders and schizophrenia. Arch Gen Psychiatry 1978; 35(7): 837-844.
2.Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet 2007; 81(3): 559-575.
3.Schumacher J, Kaneva R, Jamra RA, Diaz GO, Ohlraun S, Milanova V et al. Genomewide scan and fine-mapping linkage studies in four European samples with bipolar affective disorder suggest a new susceptibility locus on chromosome 1p35-p36 and provides further evidence of loci on chromosome 4q31 and 6q24. Am J Hum Genet 2005; 77(6): 1102-1111.
4.Kruglyak L, Daly MJ, Reeve-Daly MP, Lander ES. Parametric and nonparametric linkage analysis: a unified multipoint approach. Am J Hum Genet 1996; 58(6): 1347-1363.
5.Gusev A, Lowe JK, Stoffel M, Daly MJ, Altshuler D, Breslow JL et al. Whole population, genome-wide mapping of hidden relatedness. Genome Res 2009; 19(2): 318-326.
6.Browning SR, Browning BL. High-resolution detection of identity by descent in unrelated individuals. Am J Hum Genet 2010; 86(4): 526-539.
7.Li X, Yin X, Li J. Efficient identification of identical-by-descent status in pedigrees with many untyped individuals. Bioinformatics 2010; 26(12): i191-198.
8.Thomas A, Camp NJ, Farnham JM, Allen-Brady K, Cannon-Albright LA. Shared genomic segment analysis. Mapping disease predisposition genes in extended pedigrees using SNP genotype assays. Ann Hum Genet 2008; 72(Pt 2): 279-287.
9.Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 2009; 25(14): 1754-1760.
10.DePristo MA, Banks E, Poplin R, Garimella KV, Maguire JR, Hartl C et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat Genet 2011; 43(5): 491-498.
11.Wang K, Li M, Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res 2010; 38(16): e164.
12.Robinson JT, Thorvaldsdottir H, Winckler W, Guttman M, Lander ES, Getz G et al. Integrative genomics viewer. Nat Biotechnol 2011; 29(1): 24-26.
13.Wang K, Li M, Hadley D, Liu R, Glessner J, Grant SF et al. PennCNV: an integrated hidden Markov model designed for high-resolution copy number variation detection in whole-genome SNP genotyping data. Genome Res 2007; 17(11): 1665-1674.
14.Michaelson JJ, Sebat J. forestSV: structural variant discovery through statistical learning. Nat Methods 2012; 9(8): 819-821.
15.Purcell SM, Wray NR, Stone JL, Visscher PM, O'Donovan MC, Sullivan PF et al. Common polygenic variation contributes to risk of schizophrenia and bipolar disorder. Nature 2009; 460(7256): 748-752.
16.Large-scale genome-wide association analysis of bipolar disorder identifies a new susceptibility locus near ODZ4. Nat Genet 2012; 44(9): 1072.
17.Zhang D, Cheng L, Qian Y, Alliey-Rodriguez N, Kelsoe JR, Greenwood T et al. Singleton deletions throughout the genome increase risk of bipolar disorder. Mol Psychiatry 2009; 14(4): 376-380.
18.Iafrate AJ, Feuk L, Rivera MN, Listewnik ML, Donahoe PK, Qi Y et al. Detection of large-scale variation in the human genome. Nat Genet 2004; 36(9): 949-951.
Page 1