COVID-19 Susceptibility By Race, Blood Type

Social, demographic, and pre-existing conditions evaluated as risk factors for COVID-19
hospitalization
To explore non-genetic factors associated with hospitalization for COVID-19, age, sex, ancestry,
median household income of the residential zip code, educational attainment, body mass index,
and pre-existing conditions were evaluated as risk factors in bivariate and multivariate logistic
regression models. Ancestry was inferred via a previously described genetic ancestry
classification algorithm9,10 to ensure compatibility with the GWAS methodology described below.

As social and demographic factors are strongly associated with COVID-19 hospitalization and
many of the pre-existing conditions, final models estimating the relationship between the preexisting
condition and COVID-19 hospitalization were adjusted for age (10-year increments),
sex, ancestry, high school or below education, and median household income of the residential
zip code (in $10,000 increments). Finally, a multivariate logistic regression model was specified
to quantify the risk of the top conditions after adjusting for each other, and for sociodemographic
characteristics related to risk of hospitalization for COVID-19.

The relationship between age and COVID-19 and hospitalization was determined by
categorizing cases by 10-year age increments between 30 and 80, and then calculating the
percentage of cases in each age group and the percentage of cases which reported
hospitalization. To describe differences in hospitalization by ancestry, age-standardized
estimates were calculated by applying the percent of cases hospitalized within age strata for
European, African American, and Latino respondents and applying that percentage to the age
structure of all cases in the study population. While our data also included respondents of other
ancestries, such as East Asian and South Asian, the sample sizes for these populations were
too small for robust inferences to be performed. All analyses were conducted in R statistical
software, version 3.3.

Phenotype definitions for GWAS
Using the information derived from the surveys, we defined a set of phenotypes that aimed to
capture aspects of COVID-19 diagnosis and severity. Following preliminary analysis, we
selected one ‘diagnosis’ phenotype that contrasts positive and negative outcomes from a
COVID-19 test, and four phenotypes that capture aspects of COVID-19 disease ‘severity’,
including pneumonia, hospitalization, or the need for respiratory support (Supplementary Table
2). For the severity phenotypes, we explored different choices of how to define controls, and
found that using controls who had not reported a COVID-19 diagnosis (were neither diagnosed
with nor tested positive for COVID-19) provided a large-scale control set that appeared to
maximize power in the GWAS.

Genotyping and SNP imputation
DNA extraction and genotyping were performed on saliva samples by CLIA-certified and CAPaccredited
clinical laboratories of Laboratory Corporation of America. Samples were genotyped
on one of five genotyping platforms. The V1 and V2 platforms were variants of the Illumina
HumanHap550 + BeadChip and contained a total of about 560,000 SNPs, including about
25,000 custom SNPs selected by 23andMe. The V3 platform was based on the Illumina
OmniExpress + BeadChip and contained a total of about 950,000 SNPs and custom content to
improve the overlap with our V2 array. The V4 platform was a fully custom array of about
950,000 SNPs and included a lower redundancy subset of V2 and V3 SNPs with additional
coverage of lower-frequency coding variation. The V5 platform was based on the Illumina GSA
array, consisting of approximately 654,000 pre-selected SNPs and approximately 50,000
custom content variants. Samples that failed to reach 98.5% call rate were re-analyzed.
Individuals whose analyses failed repeatedly were re-contacted by 23andMe customer service
to provide additional samples, as is done for all 23andMe customers.

Participant genotype data were imputed using the Haplotype Reference Consortium (HRC)
panel11, augmented by the Phase 3 1000 Genomes Project panel12 for variants not present in
HRC. We phased and imputed data for each genotyping platform separately. For the nonpseudoautosomal
region of the X chromosome, males and females were phased together in
segments, treating the males as already phased; the pseudoautosomal regions were phased
separately. We then imputed males and females together, treating males as homozygous
pseudo-diploids for the non-pseudoautosomal region.

Genome-wide association study (GWAS)
Genotyped participants were included in GWAS analyses on the basis of ancestry as
determined by a genetic ancestry classification algorithm10. For each phenotype, we selected a
set of unrelated individuals so that no two individuals shared more than 700cM of DNA identical
by descent (IBD). For case-control phenotypes, if a case and a control were identified as having
at least 700cM of DNA IBD, we preferentially discarded the control from the sample.

For case-control comparisons, we tested for association using logistic regression, assuming
additive allelic effects. For tests using imputed data, we use the imputed dosages rather than
best-guess genotypes. We included covariates for age, age squared, sex, a sex:age interaction,
the top ten principal components to account for residual population structure, and dummy
variables to account for genotyping platform. The association test p-value was computed using
a likelihood ratio test, which in our experience is better behaved than a Wald test on the
regression coefficient. Results for the X chromosome were computed similarly, with men coded
as if they were homozygous diploid for the observed allele.

We ran GWAS for each phenotype separately, and combined both genotyped and imputed
data. When choosing between imputed and genotyped GWAS results, we favored the imputed
result, unless the imputed variant was unavailable or failed quality control (QC). For imputed
variants, we removed variants with low imputation quality (r2 < 0.5 averaged across batches, or
a minimum r2 < 0.3) or with evidence of batch effects (ANOVA F-test across batches, p-value <
10-50). For genotyped variants, we removed variants only present on our V1 or V2 arrays (due to
small sample size) that failed a Mendelian transmission test in trios (p-value < 10-20), that failed
a Hardy-Weinberg test in Europeans (p-value < 10-20), failed a batch effect test (ANOVA p-value
< 10-50), or had a call rate < 90%.

We repeated the GWAS analysis separately in each population cohort for which we had
sufficient data (European, Latino, African American), and then performed trans-ethnic metaanalysis
using a fixed effects model (inverse variance method13), restricting to variants of at
least 1% minor allele frequency.

Within each GWAS, we identified regions with genome-wide significant (GWS) associations. We
define the region boundaries by identifying all SNPs with p-value < 10−5 within the vicinity of a
GWS association, and then grouping these regions into intervals so that no two regions are
separated by less than 250 kb. We consider the SNP with the smallest p-value within each
interval to be the index SNP. We also annotated our findings based on linkage disequilibrium
(LD) with results from published GWAS, coding variation, and expression quantitative trait loci
(eQTL), specifically by finding annotations with r2 > 0.5 and within 500kb of the index SNP.