ARTICLE
https://doi.org/10.1038/s41467-020-15022-4 OPEN
Prioritizing disease and trait causal variants
at the TNFAIP3 locus using functional and
genomic features
John P. Ray1,12, Carl G. de Boer1,2,12, Charles P. Fulco 1,3, Caleb A. Lareau 1,4, Masahiro Kanai 1,5,6,
Jacob C. Ulirsch 1,4, Ryan Tewhey 1,7, Leif S. Ludwig1, Steven K. Reilly1,7, Drew T. Bergman 1,
Jesse M. Engreitz 1,8, Robbyn Issner1, Hilary K. Finucane 1,5, Eric S. Lander 1,3,9, Aviv Regev 1,2,9,10,13✉ &
Nir Hacohen 1,11,13✉
Genome-wide association studies have associated thousands of genetic variants with com-
plex traits and diseases, but pinpointing the causal variant(s) among those in tight linkage
disequilibrium with each associated variant remains a major challenge. Here, we use seven
experimental assays to characterize all common variants at the multiple disease-associated
TNFAIP3 locus in five disease-relevant immune cell lines, based on a set of features related to
regulatory potential. Trait/disease-associated variants are enriched among SNPs prioritized
based on either: (1) residing within CRISPRi-sensitive regulatory regions, or (2) localizing in a
chromatin accessible region while displaying allele-specific reporter activity. Of the 15 trait/
disease-associated haplotypes at TNFAIP3, 9 have at least one variant meeting one or both of
these criteria, 5 of which are further supported by genetic fine-mapping. Our work provides a
comprehensive strategy to characterize genetic variation at important disease-associated
loci, and aids in the effort to identify trait causal genetic variants.
1 Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA. 2 Klarman Cell Observatory, Broad Institute of MIT and Harvard, Cambridge, MA 02142,
USA. 3Department of Systems Biology, Harvard Medical School, Boston, MA 02115, USA. 4 Program in Biological and Biomedical Sciences, Harvard Medical
School, Boston, MA 02115, USA. 5 Analytic and Translational Genetics Unit, Massachusetts General Hospital, Boston, MA 02114, USA. 6 Program in
Bioinformatics and Integrative Genomics, Harvard Medical School, Boston, MA 02115, USA. 7Department of Organismic and Evolutionary Biology, Harvard
University, Cambridge, MA 02138, USA. 8Harvard Society of Fellows, Harvard University, Cambridge, MA 02138, USA. 9Department of Biology,
Massachusetts Institute of Technology, Cambridge, MA 02142, USA. 10 Howard Hughes Medical Institute, Cambridge, MA 02142, USA. 11 Center for Cancer
Research, Massachusetts General Hospital, Boston, MA 02114, USA. 12These authors contributed equally: John P. Ray, Carl G. de Boer. 13These authors jointly
supervised this work: Aviv Regev, Nir Hacohen. ✉email: aregev@broadinstitute.org; nhacohen@mgh.harvard.edu
NATURE COMMUNICATIONS |         (2020) 11:1237 | https://doi.org/10.1038/s41467-020-15022-4 | www.nature.com/naturecommunications 1
1234567890():,;
ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15022-4
Genome-wide association studies (GWAS) have revealed As a proof of concept, we optimize and apply seven assays to>100,000 associations of genetic variants with human characterize all known common genetic variants in the TNFAIP3traits and diseases (e.g. autoimmune disease), but it locus, a genetic locus associated with multiple autoimmune dis-
remains a challenge to pinpoint the causal variant(s) that account eases19, and where disease-associated genetic and epigenetic fea-
for the association by altering disease risk and determine their tures have been studied extensively20–24. We use cell lines derived
functions1–4. This is because they are often in tight linkage dis- from T cells, B cells, and monocytes (U937 or THP-1 monocyte
equilibrium (LD) with non-causal variants and, in the vast cell lines, GM12878 or BJAB B cell lines, or Jurkat T cell line),
majority of cases, lie in non-coding regions, where it is more representing three major cell lineages that can impact auto-
challenging to predict the impact and relevant context of variants. immunity. We find that two criteria are correlated with sig-
Most causal variants in the non-coding genome are likely to act nificant enrichment for the subset of SNPs that show disease/
through altering transcript abundance in a disease-relevant con- trait-association and, by inference, the subset of SNPs that play a
text. In the relevant context (cell type, tissue source, stimulation, causal role in these associations. These two criteria are: (i) loca-
genetic background, and disease status), experimental assays could lization within CRISPRi-sensitive regions in one of the cell types,
be used to characterize the relationship between genetic variants or (ii) localization within open chromatin regions while also
and gene regulation. However, there are several challenges in this showing allele-specific reporter activity by MPRA. We find SNPs
strategy. First, one or more aspects of the relevant context may be that fulfill at least one of these two criteria in 9 of 15 disease/trait-
unknown. Second, even in the relevant context, there are many associated TNFAIP3 haplotypes, prioritizing 18 putatively causal
possible impacts of non-coding variants (such as different effects SNPs in the locus associated to 15 diseases. By contrast, several
on gene expression or isoform usage), and each would involve a other criteria showed no enrichment for disease/trait association.
separate experimental assay, highlighting different features. Third, Our results highlight the limitations of using individual assays for
although ideally the relationship would be tested by allelic sub- implicating a variant as potentially functional, and suggests that a
stitution in the relevant context—for instance, by CRISPR- combination of assays, cell types and context will be needed to
directed base editing or homologous recombination5–8, this prioritize variants at disease loci.
approach is difficult to scale at present. As a result, various assays
have been proposed for identifying potentially causal variants,
based on the variant’s relation to or impact on different molecular Results
features in a relevant cell type. The TNFAIP3 locus harbors 15 independent disease associa-
These assays can be categorized into four classes, depending on tions. As a test case, we investigated the TNFAIP3 locus because it
(i) whether they involve observations of natural systems or has strong associations to many autoimmune diseases. TNFAIP3
engineered experimental perturbations and (ii) whether they encodes the A20 protein, which is upregulated by NF-kB upon
pertain to a region or an individual variant. immune stimulation, and dampens pathways that activate NF-kB
in a negative feedback loop (Fig. 1a)19,25,26. At least 49 GWASs
(1) Observational assays that characterize the genomic region have identified genome-wide significant SNPs in the TNFAIP3
in which the variant resides. Examples include using
ATAC-seq, DNase-I-seq, and H3K27ac ChIP-seq1,4,9,10
locus that together are associated with 16 human diseases and
, as phenotypes, including lupus (SLE), rheumatoid arthritis (RA),
well as testing whether the variant lies in spatial proximity psoriasis, inflammatory skin disorder (ISD), celiac disease,
to a target gene, based on topological assays such as 4C or
11,12 inflammatory bowel disease (IBD), and multiple sclerosis (MS).HiC . Rather than focusing only on disease-associated SNPs (that is,
(2) Observational assays that characterize the impact of those showing genome-wide-significant associations for one of
naturally occurring genetic differences at the variant. these diseases as tag SNPs or in tight LD to them), we system-
Examples include characterizing whether the variant shows atically examined all common SNPs (MAF > 0.01) in the ~300 kb
allele-specific association with expression of one or more topologically associating domain (TAD) containing TNFAIP3
nearby genes or with local chromatin features (that is, an (based on HiC data from GM12878 B cells and THP-1 monocyte
expression quantitative trait locus (eQTL) or a chromatin cell lines12,27), and 150 kb on either side of the TAD because it is
QTL, respectively), or whether the variant disrupts a known that regulatory regions can affect the expression of genes
transcription factor (TF) motif. outside of TADs28 (Fig. 1b, top; Supplementary Fig. 1). We rea-
(3) Engineered perturbational assays that test the impact of the soned that studying all common non-coding variants would allow
genomic region containing the variant. Examples include us to derive empirical null distributions for each assay because
assaying the effect of CRISPR-directed inhibition (e.g., most variants are not expected to be functional. Accordingly, we
CRISPRi13) and activation (e.g., CRISPRa14) of the region selected for analysis all 2776 common variants with minor allele
on the expression of nearby genes or on chromatin frequency > 0.01 in East Asian or European populations (in 1000
organization. Genomes, see “Methods” section).
(4) Engineered perturbational assays that test the impact of the We next analyzed the locus to estimate the number of SNPs
variant itself. Examples include testing allele-specific that contribute to disease. Of the 2776 variants, 294 were in tight
enhancer activities in massively parallel reporter assays LD (r215–18 > 0.8) to at least one of 34 ‘tag SNPs’—that is, a SNP(MPRAs) and related methods . reported as having the highest association score in one of the
These assays have been used in previous studies to suggest GWASs for the autoimmune and other diseases noted above
particular genetic variants as more likely to impact disease risk. (Fig. 1c; Supplementary Fig. 2a). Through LD analysis (r2 ≥ 0.8)
However, we do not know the extent to which each of these of the tag SNPs, we identified 15 independent haplotypes
assays actually enriches for causal variants. associated with one or more GWAS traits in Europeans (Fig. 1d;
Here, we reason that assays that usefully prioritize disease- Supplementary Fig. 2b–d); three of these haplotypes also
causal variants could be recognized by testing whether they overlapped East Asian disease-associated haplotypes, but with
effectively enrich for disease-associated variants among all var- slight differences in the associated SNPs (Fig. 1d; Supplementary
iants across a region. However, because disease causal variants for Fig. 2d). Notably, fine-mapping of immune-related UK Biobank
most associations are unknown, we use disease-associated var- phenotypes (autoimmune disease (self-reported or diagnosed),
iants (which are known and highly enriched for causal variants). self-reported allergy, and eosinophil counts) showed that, despite
2 NATURE COMMUNICATIONS |   (      2020)1 1:1237 | https://doi.org/10.1038/s41467-020-15022-4 | www.nature.com/naturecommunications
NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15022-4 ARTICLE
a Fig. 1 Disease variants in the complex autoimmune-associated TNFAIP3TNFa
locus. a TNFAIP3 encodes the A20 protein, which forms part of a negative
TNFR1
Proteasomal feedback loop to dampen NF-kB-mediated immune activation. b HiC plotsK48 Ubiquitin
degradation TRADDRIP1 for the lymphoblastoid B cell line GM12878, with color intensity
A20 proportional to the interaction frequency between genomic coordinates (x-
axis). Boxes indicate the 300 kb high-interaction domain and the 605 kb
K63 Ubiquitin
region used in this study. c, d Genetics of the TNFAIP3 locus. The positions
NF-kB (shared x-axis indicated above c) of variants with respect to the TNFAIP3
IKK
gene and a lncRNA (LOC100130476). c GWAS tag SNPs (red) and SNPs in
NF-kB nuclear translocation tight LD (greyscale boxes indicating LD to tag SNP) for many immune-
NF-kB
related phenotypes (y-axis). d GWAS haplotypes defined by combining all
SNPs in tight LD (r2 > 0.8) to GWAS tag SNPs for European (CEU; top) and
b 300 kbTNFAIP3 HiC high-interaction East Asian (CHB/JPT; bottom) populations. Colors are used to help identify
shared haplotypes between CEU and CHB/JPT populations. e Autoimmune
605 kb region
used for this study GWAS signals are enriched in open chromatin of immune cells. Heritability
enrichment (color) of disease-associated SNPs in DHS of various tissues
(x-axis) for seven autoimmune diseases (y-axis), according to LD-score
regression. Also see Supplementary Data 1, 2.
role in the autoimmune diseases with which the TNFAIP3 locus is
137,500 kb 139,000 kb associated because their localization in disease-associated tissues,
TNFAIP3 signaling, and function are correlated with disease progression in
the clinic and in animal models of disease29–34137,900 kb 138,000 kb 138,100 kb 138,200 kb 138,300 kb 138,400 kb . T cell-, B cell-,
LOC100130476 and monocyte-specific accessible chromatin and active histone
TNFAIP3
c marks (H3k27ac and H3K4me3 ChIP-seq) are also significantlyAllergy
Celiac LD (r 2)
IBD enriched (compared to other cell types) for GWAS variants (P <
ISD 1.0
MS 1 × 10−80.9 ) from studies of diseases that had associations in
Psoriasis
RA 0.8
SLE 0.7 TNFAIP3 according to stratified LD score regression35 (Fig. 1e;
Other Tag SNP Supplementary Fig. 3a–c). Moreover, deleting TNFAIP3 in these
cell types causes systemic autoimmunity in mice36–40.
d CEU We studied cell lines derived from these cell types: THP-1 and
15
14 U937 for monocytes, BJAB and GM12878 for B cells, and Jurkat
13
12 for T cells. The chromatin accessibility profiles of these cell lines
11
10
9 are enriched for autoimmune-associated risk variants similarly to
8
7 the corresponding primary cells (Supplementary Fig. 3d), and
6
5
4 among blood cell types profiled by ATAC-seq20 they were most
3
2 similar to the cell type they represent (Supplementary Fig. 4a),
1
CHB/JPT especially at the TNFAIP3 locus (Supplementary Fig. 4b),
9
6
2 suggesting that the selected cell lines could serve as models for
these cell types.
e GWAS enrichment in DNase I hypersensitive regions
Celiac disease
Crohn's disease A panel of assays to annotate genetic variation. We used both
Lupus observational and perturbational assays to characterize regulatory
Multiple sclerosis features in the regions where variants were located, and the
Rheumatoid Arthritis variants themselves (Fig. 2).
Type 1 diabetes
Ulcerative colitis Using observational assays, we first analyzed regions that
contact the TNFAIP3 promoter (primary T cell and GM12878 B
1
cell HiChIP data; ~5 kb resolution41) and regions of accessible
0.75
0.5 chromatin in any of the cells lines (using ATAC-seq in
0.25
0 unstimulated and stimulated cells (Supplementary Fig. 5a, b),
and publicly available DHS of cell types from the blood42). For
each variant, we also assessed whether it lies within a region
bound by a TF based on ChIP-seq42, and whether the variant is
predicted to affect TF binding according to its cognate motif
(Supplementary Fig. 5c).
Using perturbational assays, we sought to identify regions that
limited sample size, all but two of these separately fine-mapped can affect TNFAIP3 expression. With CRISPRi (in which KRAB-
alleles were contained on three of the 15 disease-associated dCas9 binds to a region targeted by a guide RNA and represses
haplotypes from our LD analysis (Supplementary Data 1, 2, see chromatin locally13), we identified regions whose inhibition alters
“Methods” section). Collectively, we estimate that at least 15 SNPs TNFAIP3 expression. We targeted all regions with accessible
in the locus contribute to disease. chromatin in either U937, BJAB, or Jurkat cell lines, tiled guides
While TNFAIP3 is likely to play a role in many disease-relevant across each element (and up to 100 bp on either side), and
cell types, we chose to study T cells, B cells, and monocytes. These identified guides and regions that significantly repress TNFAIP3
important innate and adaptive immune cell types likely play a expression (see the “Methods” section; Supplementary Fig. 6,
NATURE COMMUNICATIONS |       (  2020)1 1:1237 | https://doi.org/10.1038/s41467-020-15022-4 | www.nature.com/naturecommunications 3
Enrichment
Small intestine
Psoas muscle
Primary T cells (peripheral blood)
Primary T cells (cord blood)
Primary natural killer cells (peripheral blood)
Primary monocytes (peripheral blood)
Primary hematopoietic stem cells (male)
Primary hematopoietic stem cells (female)
Primary B cells (peripheral blood)
Placenta
Pancreas
Ovary
NHLF (lung fibroblasts)
NHEK (epidermal keratinocytes)
NHDF (adult dermal fibroblasts)
HMEC (mammary epithelial)
Gastric
Foreskin melanocytes
Foreskin keratinocytes
Foreskin fibroblasts (2)
Foreskin fibroblasts (1)
Fetal thymus
Fetal stomach
Fetal muscle (trunk)
Fetal muscle (leg)
Fetal lung
Fetal kidney
Fetal small intestine
Fetal large instestine
Fetal heart
Fetal brain (male)
Fetal brain (female)
Fetal adrenal gland
vHMEC
GWAS haplotype
ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15022-4
Genomic assays SNP coverage Target based on CRISPRi/CRISPRa; or (iv) displaying allele-specific
(1) H3K27ac HiChIP reporter activity using MPRA (see the “Methods” section).
2776 (100%) Region
Hits from two strategies enrich for disease-associated SNPs.
(2) Chromatin accessibility Ideally, we would assess each assay by directly testing how well it
2776 (100%) Region enriches for causal variants among the full set of variants assayed.TF
However, using metrics like ‘precision’ and ‘recall’ would require
that the causal variants be known with certainty. Because they are
(3) Transcription factor binding not, we instead tested how well the methods enrich for variants in
TF 2776 (100%) Variant tight LD with the tag SNP (as these variants are in turn enriched
TF for true causal variants), calculating a ‘pseudo-precision’ and
A G
‘pseudo-recall’. For each assay, we therefore quantified (1) the
(4) CRISPRi 
174 (6.3%) number of tested SNPs considered ‘hits’ in the assay (nH), (2) the
(variants in/near accessible chromatin) Region number of tag SNPs for which at least one SNP in tight LD was
∼30 guides/element
dCas9-KRAB tested in the assay (nT; i.e. recoverable tag SNPs), and (3) the
number of tag SNPs for which at least one SNP in tight LD was
(5) CRISPRa
3 considered an assay hit (nTH; i.e. recovered tag SNPs) (Supple-P
AI 2501 (90.1%)
NF ∼5 guides/element RegionT mentary Fig. 9a). We next calculated the pseudo-precision and
dCas9-VP64 pseudo-recall for GWAS variants for each assay. Here, we define
‘pseudo-precision’ as nTH/nH, representing the fraction of all
(6 and 7) MPRA (lentiviral and transfection) SNPs considered hits that are recovered tag SNPs, and ‘pseudo-
Reference allele 2695 (97.1%)
A
∼250 barcodes/element Variant recall’ as nTH/nT, representing the fraction of tag SNPs that
Alternate allele
G are recovered by being in tight LD with one or more hits. These
terms are similar to precision and recall except that a single causal
Fig. 2 Seven approaches for characterizing non-coding genetic variants. SNP might underlie multiple tag SNPs (by being in tight LD to
Genomic assays (left), the coverage of all common genetic variants in the each of them), making a pseudo-precision above 1 possible. By
605 kb locus (middle), and whether the assay is specific to genomic regions these measures, a highly effective approach would recover all tag
or variants (right), grouped into observational (top) and perturbational SNPs (pseudo-recall= 1) with as few SNP hits as possible (high
assays (bottom). (1) HiChIP can be used to identify active chromatin pseudo-precision). In the calculation of pseudo-precision and
regions (H3K27ac labeled) that interact with the TNFAIP3 promoter. (2) pseudo-recall, we did not consider GWAS tag SNPs that had no
DHS and ATAC-seq can be used to identify regions of accessible assayed variants in tight LD with that tag SNP (including the tag
chromatin. (3) Variants predicted to alter TF binding can be identified using SNP itself) in order to not falsely penalize the assays for technical
motif analysis in combination with evidence of TF binding by ChIP-seq42. failures (e.g., lack of PAM site for CRISPR or poor coverage in
Also see Supplementary Data 3–9. (4 and 5) Pooled CRISPRi and CRISPRa MPRA). We conducted these analyses for all variants and for the
screens can determine regulatory potential of each region by repressing subset of variants that lie in accessible chromatin in one of the
(CRISPRi) or artificially inducing (CRISPRa) each targeted region. (6 and 7) three blood cell types studied (because GWAS variants are enri-
MPRA (with lentiviral or transfection delivery strategies) can be used to ched in accessible chromatin1,4 and accessibility data is readily
test for allele-specific reporter expression. available for many cell types) (Fig. 3a, b).
To determine whether the pseudo-precision/pseudo-recall
Supplementary Data 4–6). We also applied CRISPRa (which relies performance of each method is better than expected by chance,
on dCas9-VP64 with MS2 stem loops that recruit HSF1 and p65 we created an empirical null distribution by randomly permuting
to artificially activate gene expression14), using guides that target the hit status among the assayed SNPs (1000 permutations) or by
50 bp regions surrounding each variant in the TNFAIP3 locus to shifting the hit status of each SNP to the next adjacent assayed
identify regions with the potential to induce TNFAIP3 expression SNP (Supplementary Fig. 9b, c). The shift approach preserves
(Supplementary Fig. 7, Supplementary Data 5–7). For shared positional clustering of hits inherent to LD and to some of the
guides and regions, we confirmed that CRISPRi and CRISPRa assays (e.g. CRISPRi, open chromatin). This reduces inflation of
drove the expected opposing changes in expression of TNFAIP3 positive hits within the null that may occur by permutation,
(Supplementary Fig. 7d, e). We also tested for allele-specific where the permuted hits may be in LD with many more tag SNPs
reporter expression induced by individual variants using MPRAs. than are possible given the clustered nature of the assay (thus
We synthesized all alleles for each variant, centered in 150 bp of increasing pseudo-precision and pseudo-recalls) (see the “Meth-
the surrounding reference DNA. These were cloned upstream of ods” section). Both shifting and permutation yielded similar
the TNFAIP3 promoter driving the expression of a GFP gene that results for SNPs in tight LD with GWAS tag SNPs (Fig. 3c,
contained sequence barcodes in the 3′ UTR. We used these Supplementary Fig. 9d). For each method, we compared the
barcodes to read out expression of each allele by RNA-seq. We pseudo-precision and pseudo-recall of actual data to the null
delivered them to immune cell lines by either lentivirus distribution. We did this both for all variants (Fig. 3c) and the
(L-MPRA) to integrate them into chromosomes, or transfection variants located in accessible chromatin in the three blood cell
(T-MPRA) to generate extrachromosomal reporters (Fig. 2; types (Fig. 3d).
Supplementary Fig. 8, Supplementary Data 8, 9). Variant-driven Relative to all variants, most of the methods (ATAC-seq on our
expression of the reporter was reproducible within, but not cell lines, Blood DHS+ATAC-seq on our cell lines, TF ChIP+
between, the two delivery methods (Supplementary Fig. 8b). motif, L- and T-MPRA, and CRISPRa) did not show a significant
For each assay, we determined which SNPs scored as ‘hits’ enrichment for GWAS variants (Fig. 3c). However, CRISPRi
based on SNPs being within regions annotated as: (i) interacting showed 7.5-fold enrichment for GWAS variants (95% C.I.,
with the TNFAIP3 promoter by HiChIP; (ii) accessible by ATAC- [0.9375; ∞]), albeit not significant (P= 0.087, empirical P-value
seq/DHS; (iii) within a region that modulates TNFAIP3 expression with genomic-shifts null) (Fig. 3a, c, Supplementary Fig. 9a).
4 NATURE COMMUNICATIONS |         (2020) 11:1237 | https://doi.org/10.1038/s41467-020-15022-4 | www.nature.com/naturecommunications
Perturbational Observational
TNFAIP3
NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15022-4 ARTICLE
a All variants c All variants
n CRISPRi; p = 0.087H 132 365 99 24 62 145 466
95% CI
Null median
0.6 Actual data CRISPRa; p = 0.407
nT 34 34 34 22 34 34 34
TF ChIP + motif; p = 0.504
nTH 17 29 18 15 14 23 32
0.4 ATAC cell lines; p = 0.629
e C if in
l li TA o
t R Ra RA RA
l P P P P L-MPRA; p = 0.376
 ce + 
A  m S
 P 
+ RI IS -M -M
AC S hI
C CR L T
H 0.2 Blood DHS + ATAC; p = 0.345
AT D  CTF T-MPRA; p = 0.137
0.0
0.00 0.25 0.50 0.75 1.00
Pseudo-recall
b Variants in blood accessible chromatin d Variants in blood accessible chromatin
1.5
nH 69 22 49 14 71
L-MPRA; p = 0.128
nT 29 18 29 29 29
1.0 CRISPRi; p  = 0.215nTH 16 13 12 15 26
CRISPRa; p  = 0.746
if i
ot R Ra A A
m R RSP P P P TF ChIP + motif; p  = 0.838+ I S
IP
 CR CR
I M
L- T-
M 0.5
h T-MPRA; p  = 0.011
F 
C
T
0.0
0.00 0.25 0.50 0.75
Pseudo-recall
Fig. 3 Comparison of GWAS enrichment across methods. a, b Values for nH, nT, and nTH for all methods, considering (a) all variants, and (b) only variants
in open chromatin. c, d Pseudo-precision (y axes) and pseudo-recall (x axes) for GWAS enrichment for each assay (colors), with diamonds depicting the
actual assay performance (as in Supplementary Fig. 9a), and the lines depicting the 95% CI of each assay’s null distribution (as in Supplementary Fig. 9b).
Empirical one-sided P-values derived from the genomic-shifts null are indicated next to each assay label. P-values are not corrected for multiple hypothesis
testing. c Each assay evaluated individually for all tested variants and d considering only SNPs in blood cell accessible chromatin. The relationship between
pseudo-precision and pseudo-recall is linear in the null (pseudo-precision= (nT/nH) × pseudo-recall) because both are proportional to nTH and nT and nH
are constant.
After restricting our analysis to variants located in accessible Data 10), determining, in this case, the number of credible sets (nT')
chromatin in the three blood cell types, several of the methods that were recovered (nTH') by containing one or more assay hits
(CRISPRa and TF ChIP+motif) again showed no significant (nH). Although the SNPs in a credible set are more likely to be
enrichment for GWAS variants. However, T-MPRA showed causal than when doing LD expansion, the limited availability of
significant enrichment (P= 0.011, empirical P-value with fine-mapping data restricted this analysis and reduced our
genomic-shifts null; 1.44-fold enrichment for GWAS, 95% CI statistical power. We calculated the pseudo-precision and pseudo-
[1.04; 5.2]; Fig. 3d, Supplementary Fig. 9e). recall for GWAS variants for each assay in an analogous way
Both L-MPRA and T-MPRA showed greatly increased pseudo- (Supplementary Fig. 9f–k). The rates from the credible set-based
precision with only marginally reduced pseudo-recall when analysis generally showed similar trends to the tag SNP approach,
restricting attention only to variants in accessible chromatin but were less significant due to the reduced sample size (Fig. 3c, d
(Fig. 3d, Supplementary Figs. 9e and 10). This may be because vs. Supplementary Fig. 9d–k); in addition, pseudo-precision was
many variants have the capacity to alter expression when tested in necessarily reduced for fine mapping due to reduced number of
an enhancer assay (such as MPRA), but do not reside in a region association signals, but with no change in assays hits.
of accessible chromatin in the relevant cell types and thus do not
alter disease risk. Although L-MPRA performed well for variants
in accessible chromatin, having the highest pseudo-precision of Prioritization of variants in disease-associated haplotypes.
any assay, there was limited power to evaluate L-MPRA because Finally, we used our analysis of genomic assays to prioritize SNPs
only four variants (in tight LD to 15 tag SNPs) out of the 19 on each disease-associated haplotype (Fig. 4, Supplementary
L-MPRA hits were in accessible chromatin (P= 0.128, empirical Data 3). We annotated as high-priority those variants that were
P-value with genomic-shifts null; Fig. 3d). hits in at least one of the two assays with the best performance
For CRISPRi, pseudo-precision and pseudo-recalls changed (CRISPRi for all variants and T-MPRA variants in accessible
little when focusing only on variants in accessible chromatin chromatin), finding a total of 18 such high-priority variants (Fig. 4,
(Fig. 3a–d, Supplementary Fig. 10), but pseudo-precision was asterisks). Of the 15 disease-associated haplotypes, nine included
less significant (P= 0.215, empirical P-value with genomic-shifts one or more of these 18 SNPs. These included five SNPs that had
null) because some of the SNPs tested lay just outside (within been fine-mapped in the UK Biobank, lying in 95% credible sets
100 bp) regions of accessible chromatin (Fig. 3c, d, Supplemen- representing associations with allergy, all autoimmune diseases
tary Fig. 9d, e). combined, and eosinophil counts (Fig. 4, Table 1).
We also considered another alternative proxy for causal variants, Several of these high-priority variants had other evidence
using credible sets from fine-mapping studies (Supplementary supporting a role in disease. For example, rs6927172 is the only
NATURE COMMUNICATIONS |    (     2020) 11:1237 | https://doi.org/10.1038/s41467-020-15022-4 | www.nature.com/naturecommunications 5
Pseudo-precision Pseudo-precision
ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15022-4
Associated disease: MS, IgA deficiency, RA SLE, RA, Sjogrens, IgA deficiency, ISD, Allergy ISD, SLE, RA Blood metabolites Urine metabolites
Haplotype: 1 2 3 4 5
Total assay hits
4 ** * * *
3 ** * *
2 * * * * *
1
0
Assay hits
Accessible chromatin
TF ChIP + motif
HiChIP
T-MPRA
L-MPRA
CRISPRi
CRISPRa
Fine-mapping
UKBB traits PIP10
UKBB traits CS
Haplotype composition
European (CEU)
Asian (CHB/JPT)
RA, IBD, celiac,
SLE, asthma,
IgA deficiency,
Sjogrens, ISD,
T1D, primary RA,
Associated disease: Biliary cirrhosis RA RA ISD, psoriasis SLE Psoriasis Celiac MS Celiac MS
Haplotype: 6 7 8 9 10 11 12 13 14 15
Total assay hits
4 *
3 * *
2 * * *
1
0
Assay hits
Accessible chromatin
TF ChIP + motif
HiChIP
T-MPRA
L-MPRA
CRISPRi
CRISPRa
Fine-mapping
UKBB traits PIP10
UKBB traits CS
Haplotype composition
CEU
CHB/JPT
Fig. 4 Prioritizing variants on disease-associated haplotypes. A summary of assay results and genetics data for all SNPs on each disease-associated
haplotype. Each subpanel represents a different haplotype, with associated traits and the haplotype number are indicated on the top. For each SNP (x axes),
the total number of assay hits is shown in the bar graph (top) with SNPs that are hits in CRISPRi or T-MPRA hits in accessible chromatin marked with an
asterisk. Results from each assay are shown in the middle, with hits in red, and SNPs that are assayed but were not hits in gray for each of the seven assays
(y-axis). The vertical black bars above accessible chromatin SNP status indicate SNPs that were in accessible chromatin in our tested cell lines. Fine-
mapped immune-related traits from UK Biobank (UKBB), including SNPs in the 95% credible set (CS—blue) and those that have a posterior inclusion
probability > 10% (PIP10—green) are second from the bottom. The population-specific SNPs contained within each disease-associated haplotype are
indicated (bottom) with orange for European (CEU) and purple for East Asian (CHB/JPT). Also see Supplementary Data 3.
high-priority variant on haplotype 6 (which lay in accessible haplotype had evidence of impact (with rs111710107 only being
chromatin and scored in the T-MPRA assay, but not in the in accessible chromatin, and rs111231590 having allele-specific
CRISPRi assay); this variant is associated with many diseases, reporter activity according to both T- and L-MPRA assays).
including RA, SLE, celiac, T1D, and asthma, and it is a fine- Similarly, rs643177 is one of two high-priority variants on
mapped SNP in our analysis of combined autoimmune disease in haplotype 9 (laying in accessible chromatin and a hit in T-MPRA
the UK Biobank and in previously reported studies of ulcerative assay, but not tested in CRISPRi due to the lack of a suitable
colitis, RA, and celiac)1 (posterior inclusion probability (PIP)= guide-RNAs). This variant also had evidence of interaction with
0.1343; Table 1, Supplementary Data 2). This variant also has the TNFAIP3 promoter according to HiChIP, and had allele-
evidence of allele-specific ATAC-seq and allele-specific ChIP-seq specific reporter expression in L-MPRAs. rs643177 is a fine-
for the TFs NF-kB and JunD in lymphoblastoid cell lines44,45 and mapped psoriasis SNP1 and has evidence of allele-specific binding
allele-specific ATAC-seq and allele-specific ChIP-seq for the NF- of the TF Pou2f1 (Table 1). The other high-priority variant on
kB1 p50 subunit in primary CD4 T cells20. It appears to interact haplotype 9 is rs559766217, which was a hit in the CRISPRi assay,
with the TNFAIP3 promoter by 3C, has allele-specific reporter is in accessible chromatin and contacts the TNFAIP3 promoter
activity according to a luciferase assay, and lays in a region that according to HiChIP. Four of the 17 other SNPs on the haplotype
affected TNFAIP3 expression based on 11–12 bp CRISPR- have some evidence of impact (including rs538522 and rs598493,
induced deletions46,47. Only two of the other 10 variants on the which interact with the TNFAIP3 promoter according to HiChIP;
6 NATURE COMMUNICATIONS |       (  2020) 11:1237 | https://doi.org/10.1038/s41467-020-15022-4 | www.nature.com/naturecommunications
# Assay hits # Assay hits
rs6933404 rs17066662
rs62432712 rs12524433
rs2327832 rs151268394
rs928722 rs35706049
rs35926684 rs13207033
rs6927172 rs4896289
rs11757201 rs34810205
rs17264332 rs13192841
rs6920220 rs12527282
rs111710107 rs2056224
rs111231590 rs2056223
rs35858108
rs7758499 rs13199053
rs9399221 rs66469519
rs2876370 rs73558137
rs9376287 rs68161003
rs36123502
rs62434173 rs67178054
rs62434174 rs67355614
rs62434176 rs4896295
rs3900757 rs34996815
rs62434193 rs617328rs10499194
rs80351603
rs654912 rs13205649
rs583522 rs66499821
rs598493 rs34654849
rs643177 rs12525643
rs582757 rs34324947
rs610604
rs642627 rs142373084
rs622091 rs58905141
rs1082428 rs10499197
rs1082425 rs569305282
rs612217 rs77990078
rs11356730 rs79411652
rs674451 rs78515370
rs67543742 rs112020444
rs559766217 rs80213143
rs4895498 rs111425486
rs6933987 rs113496608
rs6909442 rs112975146
rs112497003
rs61117627 rs111883038
rs58721818 rs80126770
rs60022938
rs9494883
rs654874 rs80300819
rs2307860 rs9494885
rs629953 rs7753873
rs2307859 rs7767264
rs661561 rs7774101
rs610604 rs9494886
rs622091 rs59699063
rs1082428 rs61593413
rs603904 rs59693083
rs1082425 rs5029924
rs601035 rs5029926
rs11356730 rs5029928
rs593622 rs3757173
rs592810 rs57087937
rs644340 rs7750604
rs674451 rs5029930
rs11354287 rs719149
rs9321634 rs719150
rs4896301 rs5029937
rs1561121 rs5029939
rs6570193 rs2230926
rs6903624 rs5029949
rs4896303 rs3834310
rs9494892
rs9376293 rs9494893
rs542829 rs7752903
rs12665429 rs9494894
rs606539 rs148314165
rs606103 rs200820567
rs636393 rs7749323
rs655112 rs9494895
rs661675 rs77000060
rs602414 rs6932056
rs601705 rs6570194
rs694069
rs508214 rs9494885
rs7753873
rs67297943 rs7767264
rs11970411
rs75783118 rs7774101
rs1878658 rs9494886
rs1002658 rs59699063
rs77027760 rs61593413
rs142761146 rs59693083rs5029926
rs5029928
rs9321623 rs3757173
rs9402908 rs57087937
rs9321624 rs7750604
rs611460 rs5029930
rs9321625 rs719149
rs9373197 rs719150
rs9389526 rs6570194
rs600469
rs533091 rs76227521
rs9373198 rs9494838
rs9484077 rs77080039
rs1547295 rs79731087
rs7769192 rs376951159
rs76854395 rs12213095
rs9402909 rs371949948
rs9389527 rs12198885
rs9484078 rs60046935
rs9385802 rs60781459
rs535730 rs56789590
rs9389529 rs76822624
rs9484079 rs139226197
rs1889135 rs17066649
rs928721 rs12197157
rs928720 rs12210705
rs486165 rs12199435
rs9389530 rs12201110
rs609905 rs12201389
rs688096 rs12201430
rs593115 rs12192746
rs583719 rs67664821
rs583816 rs9688558
rs2090108 rs9688555
rs522563 rs12201688
rs9376290 rs12194668
rs9376291
rs505662
rs509260 rs72989420
rs540414 rs4896275
rs686851 rs4895491
rs567915 rs72989424
rs9321626 rs72989427
rs627031 rs72989428
rs11758213 rs72989430
rs493179 rs72989434
rs678645 rs4896276
rs693699 rs4896277
rs598047 rs34979694
rs525977 rs60073400
rs529693 rs59086769
rs6904167 rs55695001
rs10782265 rs12110715
rs628791 rs12111342
rs7768363 rs72989455
rs530181 rs72989461
rs2137868 rs724371
rs481095 rs12660587
rs9402913 rs72989470
rs678385 rs12661016
rs675520 rs72989479
rs34885134 rs12664666
rs2152655 rs10484593
rs590861 rs12661817
rs605777 rs10484592
rs60221632 rs72971414
rs652401 rs17066538
rs665668
rs666584
rs666619
rs2327837
rs10607561
rs6921233
rs1518048
NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15022-4 ARTICLE
Table 1 Disease-associated variants positive for CRISPRi or chromatin accessibility with T-MPRA.
SNP ID Associated trait Tehranchi Tehranchi Fine mapped UKBB Other evidence Haplotype Hit in assays
asATAC asChIP 95% CS (SuSiE)
rs200820567 Allergy, ISD, RA, SLE, eosinophil x Eosinophil counts (PIP= Fine mapped in Adrianto et al. 2 T-MPRA+ accessible
counts, IgA deficiency, Sjogren’s 0.03); Allergy (PIP= (SLE) chromatin
0.04)
rs148314165 Allergy, ISD, RA, SLE, eosinophil Fine mapped in Adrianto et al. 2 T-MPRA+ accessible
counts, IgA deficiency, Sjogren’s (SLE) chromatin
rs112497003 Allergy, ISD, RA, SLE, eosinophil Eosinophil counts (PIP= 2 T-MPRA+ accessible
counts, IgA deficiency, Sjogren’s 0.01) chromatin
rs111883038 Allergy, ISD, RA, SLE, eosinophil Eosinophil counts (PIP= 2 L-MPRA+ accessible
counts, IgA deficiency, Sjogren’s 0.01) chromatin; T-MPRA
+ accessible
chromatin
rs6927172 Celiac, IBD, RA, Asthma, IgA x NF-kB, Combined Autoimmune Fine mapped in Huang et al. 6 T-MPRA+ accessible
deficiency, Sjogren’s, ISD, T1D, JunD (PIP= 0.13) (UC, PIP= 0.06); Farh et al. chromatin
primary biliary cirrhosis (RA, PIP= 0.11; Celiac, PIP=
0.19; UC, PIP= 0.23);
Westra et al. (RA, PIP= 0.10)
rs643177 ISD, psoriasis Pou2f1 Fine mapped in Farh et al. 9 L-MPRA+ accessible
(Psoriasis, PIP= 0.15) chromatin; T-MPRA
+ accessible
chromatin
rs59086769 Urine metabolites 5 T-MPRA+ accessible
chromatin
rs1002658 Celiac x NF-kB, 14 T-MPRA+ accessible
PU.1 chromatin
rs11758213 MS x JunD Fine mapped in Huang et al. 15 T-MPRA+ accessible
(UC, PIP= 0.075) chromatin
rs9389527 MS 15 T-MPRA+ accessible
chromatin
rs12201430 Blood metabolites x 4 T-MPRA+ accessible
chromatin
rs12192746 Blood metabolites 4 L-MPRA+ accessible
chromatin; T-MPRA
+ accessible
chromatin
rs34654849 MS, IgA deficiency, RA 1 T-MPRA+ accessible
chromatin
rs73558137 MS, IgA deficiency, RA 1 T-MPRA+ accessible
chromatin
rs5029924 Allergy, ISD, RA, SLE, eosinophil BJAB asATAC and fine 2 CRISPRi
counts, IgA deficiency, Sjogren’s mapped in Farh et al. (SLE,
PIP= 0.09)
rs5029926 Allergy, ISD, RA, SLE, eosinophil 2, 3 CRISPRi
counts, IgA deficiency, Sjogren’s
rs10499197 Allergy, ISD, RA, SLE, eosinophil 2 T-MPRA+ accessible
counts, IgA deficiency, Sjogren’s chromatin; CRISPRi
rs58905141 Allergy, ISD, RA, SLE, Eosinophil Eosinophil counts (PIP= 2 L-MPRA+ accessible
Counts, IgA deficiency, Sjogren’s 0.02) chromatin; CRISPRi
rs559766217 ISD, Psoriasis 9 CRISPRi
Variants that are positive for either chromatin accessibility with T-MPRA or CRISPRi are listed with their associated trait, and whether they were also positive in Tehranchi et al. as having allele-specific
ATAC (asATAC) or asChIP-seq for TFs in LCLs. Our fine-mapping data using UKBB traits for the 95% credible set variants are included, and other fine-mapping data or evidence for SNP functionality is
listed in Other Evidence. The haplotype for the SNP is listed in Haplotype.
rs598493 and rs610604, which are located in accessible chroma- disease-associated variants—provided that the relevant cell types
tin; and rs6909442, which has allele-specific reporter expression are known and can be studied (which remains a serious
according the T-MPRA assay). limitation).
Other examples include rs11758213 on haplotype 15, which is To study the potential utility of various genomic features for
in the 95% credible set for ulcerative colitis (PIP= 0.0074)48 and prioritizing non-coding variants, we studied seven genomic assays
had evidence of allele-specific ATAC-seq and ChIP-seq for the TF in three disease-relevant cell types to assess to the extent to which
JunD in LCLs44,45 and rs1002658 on haplotype 14, which was they enrich for disease-associated variants within a set of 2776
associated with celiac disease and had evidence of allele-specific common non-coding SNPs in the TNFAIP3 locus. We found
ATAC-seq and ChIP-seq for the TFs NF-kB and PU.144,45. significant enrichment among high-scoring SNPs for two meth-
Interestingly, haplotype 2 had 41 of 51 SNPs that scored as hits in ods: (1) variants present in CRISPRi-responsive regulatory
at least one of the seven assays, including five SNPs in accessible regions and (2) variants present in accessible chromatin that also
chromatin that score as hits in the T-MPRA assay and three SNPs showed allele-specific reporter activity by T-MPRA. These two
that scored as hits in the CRISPRi assay (Table 1). criteria identified 18 TNFAIP3 variants associated with 15 dis-
eases on 9 haplotypes; potential functional roles for these variants
in immunity were supported by additional published data (such
Discussion as allele-specific ATAC-seq, ChIP-seq, and genetic fine-map-
GWASs effectively narrow down the search for causal variants to ping). By contrast, the other genomic features did not provide
a small set of candidates, but determining which of the candidates significant enrichment.
contributes to disease risk remains a challenge. Because disease- Our data support two prioritization schemes (CRISPRi and
causal variants are likely to be correlated with functionally rele- accessible chromatin with T-MPRA) as viable methods for
vant genomic features in the cell types in which they act, it should enriching for causal variants in the TNFAIP3 locus. However,
be possible to use genomic features to help inform the search for since perturbational methods (e.g. CRISPRi, MPRA) cannot
NATURE COMMUNICATIONS |         (2020) 11:1237 | https://doi.org/10.1038/s41467-020-15022-4 | www.nature.com/naturecommunications 7
ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15022-4
currently be scaled to the same level as observational methods 3 of the 1000 genomes project (http://www.internationalgenome.org/)). We used
(ATAC-seq, ChIP-seq, HiChIP, and TF motif analysis), we could tabix (0.2.5) (tabix -h ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/
not test the generalizability of our ndings to additional disease- supporting/vcf_with_sample_level_annotation/ALL.chr6.phase3_-fi
shapeit2_mvncall_integrated_v5_extra_anno.20130502.genotypes.vcf.gz
associated loci, variants, and cell types. 6:137846078-138453052 > TNFAIP3.vcf)49 and vcftools (0.1.15) and Plink
Our focus on the TNFAIP3 locus helps to clarify a complex (v1.90b3d; vcftools --vcfg TNFAIP3.vcf --keep CEU_names.txt --out CEU --plink;
region with many genetic associations through analysis of variant plink --file CEU --out CEU; plink --bfile CEU --maf 0.01 --geno 0.01 --hwe 0.01
50,51
features and functions in three main immune cell types. While our --out CEU.filtered --make-bed) to extract all alleles in the TNFAIP3 locus
(chr6:137846078-138453052, hg19) that were MAF ≥ 0.01 from the 1000 genomes
data corroborated two reported putatively causal variants asso- phase 3 database CEU and the combination of CHB and JPT populations. For trait-
ciated with lupus (rs200820567 and rs148314165 on haplotype 2), associated variants, we reanalyzed GWAS summary statistics (www.immunobase.
it highlighted six other variants (rs58905141, rs10499197, org and refs. 52,53) for tag SNPs and those in tight LD in 1KG samples, according to
2 54
rs5029924, rs5029926, rs112497003, and rs111883038) on the same the same population in which the study was conducted (r > 0.8; 294 SNPs) .
haplotype that are also putatively causal. Whether these variants act GWAS haplotypes were defined on the basis of tight LD between GWAS tagSNPs and other genomic SNPs. We calculated LD between GWAS tag SNPs and
in concert to confer risk at this haplotype needs to be examined. other SNPs using Plink (v1.90b3d; --r2 inter-chr --ld-window-r2 0.2) for both East
While we found prioritized variants for nine haplotypes, none were Asian (EAS) and European (CEU) populations using 1000 Genomes data54. Each
found for another six haplotypes, which could be explained by lack GWAS tag SNP and all SNPs in tight LD (r2 > 0.8) within the GWAS population of
of assay sensitivity or the variants being biologically active in other study defined our initial haplotype estimates. Any of these haplotypes that sharedSNPs in tight LD (r2 > 0.8) were then merged into a single haplotype until none
cell types and conditions. Interestingly, many haplotypes contained showed any overlap, yielding 15 haplotypes associated with one or more diseases.
associations to different diseases, affirming that different auto- We found that the number of haplotypes identified was robust to this cutoff
immune diseases could have similar autoimmune genetic etiology between 0.76 and 0.89. Haplotypes identified in EAS and CEU that had any overlap
because they are presumably promoting disease through the same between the GWAS Tag SNPs were merged into a single haplotype, with
1 population-specific membership indicated in Fig. 1 and Supplementary Fig. 2, andcausal genetic variants . Our data help explain the immense genetic Supplementary Data 1. We used phased 1000 Genomes genotypes to ensure each
complexity of the locus by prioritizing 18 of the 293 disease- haplotype exists at >0.5% in each population.
associated variants, although there may be even more disease-
causal variants to be found in different contexts. Genetic fine-mapping. We performed genetic association and fine-mapping in up
Our study of common variants in the TNFAIP3 locus provides to 361,194 unrelated, white British individuals from the UK Biobank55, as deter-
a strategy to help guide future variant characterization studies at mined by the PCA-based sample selection criteria (https://github.com/Nealelab/
other loci. Increasingly accurate approaches to identify causal UK_Biobank_GWAS/blob/master/ukb31063_eur_selection.R). We restricted to allimputed variants with MAF > 0.01% (except for missense and protein-truncating
variants will require the development and integrated analysis of variants annotated by VEP56, MAF > 0.0001%), Hardy–Weinberg equilibrium P-
experimental methods that assess variant function. value > 1 × 10-10, and imputation quality (INFO) > 0.8 (https://github.com/
Nealelab/UK_Biobank_GWAS). To perform association tests for binary pheno-
types, we used a generalized linear-mixed model as implemented in SAIGE57
Methods v0.29.4 with the minimum minor allele count (MAC) threshold, MAC > 10 for
Cell culture and stimulation of immune cells. BJAB (DSMZ, cat. no. ACC 757), each GWAS. To perform association tests for quantitative phenotypes, we used a
Jurkat, Clone E6-1 (ATCC, cat. no. TIB-152), U937 (ATCC, cat. no. CRL-1593.2), linear-mixed model as implemented in BOLT-LMM58 v2.3.2 with default settings.
THP-1 (ATCC, cat. no. TIB-202), and GM12878 (Coriell, cat. no. GM12878 LCL Phenotypes for combined autoimmune disease were derived as previously
from B-Lymphocyte) cell lines were cultured using RPMI 1640 (ThermoFisher, defined58, allergy status was self-reported, and eosinophil counts were inverse rank-
21870092) containing 10% fetal bovine serum (FBS, VWR, 97068-091; 20% for based normal transformed. We included sex, age, age2, sex × age, sex × age2, and
GM12878) with 1% Penn/strep (VWR, 45000-652), 1% L-glutamine (Thermo- top 20 principal components as covariates. Genetic fine-mapping was performed
Fisher, 25030081), and 1% HEPES (Sigma, H0887-100ML). Cells were maintained using FINEMAP v1.359,60 and the summary statistics version of susieR43 v0.7.1
at a culture density between 100K and 1M cells/mL. Jurkat T cells were stimulated with the maximum number of causal variants specified as 10. LD matrices were
with 2.5 μg/mL of anti-CD3 (Biolegend, 317304) and 10 ng/mL of PMA (Sigma, calculated from imputed dosages for individuals included in each GWAS using
P1585-1MG) for 1 h prior to harvesting for CRISPRi and MPRA, and 1 and 4 h for LDstore61 v2.0b. Individual variant posterior inclusion probabilities and condi-
ATAC-seq experiments. BJAB and GM12878 B cells were stimulated with tional 95% credible sets are reported.
2.5 μg/mL of anti-IgM (Sigma-Aldrich, 86620270) and 2 μg/mL anti-CD40
(ThermoFisher, 14-0409-82) for 2 h for CRISPRi and MPRA, and 1 and 4 h for
ATAC-seq and 4C (BJAB) experiments. THP-1 and U937 monocytes were sti- GWAS immune cell enrichments. Heritability enrichments of traits (Fig. 1;
mulated with 100 ng/mL LPS (Invivogen, tlrl-peklps) for 2 h for CRISPRi and Supplementary Fig. 3) in cell lines and cell types were estimated using stratified
MPRA, and 1 and 4 h for ATAC-seq and 4C (U937) experiments. LD-score regression (s-LDSC) over accessible chromatin or histone modifications
in specific cell types as previously reported35 by interpreting the cell type-specific
repression coefficient in s-LDSC model. For hematopoietic cell types and cell lines,
Lentivirus preparation. HEK293T cells were grown using DMEM (VWR, 45000- common variants overlapping accessibility peaks from ATAC-seq data for 13
316) with 10% FBS, 1% Penn/Strep, 1% L-glutamine, 1% HEPES (10DMEM). Cells primary cell types62 were used to compute the heritability enrichment. For broad
were passaged at 80% confluence for each passage. To make lentivirus, media was tissue enrichments, DNase Hypersensitivity peaks and H3K27ac and H3K4me1
aspirated from the adherent cells and Trypsin EDTA 0.25% (VWR, 45000-664) was ChIP-seq peaks were overlapped with common variants to compute heritability
used to create a single-cell suspension; the cells were kept at 37 °C for 4 min with enrichments. The −log10 P-values for the s-LDSC regression terms for each spe-
Trypsin, and 10DMEM was added to a final concentration of 80%. The cells were cific annotation were shown as a measure of enrichment.
pipetted up and down until they were in a single cell suspension. They were then
counted and plated in a six-well plate at 500K cells/well in 2 mL 10DMEM. The
next day, when the cells were ~70% confluent, they were transfected. pVSV-G HiChIP data and analysis. H3K27ac HiChIP data previously generated41 were
(0.1 μg; Addgene, 8454), pPAX2 (1 μg; Addgene, 12260), and the donor plasmid downloaded in.fastq format from GEO accession “GSE101498”. Biological and
(1 μg), were added to 125 μL of OPTI-MEM and mixed. 6 μL of the TransIT-LT1 technical replicates of Th17, Naïve T-cell, and GM12878 H3K27ac samples were
(Mirus Bio, MIR 2300) transfection reagent was added to a separate tube of 125 μL pooled and aligned with Hi-C Pro63. Virtual 4C plots (Supplementary Fig. 5) using
OPTI-MEM (ThermoFisher, 31985062) and mixed. The OPTI-MEM LT-1 mixture a resolution of 2.5 kb and a rolling mean of 2.5 windows41. Per-fragment estimates
was then added to the OPTI-MEM plasmid mixture, mixed, and incubated at RT of interaction strength to the TNFAIP3 promoter were generated using hichipper64
for 15 min. The mixture was then added dropwise to the well. The plate was then and normalizing to the total number of unique fragments in each library. We used
swirled in order to ensure distribution of the mixture and effective transfection. a normalized interaction score of 20 to annotate regions as TNFAIP3 interacting.
The cells were put at 37 °C to incubate overnight and the media was changed at
24 h post-transfection, this time using 10DMEM with 1% BSA (Sigma, A7979). The
ATAC-seq. We used the FAST-ATAC protocol62. 10,000–20,000 cells were sorted
cells were then incubated at 37 °C for 16 h, and the supernatant was harvested. The
into RPMI 1640 containing 10% fetal bovine serum. The cells were centrifuged at
viral supernatant was spun at 500 × g for 5 min to separate cellular debris, and
500 × g for 5 min at 4 °C. All of the supernatant was aspirated, ensuring that the
stored at 4 °C for up to 3 months.
pellet was not disturbed. The pellet was then resuspended in the tagmentation
reaction mix (25 μL 2X TD Buffer (Illumina, 15027866), 2.5 μL TD Enzyme
1000 Genomes Project and GWAS catalog. We centered our study on the 2776 (Illumina, 15038061), 0.5 μL 1% Digitonin (Promega, G9441), 22 μL H2O) and
variants that lie within and 150 kb to either side of the TNFAIP3 TAD, yielding a mixed at 300 RPMs at 37 °C for 30 min on an Eppendorf Thermomixer. Imme-
605 kb locus (MAF > 0.01, combined CHB+JPT and CEU populations from phase diately after the incubation, samples were purified using a minElute kit (Qiagen,
8 NATURE COMMUNICATIONS |   (      2020)1 1:1237 | https://doi.org/10.1038/s41467-020-15022-4 | www.nature.com/naturecommunications
NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15022-4 ARTICLE
28006), eluting in 10 μL. The entire sample was PCRed (a 50 μL reaction with 25 μL U937, and THP-1 were made through lentiviral transduction of these cells with a
NEBNext, 2.5 μL F+R custom nextera primers (10 μM each; Supplementary doxycycline-inducible transactivator (ClonTech, 631363) and the TRE-dCas9-
Data 11), 10 μL of tagmented DNA, and 12.5 μL H2O) for five cycles with the KRAB-BFP construct (for CRISPRi; Addgene, 85449) or pMS2-p65-HSF1
following program (72 °C, 5 min; 98 °C, 30 s; five cycles of 98 °C, 15 s, 63 °C, 15 s, (Addgene, 73795) and dCas9-VP64-GFP (for CRISPRa; Addgene, 61422); for
72 °C, 1 min). We performed qPCR with 5 μL of the sample to determine the both, guide libraries were infected at an MOI < 0.3, and puromycin selected for
number of additional cycles required, while the rest remained on ice. The 5 μl of 4 days. Cells containing libraries were maintained in culture without doxycycline
sample was added to a qPCR mix (5 μL of PCR, 5 μl of NEBNext, 0.5 μL F+R and used for each replicate. For each replicate, cells were split and given doxycy-
custom nextera primers, 0.09 μL of 100X SYBR (Invitrogen, S7563), 4.41 μL H2O) cline 24–48 h prior to harvesting, and stimulated with relevant ligands 1–2 h prior
and qPCRed (98 °C, 30 s; 20 cycles of 98 °C, 15 s, 63 °C, 15 s, 72 °C, 1 min). The to harvesting.
number of cycles that it took to reach 1/3 the maximum fluorescence threshold in We performed FlowFISH screens69. For PrimeFlow experiments, 5 million cells
the qPCR was then applied via PCR to the original PCR sample. Libraries were were aliquoted in PBS in polypropylene tubes and centrifuged at 500 × g for 5 min.
cleaned using 1.5X Agencourt XP beads and ethanol washes per manufacturer’s All but 100 μL of the supernatant was discarded (this step is true for every
protocol. The DNA concentration of the sample was measured using Qubit and the centrifugation step in this protocol) and the cells were resuspended in the residual
average fragment size was determined using a TapeStation. Samples were then volume. Cells were then fixed according to manufacturer protocol (ThermoFisher,
multiplexed and sequenced using 50 bp paired end chemistry at an average read- 88-18005-210) using Fixation Buffer 1 for 30 min at 2–8 °C with rotating. Cells
count of 30M reads per sample. were then centrifuged at 800 × g for 5 min. and the supernatant was discarded. Cells
Paired-end ATAC-seq reads were mapped to the genome (hg19) using Bowtie2 were then permeabilized according to manufacturer protocol with addition of
(2.2.1; parameters: --maxins 2000), with duplicate reads removed using Picard RNase inhibitors through inversion, and centrifugation at 800 × g for 5 min, then
(2.20.6; MarkDuplicates REMOVE_DUPLICATES=true), and peaks (clusters of the supernatant was discarded. This step was repeated. A second fixation step was
reads representing open chromatin regions) called using Homer (4.6; findPeaks carried out using Fixation Buffer 2 according to manufacturer protocol, the
-style dnase). samples were mixed, and inverted for one hour in the dark at RT. The cells were
We calculated the ATAC-seq similarity between our cell lines and primary then centrifuged at 800 × g for 5 min at RT, and the samples were washed twice
immune cell types20 (Supplementary Fig. 4). We used pyatac (version 0.3.4) to get with PrimeFlow RNA Wash Buffer, centrifuging the samples at 800 × g between
read counts for each region previously identified as having been accessible in one or each wash for 5 min. The TNFAIP3 target probe (ThermoFisher, VA1-20723) was
more immune cell types, for GM1287865, Jurkat, BJAB, and U937. Pearson’s added at 1X in PrimeFlow RNA Target Probe Diluent, mixed thoroughly by
correlation coefficient was calculated comparing the log ATAC-seq counts (+0.5) pipetting up and down (100 μL of probe/sample), and incubated at 40 °C for 2 h,
per region to quantify the similarity between each of the primary immune cells as with inversion every 30 min. 1 mL of PrimeFlow RNA Wash Buffer was added to
well as the other cell lines, for each profiled cell line. These were sorted in each sample, the samples were inverted to mix, and centrifuged at 800 × g for
decreasing order and the top five for each cell line are displayed in Supplementary 5 min, and the supernatant was aspirated. Samples were then washed with 1 mL
Fig. 4. PrimeFlow RNA Wash Buffer containing RNase inhibitors twice followed by
centrifugation at 800 × g for 5 min. 100 μL of PrimeFlow RNA PreAmp Mix was
CRISPR screens. The guide libraries targeting the TNFAIP3 locus for CRISPRi then added to each sample and briefly vortex to mix, and the samples were then
and CRISPRa are available in Supplementary Data 4 and 7. To design the guide incubated for 1.5 h at 40 °C with mild vortexing once every 30 min. Samples were
library, all possible 20 bp sgRNAs with the Cas9 protospacer adjacent motif NGG washed three times with 1 mL of PrimeFlow RNA Wash Buffer, and they were
within the region surrounding TNFAIP3 (chr6:13784700–138453100, hg19) were centrifuged at 800 × g for 5 min, and the supernatant was aspirated. 100 μL of
considered. On-target scores for each guide were determined using the Rule Set 2 PrimeFlow RNA Amp Mix was then added to each sample, the samples were mixed
method described in ref. 66. To determine the number of off-target locations, by votexing, and were incubated for 1.5 h at 40 °C with mild vortexing once every
bowtie (0.12.7)67 was used to map guides to the human reference (hg19) with a 30 min. The cells were then washed twice in 1 mL of PrimeFlow RNA Wash Buffer
maximum 10,000 matches, with up to three mismatches (parameters: -n 3 -l 15 -e and centrifuged at 800 × g for 5 min. Each sample received 100 μL of PrimeFlow
10000 -y --all -S). Using this set of potential mapping locations in the genome, off- RNA Label Probe diluted in PrimeFlow RNA Label Probe Diluent and incubated
target score was calculated using the method of Hsu et al.68. Brie for 1 h at 40 °C with mild vortexing once at 30 min. Samples were then washed withfly, single off
targets were calculated as e moves over positional mismatches between guide and 1 mL of PrimeFlow RNA Wash Buffer at RT followed by centrifugation at 800 × g
off-target, where the m is as below and d is mean pairwise distance between for 5 min. The samples were then washed five times with 35 °C PrimeFlow RNA
mismatches: Wash Buffer following each wash with centrifugation at 800 × g for 5 min. Samples
were then left in 100 μL of PBS and stored in the dark at 4 °C until sorting.
ð  ½ Þ 
 1  1 Cells expressing CRISPRi or CRISPRa constructs along with sgRNA librariesΠ 1 W e ´
e2M ð19d
Þ
´
2
´ 4þ 1 n were sorted into six 10% bins, sorting on the extremes of expression (30% on eithermm19 the low or high portion of the expression distribution, each divided into three
contiguous bins each comprised of ~10% of the overall distribution). For each
M ¼ ½0; 0; 0:014; 0; 0; 0:395; 0:317; 0; 0:389; 0:079; 0:445; 0:508; 0:613; experiment and cell type, between 300K and 1M cells were sorted per bin. Genomic
0:851; 0:732; 0:828; 0:615; 0:804; 0:685; 0:583 DNA for each sample was then reverse-crosslinked using ChIP Lysis Buffer (1%
SDS, 0.01 M EDTA, 0.05 M Tris–HCl pH 7.5). Briefly, sorted cells were spun at
Individual off-targets are aggregated into a single guide using: 800 × g for 10 min at 4 °C, the supernatant was aspirated, and the cells were
100 resuspended in 50 μL of ChIP Lysis Buffer, and incubated at 65 °C for 10 min. The
Sguide ¼ þPn ð Þ samples were then cooled to 37 °C and 2 μL of RNase Cocktail (ThermoFisher,100 mmi¼1 Shit hi AM2286) was added to each sample and the sample was mixed well by pipetting,
On-target scores range from 0 to 100, with 100 being optimal. Off-target scores followed by incubation at 37 °C for 30 min. 10 μL of Proteinase K (NEB, P8107S)
range from 0 to 100 with 100 being no off-target effects predicted. CRISPRi guides was added to each sample and the sample was mixed well by pipetting, followed by
were selected to target the locations of ATAC-seq peaks from Jurkat, BJAB, or incubation at 37 °C for 2 h and then 95 °C for 20 min. gDNA was extracted using
U937—with or without stimulation (overlapping peaks merged), and aimed to tile Agencourt XP beads at 0.7X following the manufacturers protocol, and the sample
the region uniformly, with an average of ~30 guides per element. For CRISPRa, the was eluted at 100 μL. Libraries were prepared by PCR of each sample, splitting each
targeted elements were the locations of SNPs (±25 bp) and guides were selected to into four 50 μL reactions (25 μL NEB Next Master Mix, 2.5 μL barcoded sequencing
get ~5 guides per SNP; most SNPs with at least one guide (2501/2776). In both forward and reverse primers (Supplementary Data 11), 11.5 μL gDNA, and 11.5 μL
cases, we excluded guides for which there were off-target matches near the ddH2O; program: 98 °C for 30 s, 25 cycles of 98 °C for 15 s, 62 °C for 15 s, 72 °C for
TNFAIP3 locus, as well as any that had more than three off-target perfect matches 16 s, then 72 °C for 2 min. The libraries were then gel purified using a 2% gel
anywhere in the genome. We included 770 non-targeting guides in CRISPRa (expected band size of 206 bp). Samples were sequenced aiming to get >1,000,000
library and 6282 in the CRISPRi library, which were created by reversing (but not reads per bin, on either an Illumina HiSeq 2500 or MiSeq using a custom
complementing) selected targeting guides. Prior to synthesis, Gs were added to all sequencing and index primers for CRISPRi and CRISPRa (Supplementary
sgRNAs not starting with a G to aid in transcription efficiency. The sgOPTI vector Data 11).
(for CRISPRi; Addgene, 71409) or the sgSAM vector (for CRISPRa; made in house, For CRISPRi/a analysis, reads covering the guide sequences from each bin were
available upon request) was digested with BsmbI (NEB, R0580S) overnight, PCR aligned to the designed guide sequences using Bowtie2 (2.2.1; default settings)70,
cleaned, and the digest was repeated for two hours with thermostable alkaline and the total number of each guide observed in each bin counted. Read counts
phosphatase (Promega, M9910) added during the final hour of digestion. The cut from each bin were modeled as if originating from a negative binomial distribution,
vector was then gel purified using a 0.7% agarose gel. Guides for CRISPRi and where the underlying distribution of cells targeted by each guide had a log
CRISPRa (Supplementary Data 4, 7) were synthesized using Agilent Technologies (expression level) that was normally distributed for each guide, with the same
100K arrays, with common PCR priming sequences on each element. The oligos variance as the entire distribution (since most guides are expected to have no effect)
were amplified to add Gibson assembly homology arms, and inserted into the and different means (that varied based on the effect of the guide). The percent of
sgOPTI vector using Gibson assembly using 500 ng of vector and 70 ng of insert. cells that were sorted into each bin was used to determine which part of the normal
Lentivirus (protocol in methods above) was then made for all guide libraries and distribution each bin corresponded to, assuming that the leftmost and rightmost
CRISPR-associated vectors (see below). Stable CRISPRi-expressing GM12878, expression bins each exclude the most extreme 0.1% of cells. The guide abundance
BJAB, Jurkat, U937, and THP-1 cell lines and CRISPRa-expressing BJAB, Jurkat, within unsorted cells was quantified and used to estimate guide abundance within
NATURE COMMUNICATIONS |         (2020) 11:1237 | https://doi.org/10.1038/s41467-020-15022-4 | www.nature.com/naturecommunications 9
ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15022-4
the library. A pseudocount was added to each guide count consisting of one read reaction for 6 h at 37 °C followed by Ampure XP purification and elution with
for every 100,000 total reads sequenced in that bin, corresponding to a prior that 55 μL of Buffer EB.
there is no expression difference for cells containing the guide. For each guide, the 10 μL of the mpra:pTNFAIP3:gfp plasmid was electroporated (2 kV, 200Ω,
mean expression for that guide was estimated by maximizing the likelihood of the 25 μF) into 220 μL of 10-beta cells. Electroporated bacteria was split across six tubes
observed guide counts for each bin under this model, given that guide’s overall and each recovered in 2 mL of SOC for 1 h at 37 °C then added to 500 mL of LB
abundance. A z-score was estimated for each guide corresponding to how much the with 100 μg/mL of carbenicillin and grown for 9 h at 37 °C prior to plasmid
mean TNFAIP3 expression of cells containing that guide differed from those purification (Qiagen, 12991). The plasmid prep was then normalized to 1 μg/μL to
containing non-targeting guides by subtracting the mean of the non-targeting generate our final mpra:pTNFAIP3:gfp library used for transfection and lentiviral
guides. delivery.
In order to get element-level statistics, the z-scores for each guide were For all transfections, cells were grown to a density of ~1 × 106 cells/mL and 5 ×
combined in two ways: a significance z-score (proportional to a signed P-value), 107 cells were used for each experiment. Cells were collected by centrifugation at
and an effect-size z-score (the average z-score of guides targeting the element). 300 × g and eluted in 550 μL of RPMI with 55 μg of mpra:pTNFAIP3:gfp library.
Significance z-scores were calculated by applying Stouffer’s method to the Electroporation was performed in 100 μL volumes with the Neon transfection
individual guide’s z-scores. In order to correct these significance z-scores for the system (Life Technologies) applying three pulses of 1200 V for 20 ms each
noise of the assay, they were scaled by the standard deviation of Stouffer z-scores (GM12878) and three pulses of 1325 V for 10 ms each (Jurkat). Using separate
calculated from the non-targeting guides. These scaling factors were calculated control transfections, we achieved transfection efficiencies of 40–60% for all
independently for every number of guides per targeted element n (since the noise replicates. Cells were allowed to recover in 180 mL in RPMI with 15% FBS for 24 h
in the Stouffer z-score depends on the number of guides used to calculate it). For then collected by centrifugation, washed once with PBS, collected and frozen at
example, Stouffer z-scores for elements targeted with n= 5 guides were normalized −80 °C.
by the standard deviation of non-targeting Stouffer z-scores, each calculated from For all transductions, 500 × 106 cells were split into 24-well plates (2M per well
randomly sampled groups of five non-targeting guides. Here, non-targeting in 1 mL of media, 10 plates) infected with lentivirus at an MOI > 1 using polybrene
Stouffer z-scores were calculated by sampling the non-targeting guides into groups (8 μg/mL) using spin transduction (1760 × g, 90 min, 32 °C). Cells were then pooled
of size n, including each non-targeting guide 10 times total, and calculating a set of and centrifuged at 500 × g, the viral supernatant was aspirated, and the cells were
Stouffer z-scores from each sampling, and using the standard deviation of these z- resuspended in fresh media at 5 × 105 cells/mL, and cultured for 4 days
scores to scale the significance z-scores for each element for that n. P-values were maintaining a density between 2 and 10 × 105/mL. Cells were then harvested
then calculated from these z-scores, considering only one-tailed tests through centrifugation at 500 × g, washed with PBS, centrifuged again, and cell
(downregulation for CRISPRi and upregulation for CRISPRa). For an element to be pellets were frozen at −80 °C.
considered significantly regulating TNFAIP3, we required that both replicates’ Total RNA was extracted from cells using Qiagen Maxi RNeasy (Qiagen, 75162)
Benjamini–Hochberg FDRs were less than sqrt(0.1) (i.e. combined FDR < 0.1, and following the manufacturer’s protocol including the on-column DNase digestion. A
both replicates close to significant independently) and for which the direction of second DNase treatment was performed on the purified RNA using 5 μL of Turbo
expression change was identical. In cases where there were more than two DNase (Life Technologies, AM2238) with buffer, in 750 μL of total volume for 1 h
replicates, we included only the two replicates for which the TNFAIP3 promoter at 37 °C. The digestion was stopped with the addition of 7.5 μL 10% SDS and 75 μL
positive control guides showed the strongest effect. Element- and guide-level data of 0.5 M EDTA followed by a 5 min incubation at 70 °C. The total reaction was
are available in Supplementary Data 5, 6. then used for pulldown of GFP mRNA. Water was added to the DNase digested
RNA to bring the total volume to 898 μL to which 900 μL of 20X SSC (Life
Technologies, 15557-044), 1800 μL of Formamide (Life Technologies, AM9342)
16,71 and 2 μL of 100 μM biotin-labeled GFP probe (GFP_BiotinCapture_1-3, IDT,MPRA. MPRA oligosynthesis and cloning was adapted from refs. , tagging Supplementary Data 11) were added and incubated for 2.5 h at 65 °C. Biotin probes
each allele with an average of ~250 DNA barcodes. Oligos were synthesized by were captured using 400 μL of pre-washed Streptavidin beads (Life Technologies,
Agilent Technologies containing 150 bp of genomic context and 15 bp of adapter 65001) eluted in 500 μL of 20X SSC. The hybridized RNA/probe bead mixture was
sequence at either end (5′-ACTGGCCGCTTGACG[150 bp oligo]CACTGCGGC agitated on a nutator at room temperature for 15 min. Beads were captured by
TCCTGC-3′; Supplementary Data 8; 180 bp total). 20 bp barcodes and additional magnet and washed once with 1× SSC and twice with 0.1× SSC. Elution of RNA
adapter sequences were added by performing 28 emulsion PCR reactions each 50 was performed by the addition of 25 μL water and heating of the water/bead
μL in volume containing 1.86 ng of oligo, 25 μL of Q5 NEBNext MasterMix (NEB, mixture for 2 min at 70 °C followed by immediate collection of eluent on a magnet.
M0541S), 1 unit Q5 HotStart polymerase (NEB, M0493S), 0.5 μM MPRA_v3_F A second elution was performed by incubating the beads with an additional 25 μL
and MPRA_v3_20I_R primers (Supplementary Data 11) and 2 ng BSA (NEB, of water at 80 °C. A final DNase treatment was performed in 50 μL total volume
B9000). PCR master mix was emulsified by vortexing with 220 μL Tegosoft DEC using 1 μL of Turbo DNase with addition of the buffer incubated for 60 min at
(Evonik), 60 μL ABIL WE (Evonik) and 20 μL mineral oil (Sigma, M5904) per 50 37 °C followed by inactivation with 1 μL of 10% SDS and purification using RNA
μL PCR reaction at 4 °C for 5 min. 50 μL of Emulsion mixture was added to each clean SPRI beads (Beckman Coulter, A63987).
well of a 96-well plate and cycled with the following conditions; 95 °C for 30 s, 15 First-strand cDNA was synthesized from half of the DNase-treated GFP mRNA
cycles of (95 °C for 20 s, 60 °C for 10 s, 72 °C for 15 s), 72 °C for 5 min. Amplified with SuperScript III and a primer specific to the 3′ UTR (MPRA_v3_Amp2Sc_R,
emulsion mixture was broken and purified by adding 1 mL of 2-butanol (VWR, Supplementary Data 11) using the manufacturer’s recommended protocol,
AA43315-AK), 50 μL of AMPure XP SPRI (Beckman Coulter, A63881) and 80 μL modifying the total reaction volume to 40 μL and performing the elongation step
of binding buffer (2.5 M NaCl, 20% PEG-8000) per 350 μL of Emulsion mix and at 47 °C for 80 min. Single-stranded cDNA was purified by SPRI and eluted in
vigorously vortexed followed by incubation for 10 min at room temperature. 30 μL EB.
Broken emulsion/butanol mixture was spun at 2900 × g for 5 min and the butanol To minimize amplification bias during the creation of cDNA tag sequencing
phase was discarded. The aqueous phase was placed on a magnetic rack for 20 min libraries, samples were amplified by qPCR to estimate relative concentrations of GFP
prior to aspiration. Remaining beads were washed once with 2-butanol, three times cDNA using 1 μL of sample in a 10 μL PCR reaction containing 5 μL Q5 NEBNext
with 80% EtOH and eluted in EB (Qiagen, 19086) to yield our barcoded oligo pool. master mix, 1.7 μL Sybr green I diluted 1:10,000 (Life Technologies, S-7567) and
To create our mpraΔorf library, barcoded oligos were inserted into SfiI digested 0.5 μM of TruSeq_Universal_Adapter and MPRA_Illumina_GFP_F primers
pMPRA-lenti2 (pMPRA-lenti1ΔSfi1; pMPRA-lenti1: Addgene, 61600) by Gibson (Supplementary Data 11). Samples were amplified with the following qPCR
Assembly (NEB, E2611) using 1.1 μg of oligos and 1 μg of digested vector in a 40 μL conditions: 95 °C for 20 s, 40 cycles (95 °C for 20 s, 65 °C for 20 s, 72 °C for 30 s),
reaction incubated for 60 min at 50 °C followed by AMPure XP SPRI purification 72 °C for 2 min. The number of cycles for sample amplification was 1−n (the
and elution in 20 μL of EB. Half of the ligated vector was then transformed into number of cycles it took for each sample to pass the threshold) from the qPCR. To
100 μL of 10-beta e.coli (NEB, C3020K) by electroporation (2 kV, 200Ω, 25 μF). add Illumina sequencing adapters, 10 μL of cDNA samples and mpra:pTNFAIP3:gfp
Electroporated bacteria were immediately split into eight 1 mL aliquots of SOC plasmid control (diluted to the qPCR cycle range of the samples) were amplified
(NEB, B9020S) and recovered for 1 h at 37 °C then independently expanded in using the reaction conditions from the qPCR scaled to 50 μL, excluding Sybrgreen I.
20 mL of LB supplemented with 100 μg/mL of carbenicillin (EMD, 69101-3) on a Amplified cDNA was SPRI purified and eluted in 40 μL of EB. Individual sequencing
floor shaker at 37 °C for 6.5 h. After outgrowth aliquots were pooled prior to barcodes were added to each sample by amplifying the entire 40 μL elution in a
plasmid purification (QIAGEN, 12963). For each of the aliquots we plated serial
6 100 μL Q5 NEBNext reaction with 0.5 μM of TruSeq_Universal_Adapter primer anddilutions after SOC recovery and estimated a library size of ~3.2 × 10 CFUs, a reverse primer containing a unique 8 bp index (Illumina_Multiplex, Supplementary
representing ~250 barcodes per allele. Data 11) for sample demultiplexing post-sequencing. Samples were amplified at
To insert the TNFAIP3 promoter and GFP ORF, 20 μg of mpra:Δorf plasmid 95 °C for 20 s, six cycles (95 °C for 20 s, 64 °C for 30 s, 72 °C for 30 s), 72 °C for 2 min.
was linearized with XbaI (NEB, R0145S) and KpnI-HF (NEB, R3142S) and 1x Indexed libraries were SPRI purified and pooled according to molar estimates from
cutsmart buffer (NEB, B7204S) in a 500 μL volume for 3.5 h at 37 °C, followed by Agilent TapeStation quantifications. Samples were sequenced using 1 × 30 bp
SPRI cleaning. An amplicon containing 165 bp of the TNFAIP3 ORF, GFP open- chemistry on an Illumina HiSeq 2500 or NextSeq.
reading frame and a partial 3′ UTR was then inserted by Gibson assembly using To determine oligo/barcode combinations within the MPRA pool, Illumina
10 μg of XbaI and KpnI linearized mpraΔorf plasmid, 33 μg of the pTNFAIP3/GFP libraries were prepared from the mpraΔorf plasmid library by performing four
amplicon in 400 μL of total volume for 90 min at 50 °C followed by a 1.5× beads/ separate amplifications with 200 ng of plasmid in a 100 μL Q5 NEBNext PCR
sample SPRI purification. The total recovered volume was digested a second time reaction containing 0.5 μM of TruSeq_Universal_Adapter and
to remove remaining uncut vectors by incubation with KpnI and XbaI in a 100 μL MPRA_v3_TruSeq_Amp2Sa_F primers (Supplementary Data 11) with the
10 NATURE COMMUNICATIONS |         (2020)1 1:1237 | https://doi.org/10.1038/s41467-020-15022-4 | www.nature.com/naturecommunications
NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15022-4 ARTICLE
following conditions: 95 °C for 20 s, 6 cycles (95 °C for 20 s, 62 °C for 15 s, 72 °C for GWAS positive results in a null permutation test. These null distributions form a
30 s), 72 °C for 2 min. Amplified material was SPRI purified using a 0.6× bead/ straight line in Fig. 3 because the numerator for both is the number of GWAS tag
sample ratio and eluted with 30 μL of EB. Sequencing indexes were then attached SNPs recovered (nTH) and the denominators for both pseudo-precision and
using 20 μL of the eluted product and the same reaction conditions as for the tag- pseudo-recall are invariant across the randomization (nH for pseudo-precision and
seq protocol, except the number of enrichment cycles was lowered to 5. Samples nT for pseudo-recall).
were molar pooled and sequenced using 2 × 150 bp chemistry on Illumina HiSeq
2500 and NextSeq. Reporting summary. Further information on research design is available in
MPRA RNA output and DNA input sequencing reads were mapped to the the Nature Research Reporting Summary linked to this article.
known tag sequences using a custom python script (quantifyRNATags.py; available
from https://github.com/Carldeboer/MPRAs), allowing for up to four mismatches
within the constant region (the common sequence before the tag sequence) and no Data availability
mismatches within the tag sequence. The barcode counts were input, and tags Raw and processed sequencing data for this study are available on NCBI GEO, under
having fewer than 30 reads in the input (DNA) or 4 reads in the output (RNA) accession “GSE136703”. Other sources for data that support our findings are available
were excluded from subsequent analysis. The log(DNA/RNA) ratio (expression) from: 1000 Genomes, ENCODE, ChIP-Atlas, Immunobase, and GWAS Catalog.
was calculated using raw counts, scaled so that the median expression is 0, and the
expression levels G+C-content normalized such that the mean expression for every
%G+C was 0. Finally, to eliminate instances where the tag sequence modifies the Code availability
apparent expression level, any tags containing any one of eight blackballed 5-mer CRISPR analysis software is available at the following link: https://github.com/
DNA sequences were excluded. Blackballed 5-mers were defined as those for which Carldeboer/MAUDE. Both code and motifs for TF binding motif analysis are available at
the absolute value of the average expression level of all tags containing that 5-mer the following link: https://github.com/Carldeboer/VEP_GOMER. Code for processing
was >0.15. MPRA data is available at the following link: https://github.com/Carldeboer/MPRAs.
SNPs were tested for allele-specific reporter activity by a two-sided Student’s t-
test, comparing the normalized log(RNA/DNA) expression values for the tags for
allele A compared to the tags for allele B. Only SNPs for which we had at least 80 Received: 3 February 2020; Accepted: 17 February 2020;
good tags between the two alleles were tested. P-values were corrected for multiple
hypothesis testing by Benjamini–Hochberg FDR correction. Only SNPs that had an
FDR < 0.1 for at least two of the replicates and where the direction of allele-specific
reporter activity was consistent between all replicates were considered to be
significant.
References
Predicted TF-binding perturbation. In order to find TFs whose motifs were 1. Farh, K. K.-H. et al. Genetic and epigenetic fine mapping of causal
disrupted, both alleles were scanned for each SNP with human and selected mouse autoimmune disease variants. Nature 518, 337–343 (2015).
motifs from CIS-BP72 using VEP56 and a custom VEP module implementing the 2. Gusev, A. et al. Partitioning heritability of regulatory and cell-type-specific
GOMER approach73 for motif scanning (https://github.com/Carldeboer/ variants across 11 common diseases. Am. J. Hum. Genet. 95, 535–552 (2014).
VEP_GOMER). In order to be considered a motif disruption, the region sur- 3. Ulirsch, J. C. et al. Interrogation of human hematopoiesis at single-cell and
rounding the SNP must both be bound on one or both alleles (>95% of >1% MAF single-variant resolution. Nat. Genet. 51, 683–693 (2019).
SNPs), and the binding score difference in the log binding score between the alleles 4. Maurano, M. T. et al. Systematic localization of common disease-associated
must be at least 0.1 (roughly corresponding to about ~1% of SNPs being perturbed variation in regulatory DNA. Science 337, 1190–1195 (2012).
per motif). Both code and motifs for this analysis are available here: https://github. 5. Claussnitzer, M. et al. FTO obesity variant circuitry and adipocyte browning in
com/Carldeboer/VEP_GOMER. humans. N. Engl. J. Med. 373, 895–907 (2015).
6. Gupta, R. M. et al. A genetic variant associated with five vascular diseases is a
Data integration and analysis. In order to gauge how much each assay enriched distal regulator of endothelin-1 gene expression. Cell 170, 522–533.e515
for GWAS signal (as in Fig. 3), we considered all GWAS tag SNPs. Since the set of (2017).
causal SNPs remains unknown, we must use the set of potentially causal SNPs as an 7. Mohanan, V. et al. C1orf106 is a colitis risk gene that regulates stability of
enriched gold standard (e.g. fine mapped variants, or SNPs in tight LD (r2 > 0.8), as epithelial adherens junctions. Science 359, 1161–1166 (2018).
used here). However, a single causal variant could underlie multiple GWAS tag 8. Sankaran, V. G. et al. Human fetal hemoglobin expression is regulated by the
SNPs, for instance, if the causal SNP is in tight LD with both tag SNPs in the developmental stage-specific repressor BCL11A. Science 322, 1839–1842
GWAS population. Although a single GWAS tag SNP could represent multiple (2008).
underlying causal SNPs, we expect this to be uncommon, and a potential expla- 9. Buenrostro, J. D., Giresi, P. G., Zaba, L. C., Chang, H. Y. & Greenleaf, W. J.
nation featuring fewer causal SNPs should be favored. With these considerations in Transposition of native chromatin for fast and sensitive epigenomic profiling
mind, we evaluated each assay for its ability to identify GWAS tag SNPs by being in of open chromatin, DNA-binding proteins and nucleosome position. Nat.
tight LD to hits. Since some assays could not assay every variant, only assayed Methods 10, 1213–1218 (2013).
variants are included. Similarly, if a tag SNP had no assayed SNPs in tight LD, that 10. Degner, J. F. et al. DNase I sensitivity QTLs are a major determinant of human
tag SNP was not included in evaluation of the assay since it could not have been expression variation. Nature 482, 390–394 (2012).
recovered by the assay. The pseudo-precision and pseudo-recall were calculated for 11. Zhao, Z. et al. Circular chromosome conformation capture (4C) uncovers
each assay as described in the “Results” section. The enrichment analysis using extensive networks of epigenetically regulated intra- and interchromosomal
credible sets instead of tag SNPs was performed identically, but instead of evalu- interactions. Nat. Genet. 38, 1341–1347 (2006).
ating the recovery of tag SNPs by being in tight LD to hits, we evaluated the 12. Rao, S. S. P. et al. A 3D map of the human genome at kilobase resolution
recovery of credible sets by having one or more hits within each credible set. reveals principles of chromatin looping. Cell 159, 1665–1680 (2014).
In order to gauge the significance of enrichment for each assay with limited tag 13. Hsu, J. Y. et al. CRISPR-SURF: discovering regulatory elements by
SNPs, we created an empirical null distribution by randomizing the data. Since
deconvolution of CRISPR tiling screen data. Nat. Methods 15, 992–993
some of the assays (e.g. DHS, CRISPRi/a) have an inherent clustering of their hits
(i.e. SNPs within the same enhancer will share the same hit status), our null aimed (2018).
to preserve this clustering. Specifically, the null was derived by ordering the assayed 14. Simeonov, D. R. et al. Discovery of stimulation-responsive immune enhancers
SNPs by genomic position and reassigning hit status Ha,s’=Ha,((s+i) mod n) for every
with CRISPR activation. Nature 549, 111–115 (2017).
possible i (0 < i < n; where n is the number of assayed SNPs and mod is the modulo 15. Mogno, I., Kwasnieski, J. C. & Cohen, B. A. Massively parallel synthetic
operation), and, each time (i.e. for each value of i), calculating pseudo-precision promoter assays reveal the in vivo effects of binding site variants. Genome Res.
and pseudo-recall. P-values represent the fraction of this empirical null with at least 23, 1908–1915 (2013).
as high a pseudo-precision and pseudo-recall as that observed from the actual data. 16. Tewhey, R. et al. Direct identification of hundreds of expression—modulating
We also tested random permutation of the SNP hit status as an alternative to variants using a multiplexed reporter assay. Cell 165, 1519–1529 (2016).
shifting. Here, we used 1000 independent random permutations of the SNP hits to 17. Ulirsch, J. C. et al. Systematic functional dissection of common genetic
create the null model. Although both approaches yielded similar results (Fig. 3 variation affecting red blood cell traits. Cell 165, 1530–1545 (2016).
genomic shifts and Supplementary Fig. 9d, e permutation), we opted to focus on 18. Vockley, C. M. et al. Massively parallel quantification of the regulatory effects
the random null created by shifting hit status; the random permutations fail to of noncoding genetic variation in a human cohort. Genome Res. 25,
capture the clustering of hits that results from genomic proximity and shared hit 1206–1214 (2015).
origins (e.g. adjacent SNPs in the same open chromatin region). For example, if 19. Ma, A. & Malynn, B. A. A20: linking a complex regulator of ubiquitylation to
there was only a single functional enhancer with CRISPRi which contained 10 immunity and human disease. Nat. Rev. Immunol. 12, 774–785 (2012).
SNPs, and the clustering of these SNPs prioritizes only one GWAS signal in one 20. Calderon, D. et al. Landscape of stimulation-responsive chromatin across
region, the distribution of these SNPs randomly could result in as many as 10 diverse human immune cells. Nat. Genet. 51, 1494–1505 (2019).
NATURE COMMUNICATIONS |         (2020) 11:1237 | https://doi.org/10.1038/s41467-020-15022-4 | www.nature.com/naturecommunications 11
ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15022-4
21. Sokhi, U. K. et al. Dissection and function of autoimmunity-associated 50. Danecek, P. et al. The variant call format and VCFtools. Bioinformatics 27,
TNFAIP3 (A20) gene enhancers in humanized mouse models. Nat. Commun. 2156–2158 (2011).
9, 658 (2018). 51. Purcell, S. et al. PLINK: a tool set for whole-genome association and
22. Adrianto, I. et al. Association of a functional variant downstream of TNFAIP3 population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575 (2007).
with systemic lupus erythematosus. Nat. Genet. 43, 253–258 (2011). 52. Leslie, R., O’Donnell, C. J. & Johnson, A. D. GRASP: analysis of
23. Wang, S., Wen, F., Tessneer, K. L. & Gaffney, P. M. TALEN-mediated genotype–phenotype results from 1390 genome-wide association studies and
enhancer knockout influences TNFAIP3 gene expression and mimics a corresponding open access database. Bioinformatics 30, i185–i194 (2014).
molecular phenotype associated with systemic lupus erythematosus. Genes 53. Welter, D. et al. The NHGRI GWAS Catalog, a curated resource of SNP–trait
Immun. 17, 165–170 (2016). associations. Nucleic Acids Res. 42, D1001–D1006 (2014).
24. Wang, S., Wen, F., Wiley, G. B., Kinter, M. T. & Gaffney, P. M. An enhancer 54. Altshuler, D. M. et al. A global reference for human genetic variation. Nature
element harboring variants associated with systemic lupus erythematosus 526, 68–74 (2015).
engages the TNFAIP3 promoter to influence A20 expression. PLoS Genet. 9, 55. Bycroft, C. et al. The UK Biobank resource with deep phenotyping and
e1003750–1003710 (2013). genomic data. Nature 562, 203–209 (2018).
25. Coornaert, B. et al. T cell antigen receptor stimulation induces MALT1 56. McLaren, W. et al. The ensembl variant effect predictor. Genome Biol. 17, 122
paracaspase-mediated cleavage of the NF-κB inhibitor A20. Nat. Immunol. 9, (2016).
263–271 (2008). 57. Zhou, W. et al. Efficiently controlling for case-control imbalance and sample
26. Housley, W. J. et al. Genetic variants associated with autoimmunity drive relatedness in large-scale genetic association studies. Nat. Genet. 50,
NFκB signaling and responses to inflammatory stimuli. Sci. Transl. Med. 7, 1335–1341 (2018).
291ra293–291ra293 (2015). 58. Loh, P.-R., Kichaev, G., Gazal, S., Schoech, A. P. & Price, A. L. Mixed-model
27. Phanstiel, D. H. et al. Static and dynamic DNA loops form AP-1-bound association for biobank-scale datasets. Nat. Genet. 50, 906–908 (2018).
activation hubs during macrophage development. Mol. Cell 67, 1037–1048. 59. Benner, C., Havulinna, A., Salomaa, V., Ripatti, S. & Pirinen, M. Refining fine-
e1036 (2017). mapping: effect sizes and regional heritability. bioRxiv 1–18, https://doi.org/
28. Luppino, J. M. et al. Cohesin promotes interactions across domain boundaries 10.1101/318618 (2018).
to regulate gene expression. bioRxiv 1663, 231–265 (2019). 60. Benner, C. et al. FINEMAP: efficient variable selection using summary data
29. Weber, M. S. et al. Type II monocytes modulate T cell-mediated central from genome-wide association studies. Bioinformatics 32, 1493–1501 (2016).
nervous system autoimmune disease. Nat. Med. 13, 935–943 (2007). 61. Benner, C. et al. Prospects of fine-mapping trait-associated genomic regions
30. Bradshaw, E. M. et al. Monocytes from patients with type 1 diabetes by using summary statistics from genome-wide association studies. Am. J.
spontaneously secrete proinflammatory cytokines inducing Th17 cells. J. Hum. Genet. 101, 539–551 (2017).
Immunol. (Baltim., Md.: 1950) 183, 4432–4439 (2009). 62. Corces, M. R. et al. Lineage-specific and single-cell chromatin accessibility
31. Udalova, I. A., Mantovani, A. & Feldmann, M. Macrophage heterogeneity in charts human hematopoiesis and leukemia evolution. Nat. Genet. 48,
the context of rheumatoid arthritis. Nat. Rev. Rheumatol. 12, 472–485 (2016). 1193–1203 (2016).
32. Rawlings, D. J., Metzler, G., Wray-Dutra, M. & Jackson, S. W. Altered B cell 63. Servant, N. et al. HiC-Pro: an optimized and flexible pipeline for Hi-C data
signalling in autoimmunity. Nat. Publ. Group 17, 421–436 (2017). processing. Genome Biol. 16, 259 (2015).
33. Bluestone, J. A., Bour-Jordan, H., Cheng, M. & Anderson, M. T cells in the 64. Lareau, C. A. & Aryee, M. J. hichipper: a preprocessing pipeline for calling
control of organ-specific autoimmunity. J. Clin. Investig. 125, 2250–2260 DNA loops from HiChIP data. Nat. Methods 15, 155–156 (2018).
(2015). 65. Schep, A. N. et al. Structured nucleosome fingerprints enable high-resolution
34. Rosenblum, M. D., Gratz, I. K., Paw, J. S. & Abbas, A. K. Treating human mapping of chromatin architecture within regulatory regions. Genome Res. 25,
autoimmunity: current practice and future prospects. Sci. Transl. Med. 4, 1757–1770 (2015).
125sr121–125sr121 (2012). 66. Doench, J. G. et al. Optimized sgRNA design to maximize activity and
35. Finucane, H. K. et al. Heritability enrichment of specifically expressed genes minimize off-target effects of CRISPR-Cas9. Nat. Biotechnol. 34, 184–191
identifies disease-relevant tissues and cell types. Nat. Genet. 50, 621–629 (2016).
(2018). 67. Langmead, B., Trapnell, C., Pop, M. & Salzberg, S. L. Ultrafast and memory-
36. Matmati, M. et al. A20 (TNFAIP3) deficiency in myeloid cells triggers erosive efficient alignment of short DNA sequences to the human genome. Genome
polyarthritis resembling rheumatoid arthritis. Nat. Genet. 43, 908–U129 Biol. 10, R25 (2009).
(2011). 68. Hsu, P. D. et al. DNA targeting specificity of RNA-guided Cas9 nucleases. Nat.
37. Song, X.-T. et al. A20 is an antigen presentation attenuator, and its inhibition Biotechnol. 31, 827–832 (2013).
overcomes regulatory T cell–mediated suppression. Nat. Med. 14, 258–265 69. Fulco, C. P. et al. Activity-by-contact model of enhancer-promoter regulation
(2008). from thousands of CRISPR perturbations. Nat. Genet. 51, 1664–1669 (2019).
38. Maelfait, J. et al. A20 (Tnfaip3) deficiency in myeloid cells protects against 70. Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat.
influenza A virus infection. PLoS Pathog. 8, e1002570–1002511 (2012). Methods 9, 357–359 (2012).
39. Hammer, G. E. et al. Expression of A20 by dendritic cells preserves immune 71. Melnikov, A. et al. Systematic dissection and optimization of inducible
homeostasis and prevents colitis and spondyloarthritis. Nat. Immunol. 12, enhancers in human cells using a massively parallel reporter assay. Nat.
1184–1193 (2011). Biotechnol. 1–9, https://doi.org/10.1038/nbt.2137 (2012).
40. Chu, Y. et al. B cells lacking the tumor suppressor TNFAIP3/A20 display 72. Weirauch, M. T. et al. Determination and inference of eukaryotic
impaired differentiation and hyperactivation and cause inflammation and transcription factor sequence specificity. Cell 158, 1431–1443 (2014).
autoimmunity in aged mice. Blood 117, 2227–2236 (2011). 73. Granek, J. A. & Clarke, N. D. Explicit equilibrium modeling of transcription-
41. Mumbach, M. R. et al. Enhancer connectome in primary human cells factor binding and gene regulation. Genome Biol. 6, R87 (2005).
identifies target genes of disease-associated DNA elements. Nat. Genet. 49,
1602–1612 (2017).
42. Oki, S. et al. ChIP-Atlas: a data-mining suite powered by full integration of
public ChIP-seq data. EMBO Rep. 19, e46255 (2018). Acknowledgements
43. Wang, G., Sarkar, A. K., Carbonetto, P. & Stephens, M. A simple new Work was funded by NHGRI R01HG008131-01 (NH), the Klarman Cell Observatory
approach to variable selection in regression, with application to genetic fine- and HHMI (AR). JPR is funded by NIH 5 F32 AI129249. CGD is funded by a CIHR and
mapping. bioRxiv 1–41, https://doi.org/10.1101/501114 (2018). NIH K99-HG009920-01. NH is the David P. Ryan endowed chair, and thanks Sandra,
44. Tehranchi, A. et al. Fine-mapping cis-regulatory variants in diverse human Sarah and Arthur Irving for support. RT is supported by K99HG008179. UK Biobank
populations. eLife 8, 1330 (2019). analyses were conducted via application 31063. We thank BM Neale and his colleagues
45. Tehranchi, A. K. et al. Pooled ChIP-Seq links variation in transcription factor who provided scripts and resources for the UK Biobank analyses.
binding to complex disease risk. Cell 165, 730–741 (2016).
46. McGovern, A. et al. Capture Hi-C identifies a novel causal gene, IL20RA, in
the pan-autoimmune genetic susceptibility region 6q23. Genome Biol. 17, 212
(2016). Author contributions
47. Wu, J. et al. CRISPR/cas9 mediated knockout of an intergenic variant JPR, CGB, NH conceived the project. JPR and CGB planned experiments. JPR, LL, SKR,
rs6927172 identified IL-20RA as a new risk gene for multiple autoimmune DTB, RI carried out experiments. RT assisted with MPRA experimental protocol. CPF
diseases. Genes Immun. 17, 160 (2018). and JME assisted with noncoding CRISPR screening and FlowFISH protocol. CGB
48. Huang, H. et al. Fine-mapping inflammatory bowel disease loci to single- carried out experimental analyses. HKF, CAL, MK, JCU carried out heritability
variant resolution. Nature 547, 173–178 (2017). enrichment and fine-mapping genetic analyses. NH and AR supervised work. JPR, CGB,
49. Li, H. Tabix: fast retrieval of sequence features from generic TAB-delimited AR, and NH wrote manuscript with feedback from ESL. All authors provided feedback
files. Bioinformatics 27, 718–719 (2011). on the manuscript prior to submission.
12 NATURE COMMUNICATIONS |         (2020) 11:1237 | https://doi.org/10.1038/s41467-020-15022-4 | www.nature.com/naturecommunications
NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15022-4 ARTICLE
Competing interests Reprints and permission information is available at http://www.nature.com/reprints
AR is a co-founder and equity holder of Celsius Therapeutics, a founder of Immunitas,
and SAB member of ThermoFisher Scientific, Asimov Neogene Therapeutics, and Syros Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in
Pharmaceuticals. NH is a co-founder and equity holder of Neon Therapeutics. ESL serves published maps and institutional affiliations.
on the Board of Directors for Codiak BioSciences and Neon Therapeutics, and serves on
the Scientific Advisory Board of F-Prime Capital Partners and Third Rock Ventures; he is
also affiliated with several non-profit organizations including serving on the Board of Open Access This article is licensed under a Creative Commons Attri-
Directors of the Innocence Project, Count Me In, and Biden Cancer Initiative, and the bution 4.0 International License, which permits use, sharing, adaptation,
Board of Trustees for the Parker Institute for Cancer Immunotherapy. ESL has served distribution and reproduction in any medium or format, as long as you give appropriate
and continues to serve on various federal advisory committees.
credit to the original author(s) and the source, provide a link to the Creative Commons
license, and indicate if changes were made. The images or other third party material in this
Additional information article are included in the article’s Creative Commons license, unless indicated otherwise in
Supplementary information is available for this paper at https://doi.org/10.1038/s41467- a credit line to the material. If material is not included in the article’s Creative Commons
020-15022-4. license and your intended use is not permitted by statutory regulation or exceeds the
permitted use, you will need to obtain permission directly from the copyright holder. To
Correspondence and requests for materials should be addressed to A.R. or N.H. view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
Peer review information Nature Communications thanks the anonymous reviewer(s) for
their contribution to the peer review of this work. © The Author(s) 2020
NATURE COMMUNICATIONS |         (2020) 11:1237 | https://doi.org/10.1038/s41467-020-15022-4 | www.nature.com/naturecommunications 13