Genome folding in evolution and disease

4.2 Materials and Methods

4.2.1 Selection of subjects with apparently balanced chromosome abnormalities

BCA breakpoints and clinical data were obtained from DGAP cases for which whole-genome sequencing was performed using a previously described large-insert jumping library approach.(Higgins et al. 2008; Ligon et al. 2005; Kim and Marcotte 2008; Lu et al. 2007; Redin et al. 2017; Talkowski et al. 2011) A total of 151 cases were filtered to select only subjects whose translocation or inversion breakpoints fall within intergenic regions (GRCh37) and did not overlap known long intergenic non-coding RNAs (lincRNAs) or pseudogenes, as these elements have been shown to exert functional roles (reviewed in (Quinn and Chang 2016) (Pink et al. 2011; Muro and Andrade-Navarro 2010)). Of 151 DGAP subjects, only 17 fulfilled our selection criteria, 12 of whom had available and reportedly normal clinical array results, suggesting lack of large duplications or deletions.

4.2.2 Clinical descriptions of DGAP cases

The clinical presentation of the 17 subjects varied, ranging from developmental delay to neurological conditions, offering the opportunity to assess long-range position effects in different phenotypes. Subjects’ karyotypes are presented in the main text using the International System for Human Cytogenetic Nomenclature (ISCN2016) (Table 4.1). Detailed case descriptions are included in the Supplemental Note: Case Reports, as well as a nomenclature developed to describe chromosome rearrangements using next-generation sequencing.(Ordulu et al. 2014) Reported ages of DGAP subjects are from time of enrollment. All reported genomic coordinates use GRCh37.

Table 4.1: **Description of the 17 analyzed DGAP cases with non-coding BCAs.** Corresponding clinical karyotypes are reported, with overlap of breakpoints with regulatory elements (E = enhancer, DHS = DNaseI hypersensitive sites, CTCF = CTCF binding sites), and TADs from H1-hESC, IMR90, and GM12878 (1= one breakpoint within TAD, 2=both BCA breakpoints are located within TAD). Top-ranking position effect genes are provided for the +-1 Mb windows surrounding the BCA breakpoints; each gene is highlighted with different evidence supporting its inclusion (a = ClinGen known recessive genes, b= ClinGen genes with emerging and sufficient evidence suggesting haploinsufficiency is associated with clinical phenotype, c = HI scores less than 10, d = within H1-ESC TAD, e = DHS enhancer-promoter disrupted interactions).
Subject ID and Reported Karyotype	Disruption of Functional Elements	Breakpoints within TADs (hESC / IMR90 / GM12878)	Top-ranking Candidates +-1 Mb
DGAP017 46,X,t(X;10)(p11.2;q24.3)	DHS	2/2/1	-
DGAP111 46,XY,t(16;20)(q11.2;q13.2)dn	CTCF	1/1/2	ORC6^a
DGAP113 *46,XY,t(1;3)(q32.1;q13.2)dn	-	2/2/2	ASPM^a
DGAP126 46,XX,t(5;10)(p13.3;q21.1)dn	-	2/1/2	-
DGAP138 46,XY,t(1;6)(q23;q13)dn	-	2/2/2	GRIK2^ac
DGAP153 46,X,t(X;17)(p11.23;p11.2)dn	-	1/1/1	-
DGAP163 46,XY,t(2;14)(p23;q13)dn	-	2/2/2	SOS1^cde, COCH^de
DGAP176 46,Y,inv(X)(q13q24)mat	DHS, CTCF	2/1/2	ACSL4^bd, COL4A5^bcde
DGAP249 46,XX,t(2;11)(q33;q23)dn	E, DHS	2/2/2	SATB2^bcde, SORL1^e
DGAP252 46,XY,t(3;18)(q13.2;q11.2)dn	-	2/2/2	RBBP8^a,GATA6^bcde
DGAP275 46,XX,t(7;12)(p13;q24.33)dn	DHS	1/1/2	ANKLE2^e, POLE^e
DGAP287 46,XY,t(10;14)(p13;q32.1)dn	CTCF	2/2/2	-
DGAP288 46,XX,t(6:17)(q13;q21)dn	DHS	2/2/2	SOX9^bcd
DGAP315 46,XX,inv(6)(p24q11)dn	-	1/1/2	-
DGAP319 46,XX,t(4;13)(q31.3;q14.3)dn	-	2/1/2	-
DGAP322 46,XY,t(1;18)(q32.1;q22.1)	DHS	1/2/2	IRF6^bcd
DGAP329 46,XX,t(2;14)(q21;q24.3)dn	-	1/2/2	ZEB2^bcde

4.2.3 Analysis of genes bordering the rearrangement breakpoints

The presence of annotated genes or pseudogenes and lincRNAs was assessed in windows of +-3 and +-1 Mb neighboring each subject’s translocation and inversion breakpoints, and within reported H1-hESC topologically associated domains (TADs)(Dixon et al. 2012) where the breakpoints were located. The gene annotation file was obtained from Ensembl GRCh37 archive,(Flicek et al. 2014) and we used the Human Body Map lincRNAs catalog.(Cabili et al. 2011) Haploinsufficiency (HI) and triplosensitivity scores were assigned using Huang et al., 2010(Huang et al. 2010) and version hg19 of ClinGen(Rehm et al. 2015) data downloaded on 9/20/2016.

4.2.4 Assessment of disrupted functional elements and chromatin interactions bordering rearrangement breakpoints

The disruption of regulatory elements such as enhancers, promoters, locus control regions, and insulators can lead to disease-related gene expression changes; DNase I hypersensitive (DHS) sites have been used as markers for the identification of such elements.(Thurman et al. 2012) In addition, the alteration of TAD boundaries has been previously shown to cause a rewiring of enhancers with pathological consequences;(Lupiáñez et al. 2015; Giorgio et al. 2015; Narendra et al. 2015) CCCTC-Binding Factor (CTCF) binding sites have been found to be enriched in TAD boundaries,(Dixon et al. 2012) and several mutations of boundary-defining sites have been associated with cancer.(Flavahan et al. 2016; Hnisz et al. 2016 b) Based on these observations, we assessed the number of regulatory elements that were potentially disrupted by the analyzed DGAP breakpoints. We compared the breakpoint positions of the selected DGAP subjects against data corresponding to CTCF binding sites, DHS sites, and chromatin segmentation classifications (Broad ChromHMM) derived from a lymphoblastoid cell line (GM12878) and human stem cells (H1-hESC), obtained from the Encyclopedia of DNA Elements (ENCODE) project(Dunham et al. 2012) and accessed through the University of California Santa Cruz Genome Browser.(Kent et al. 2002) Enhancer positions were additionally obtained from Andersson et al., 2014(Andersson et al. 2014) for tissue and primary cells, and the VISTA Enhancer browser, human version hg19.(Visel et al. 2007) Finally, lists of transcription factor (TF) binding sites and gene promoters were obtained from the Ensembl database human version GRCh37.(Flicek et al. 2014) Hi-C interaction data and TAD positions for H1-hESC, GM06990, and IMR90 at 20 Kb, 40 Kb, 100 Kb, and 1 Mb resolution were obtained from Dixon et al., 2012(Dixon et al. 2012) and the WashU EpiGenome Browser.(Zhou and Wang 2012) A high-resolution dataset of chromatin loops and domains was obtained from Rao et al., 2014 for IMR90 and GM12878 cells.(Rao et al. 2014) Lastly, distal DHS/enhancer–promoter connections(Thurman et al. 2012) were used to assess disrupted predicted cis-regulatory interactions by the BCAs. Genomic overlaps between the rearrangement breakpoints, functional elements and disrupted chromatin interactions were calculated using custom Perl scripts, the BEDtools suite(Quinlan and Hall 2010) and the genomic association tester (GAT) tool.(Heger et al. 2013)

4.2.5 Ontological analysis of genes neighboring breakpoints

Phenotype similarity between potential position effect genes and DGAP cases was calculated by converting the phenotypes of the 17 subjects to Human Phenotype Ontology (HPO)(Köhler et al. 2014) terms and calculating their phenomatch score as described in Ibn-Salem et al., 2014.(Ibn-Salem et al. 2014) The phenomatch score quantiﬁes the information content of the most speciﬁc HPO term that is part of or a common ancestor (more general term) of a set of phenotypes. Our set of phenotypes is constituted by the HPO terms associated to DGAP cases and the ones annotated to candidate position effect genes within windows of +-3 and +-1 Mb of sequence in proximity to the breakpoints. We used two background models to assess signiﬁcance of this similarity. The ﬁrst is based on randomly permuting the associations of phenotypes to genes; to this effect, the phenotype-gene associations are shuffled 100 times randomly and the similarity of these random phenotypes to the studied case clinical findings is calculated. The second background control is based on shifting the breakpoint location along the chromosome; each breakpoint is shifted by -9, -6, -3, +3, +6, and +9 Mb and the similarity of genes in proximity to the shifted breakpoints is computed.

4.2.6 Quantitative real-time PCR

LCLs derived from DGAP236-02m, DGAP244-02m and DGAP245-02m were used as karyotypically normal male controls. These are karyotypically normal fathers of enrolled DGAP cases with no history of disease. LCL 17402 (DGAP163) was used to test differential gene expression for SOS Ras/Rac guanine nucleotide exchange factor 1 (SOS1 [MIM: 182530]), and LCL 18060 (DGAP176) was used to test midline 2 (MID2 [MIM: 300204]), p21 (RAC1) activated kinase 3 (PAK3 [MIM: 300142]), and POU class 3 homeobox 4 (POU3F4 [MIM: 300039]) expression using quantitative polymerase chain reaction (qPCR). Glucuronidase beta (GUSB [MIM: 611499]) was used as a housekeeping control. qPCR experiments were performed by the Harvard Biopolymers Facility using TaqMan probes Hs00264887_s1 (POU3F4), Hs00201978_m1 (MID2), Hs00176828_m1 (PAK3), Hs00893134_m1 (SOS1), and Hs00939627_m1 (GUSB). Data were analyzed using the ΔCT method.

4.2.7 Assessment of DGAP breakpoints overlapping with non-coding structural variants in public databases

To find similar non-coding structural rearrangement subjects and compare their annotated clinical phenotypes to those observed in DGAP cases, we searched the DatabasE of genomiC varIation and Phenotype in Humans using Ensembl Resources (DECIPHER)(Firth et al. 2009) version 2015-07-13, as well as the dbVar database from the National Center for Biotechnology Information (NCBI) Variation Viewer 1.5.(Lappalainen et al. 2012) Both databases are comprehensive community-supported repositories of clinical cases with novel and extremely rare genomic variants.