Identification of Immune complement function as a determinant of adverse SARS-CoV-2 infection outcome

Ethics and Data Governance Approval

The study is approved by the Columbia University Irving Medical Center Institutional Review Board (IRB# AAAL0601) and the requirement for an informed consent was waived. A data request associated with this protocol was submitted to the Tri-Institutional Request Assessment Committee (TRAC) of New- York Presbyterian, Columbia, and Cornell and approved. The research on the UK Biobank data has been conducted using the UK Biobank Resource under Application Number 41039. The transcriptomics analysis samples were collected and processed through the Weill Cornell Medicine Institutional Review Board (IRB) Protocol 19-11021069.

Cohort and Study Description

In this observational cohort study, we used a data warehouse derived from electronic health records (EHRs) from 11,116 patients treated at New York-Presbyterian/Columbia University Irving Medical Center for suspected cases of SARS-CoV-2 infection. For these patients we collected contemporary data from their current encounter (i.e. the encounter associated with their suspected SARS-CoV-2 infection) as well as historical data, if available, from their previous encounters. Contemporary data (data collected between February 1, 2020 and April 12, 2020) included insurance billing information, laboratory measurements, procedures, and SARS-CoV-2 diagnostic test results. These data were derived from the data warehouse tables in Epic. 6,927 patients have historical data (data collected prior to September 24, 2019) available from an OMOP v5 instance stored using MySQL, which included all of the standard tables for recording condition, procedure, medication, and measurement data (among others). Of these we used the insurance billing information from the condition occurrence table and demographics from the person table. See Preparation of data for modeling for further details on data preparation.

We used the contemporary data to define inclusion criteria and outcomes (requiring mechanical respiration and mortality) and used historical data to define patient comorbidities. We defined three hypothesized comorbidity covariates, macular degeneration, complement deficiency disorders, and disorders of coagulation. We used historical data to define these comorbidities, age, and sex. We did not include race and ethnicity data in the modeling as we have previously found issues with the data quality²⁸. The race/ethnicity data we do have is included in the tables for reference. We also modeled other comorbidities previously associated with morbidity and mortality (Zhou et al and others), including history of cardiovascular disease, hypertension, obesity, and diabetes (Table 1, Table S1) -- all derived from the historical data. Coded covariate definitions, as well as lists of which diagnosis codes are most common in each group, are available in the supplemental materials and methods. We used established institutional procedures and an institutional clinical data warehouse to extract all data from the EHR.

Defining patient outcomes

Outcome definitions were defined by data derived from the electronic health record between February 1, 2020 and April 12, 2020. Mortality is derived from a death note filed by a resident or primary provider that records the date and time of death. Intubation was used as an intermediary endpoint and is a proxy for a patient requiring mechanical respiration. We used note types that were developed for patients with SARS-CoV-2 infection to record that this procedure was completed. We validated outcome data derived from notes against the patient’s medical record using manual review.

Preparation of data for modeling

We used MySQL and python libraries (pymysql, pandas) to extract and prepare the data for modeling. The code for data preparation is available in the github (https://github.com/tatonettilab/complementcovid) as a Jupyter Notebook titled Data Setup. We begin by creating a master list of suspected covid patients. These are patients that are either diagnosed with the disease, as indicated by a ICD10 code for SARS-CoV-2 infection, in their billing data or a patient that was tested for the presence of the virus using RT-PCR as indicated by a “lab” order for the test. We found 2,821 using the former method and 11,116 patients using the latter. We then extracted birthdates, death dates (if the patient had died or a null value otherwise), and sex codes (1 for female, 2 for male). Patients which had sex codes for non-binary genders were excluded from our analysis. We then define a “first diagnosis date” for each patient as either their first diagnosis date (by billing code) or the first date that they tested positive for SARS-CoV-2, whichever comes first. Next, we calculate each patient’s age at the time of this “first diagnosis date.” Each of the outcomes and covariates are extracted from their respective tables as detailed in the github. Whenever possible, we use the highest-level ancestor code (from the structured vocabulary in OMOP) that represents the concept we want to model. We then use the concept ancestor tables to grab all the descendant codes. Note that diabetic kidney disease was considered for inclusion and so is represented in the data preparation script, however, it was never modeled. Cough is included as a covariate as a reference symptom for comparison. The last step in the preparation process was to compute the censor dates. To do, we iterated through each patient in our master list and computed their time (in days) to intubation (if they required mechanical respiration) or death (if they died). If not, then the study end date (April 25, 2020) was used as the patient’s censored time (in days). Finally, for any patients that were not SARS-CoV-2 positive, their time-to-event values were set to a null indicator to be dropped from the dataset later. Finally, the data are all combined in a pandas (version 1.0.3) dataframe and saved to disk as a pickle file for efficient loading.

Statistical Model

Our patient timelines may be censored since our study cohort included patients that were being treated at the time of analysis. We performed survival analysis on the intubation orders and death using a Cox proportional-hazards model and visualized the risk using Kaplan-Meier curves using the lifelines python package (version 0.24.4). Error estimates on the Kaplan-Meier curves are estimated using Greenwood’s Exponential Formula²⁹. We fit both univariate models and models fit on the covariate, age, and sex and used log-likelihood to assess significance. We reported Cox proportional hazards coefficients and their 95% confidence intervals (Table 1). We modeled whether or not a patient had macular degeneration, a complement deficiency disorder, or a coagulation disorder as binary variables (1=yes, 0=no). Code definitions provided in Table S1. We also included other significant comorbidities suggested by previous studies, CAD, hypertension, T2DM, obesity, or smoking status as binary variables (1=yes, 0=no), sex as a binary variable (0=female, 1=male), age as quantitative variable, older age over 65 (note that age over 65 is used only for illustrative purposes and is not used in multivariate modeling -- in the multivariate model age as a quantitative variable is used), and outcome as a binary variable (1=yes, 0=no). The outcome of interest was coded as 0 until the day it occurred (the date of the first intubation order following admission or the death date) or the date of analysis, whichever occurred first. Survival curves are generated for the indicated variables by setting all other variables to their respected averages within the training data. Note that we dropped patients who experienced the outcome before their initial diagnosis. This is either due to patients being hospitalized prior to infection (in the case of intubation) or errors in the coded data. We dropped 121 patients for intubation prior to infection and 12 patients for prior death. We also restricted the study to 90 days from the start date. One patient was removed for having an event outside of this range.

Covariate Correlations

Using the data prepared as discussed above, we computed pairwise statistical correlations between age, sex as well as history of macular degeneration, complement deficiency disorders, coagulation disorders, HTN, T2DM, obesity, and CAD. We computed them using data from all suspected patients (tested both positive and negative) as well as only those patients who tested positive. We used spearman rho and the tanimoto coefficients (1-Jaccard distance) as our measures of correlation. For the comparison using the tanimoto coefficient we binarized age as greater than or equal to 65.

Statistical Software

We used Jupyter Notebooks (jupyter-client version 5.3.4 and jupyter-core version 4.6.1) running Python 3.7 and all fit models using the python lifelines package (version 0.24.4).

Sample Collection and Processing

Patient specimens were collected with patients’ consent at New York Presbyterian Hospital (NYPH) and then processed for RT-PCR as described previously³⁰. Nasopharyngeal (NP) swab specimens were collected using the BD Universal Viral Transport Media system (Becton, Dickinson and Company, Franklin Lakes, NJ) from symptomatic patients.

Extraction of Viral RNA and RT-PCR detection

Total viral RNA was extracted from deactivated samples using automated nucleic acid extraction on the QIAsymphony and the DSP Virus/Pathogen Mini Kit (QIAGEN). One step reverse transcription to cDNA and real-time PCR (RT-PCR) amplification of viral targets, E (envelope) and S (spike) genes and internal control, was performed using the Rotor-Gene Q thermocyler (QIAGEN).

Human Transcriptome Analysis

RNA-seq reads that mapped unambiguously to the human reference genome via Kraken2 were used to detect transcriptional responses to SARS-CoV-2 infection as described previously³⁰. Briefly, reads were trimmed with TrimGalore, aligned with STAR (v2.6.1d) to the human reference build GRCh38 and the GENCODE v33 transcriptome reference, gene expression was quantified using featureCounts, stringTie and salmon using the nf-core RNAseq pipeline. Sample QC was reported using fastqc, RSeQC, qualimap, dupradar, Preseq and MultiQC. Reads, as reported by featureCounts, were normalized using variance- stabilizing transform (vst) in DESeq2 package in R and DESeq2 was used to call differential expression with either Positive cases vs Negative, or viral load (High/Medium/Low/None) as reported by RT-PCR cycle threshold (Ct) values. Transcript counts (per million) were used to rank genes and perform gene set enrichment analysis (GSEA).

Reverse Transcriptase, quantitative real-time PCR (RT-PCR)

The presence of SARS-CoV-2 in clinical samples was determined by RT-PCR. Briefly, primers for the E (envelope) gene (which detects all members of the lineage B of beta-CoVs), and the S (spike) gene (which specifically detect SARS-CoV-2). Samples were annotated using RT-PCR cycle threshold (Ct) value for SARS-CoV-2 primers as follows: Ct < 18 were assigned “high viral load"; Ct 18 - 24 were assigned “medium viral load”; and Ct 24 - 40 were assigned “low viral load” stratifications; Ct > 40 was classified as negative (-).

Data Source

UK Biobank subjects that were of White British descent, in the UK Biobank PCA calculations and therefore without 3rd degree and above relatedness and without aneuploidy, were used in this study, totaling 337,147 subjects (181,032 females and 156,115 males) (Bycroft 2018). Of the nearly 500,000 participants, approximately 50,000 subjects were genotyped on the UK BiLEVE Array by Affymetrix while the rest were genotyped using the Applied Biosystems UK Biobank Axiom Array, with over 800,000 markers using build GRCh37 (hg19). The arrays share 95% marker coverage. We extracted markers with a minor allele frequency greater than 0.005, INFO score greater than 0.3, and Hardy- Weinberg equilibrium test mid-p value greater than 10-10 using PLINK2³¹. UKBB version 3 Imputation combined the Haplotype Research Consortium with the UK10K haplotype resource using the software IMPUTE4 (UK Biobank White paper). Association analyses were performed using a logistic regression model with additive gene dosage and covariates including age at 2018, sex, first 10 principal components (provided by the UK Biobank), and the genotyping array the sample was carried out on. We determined the alpha threshold for study-wide significance using an empirical permutation analysis (see Empirical Permutation Evaluation to set Alpha Thresholds). We performed a study-wide association analysis comparing variants for subjects that were SARS-CoV-2 positive and required hospitalization against the entire population of 337,147 subjects

Targeted Gene Set Definition

The union of coagulation and complement related gene sets (with immunoglobulin genes removed) that are part of MsigDB was used to define the set of 102 genes used in this study. For each gene, we used the transcriptional start and stop site from the hg19 build of the human genome to define a catchment window of 80kbp. From the 805,426 variants profiled in the UK Biobank genotyping data after quality control and QC filters using PLINK2 (see above), 3,540 variants within the transcribed region of the genes of interest or within 80kbp flanking the transcribed region, 2,888 are within 60kbp, 2,292 are within 40kbp, and 936 are located in haplotype blocks with study genes.

Empirical Permutation Evaluation to set Study-wide Alpha Thresholds

We used permutation to estimate null distributions of the number of hits expected at 9 alpha thresholds varying from (5×10^-5 to 0.05) and by varying the distance threshold from 40kb to 80kb. As shown previously, 80% of GWAS hits are within 60Kb of the nearest gene³². Further, as shown in Supplemental Figure S4, we empirically determined that the majority of eQTLs (>70%) are within 60kb of gene bodies. We performed two sets of permutation analyses: (i) permuted the initial set of genes on which the included variant loci were chosen and (ii) permuted the case/control labels. We repeated each 100 times and used the resulting data to fit a negative binomial distribution as our estimate of the null. Additionally, we evaluated each of the sampled SNP variant sets from (i) and compared their MAF distribution with the MAF distribution of the Complement and Coagulation set. We removed any sets that were significantly different (nominal p-value < 0.05) according to a Mann-Whitney U test (52 of 100 sets were removed due to this criterion; see Supplemental Figure S5). We found that the negative binomial fit the data the best according to a goodness of fit test (Supplemental Figure S2). We used this distribution to assess statistical significance for each combination of alpha and distance values. The result is two estimates of the significance for each alpha (α), distance (d) pair, and , from permutation analyses (i) and (ii) above, respectively. For example: where is the number of permutation loci with a p-value under the threshold, α. The parameters r and p of the negative binomial represent the number of successes/failures and the probability of success, respectively. Both r and p are fit using non-linear least squares (the curve_fit function in scipy.optimize) on , the count data from the permutation analyses for the given α and d. The P is then calculated using the CDF of the fitted negative binomial distribution.

For the gene set permutation analysis (i.e. (i) above) we evaluated each of the 100 replicates to confirm that the minor allele frequency distribution was statistically indistinguishable from that of the complement and coagulation gene set variants. We did so by performing a Mann-Whitney U test between the two distributions and excluded any replicates that showed a significant difference (nominal p-value < 0.05). 52 replicates were excluded because of this requirement (Figure SX). This MAF distribution analysis is not necessary for the case/control permutation analysis (i.e. (ii) above) as the loci are the same in each replicate and it is the case/control labels that are permuted.

Finally, to set the study-wide alpha for each study we chose the greatest threshold value that was gave a P of 0.05 or less for both permutation analysis method:

Finally, this entire process was repeated for two cohorts of patients, (a) the initial COVID cohort released by the UK Biobank in April 2020 and (b) the updated COVID cohort released in May 2020. The chosen a for April was 0.001 and the chosen α for May was 0.0025. A data file of all of the distribution fit results and their resulting chi-squared goodness-of-fit statistics is made available in the supplemental materials.

We also performed this permutation significance estimation for the haplotype-derived SNP sets although the distances for all loci chosen using that method are below the minimum in this analysis of 40Kb so those results are constant with regards to distance (Figures S3a-b). The chosen α for the LD-derived SNP sets is 0.01 and 0.0075 for April and May, respectively.

Haplotype block-based selection of SNPs

We identified haplotype blocks based on linkage disequilibrium within the UK Biobank data genotype data of the 337,147 subjects using PLINK1.9, where the lower 90% confidence interval is greater than 0.70 and the upper 90% confidence interval is at least 0.98. We identified blocks of interests, and subsequently the variants within those blocks, as those that contain any part of the genes of interest as denoted by the transcriptional start and end sites from the hg19 build of the human genome. From the 805,426 variants profiles in the UK Biobank genotype data, we identified 7,281 variants within the genes of interest. After applying additional QC filters using PLINK2, 936 variants remained for analysis.

Software

We used PLINK v2.00a2LM 64-bit Intel (26 Aug 2019) to run the genetic association analysis. We used PLINK v1.90b6.10 64-bit (17 Jun 2019) to identify haplotype blocks based on linkage disequilibrium. We used Jupyter Notebooks (jupyter-client version 5.3.4 and jupyter-core version 4.6.1) running Python 3.7, numpy 1.18.1, and scipy 1.4.1 for the permutation analyses.