A systems biology approach identifies candidate drugs to reduce mortality in severely ill patients with COVID-19

  1. Vinicius M. Fava
  2. Mathieu Bourgey
  3. Pubudu M. Nawarathna
  4. Marianna Orlova
  5. Pauline Cassart
  6. Donald C. Vinh
  7. Matthew Pellan Cheng
  8. Guillaume Bourque
  9. Erwin Schurr
  10. David Langlais

Reviewed by

Read the full article

Listed in

Log in to save this article

Abstract

Despite the availability of highly efficacious vaccines, coronavirus disease 2019 (COVID-19) caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) lacks effective drug treatment, which results in a high rate of mortality. To address this therapeutic shortcoming, we applied a systems biology approach to the study of patients hospitalized with severe COVID. We show that, at the time of hospital admission, patients who were equivalent on the clinical ordinal scale displayed significant differential monocyte epigenetic and transcriptomic attributes between those who would survive and those who would succumb to COVID-19. We identified messenger RNA metabolism, RNA splicing, and interferon signaling pathways as key host responses overactivated by patients who would not survive. Those pathways are prime drug targets to reduce mortality of critically ill patients with COVID-19, leading us to identify tacrolimus, zotatifin, and nintedanib as three strong candidates for treatment of severely ill patients at the time of hospital admission.

Article activity feed

  1. Version published to 10.1126/sciadv.abm2510
  2. ScreenIT

    SciScore for 10.1101/2021.09.14.21262309: (What is this?)

    Please note, not all rigor criteria are appropriate for all manuscripts.

    Table 1: Rigor

    NIH rigor criteria are not applicable to paper type.

    Table 2: Resources

    Software and Algorithms
    SentencesResources
    The DoubletFinder package was used to identify and filter out the remaining cell doublets (66).
    DoubletFinder
    suggested: (DoubletFinder, RRID:SCR_018771)
    Single-cell RNA-seq differential gene expression (DGE) analysis: Gene expression changes across the two comparisons “Deceased vs Alive” at day 0 and “Critically vs Moderate” ill at follow-up (days 5 and 15) was evaluated in each lineage using the MAST approach adjusting for Age and Sex covariates (69).
    MAST
    suggested: (MAST, RRID:SCR_016340)
    Single-cell RNA-seq pathway enrichment analysis: Enrichment analyses for GO biological processes were performed separately on up- and down-regulated DGE genes and for each cell sub-population using the enrichGO function of the clusterProfiler R package v3.13 (70).
    clusterProfiler
    suggested: (clusterProfiler, RRID:SCR_016884)
    Enrichment heatmaps were generated based on these selected results using the pheatmap R package v1.0.12.
    pheatmap
    suggested: (pheatmap, RRID:SCR_016418)
    ATAC-seq data analyses: The fifteen ATAC-seq libraries were checked for reads quality and low-quality reads were removed and adaptor sequences trimmed with Cutadapt v1.18 (73).
    Cutadapt
    suggested: (cutadapt, RRID:SCR_011841)
    Reads were then aligned to the human genome (hg38) using BWA v0.7.17 default parameters (74).
    BWA
    suggested: (BWA, RRID:SCR_010910)
    Next, samtools v1.9 was used to remove reads aligned to mitochondria or alternative contigs (76).
    samtools
    suggested: (SAMTOOLS, RRID:SCR_002105)
    AlignmentSieve from DeepTools v3.5.0 was used to filter for unique paired reads (with --samFlagExclude 1804), to remove ENCODE hg38 blacklisted regions and to select fragments between 40 and 2000bp in length (77, 78).
    DeepTools
    suggested: (Deeptools, RRID:SCR_016366)
    Library quality metrics were summarized with MultiQC v1.8 (79).
    MultiQC
    suggested: (MultiQC, RRID:SCR_014982)
    Chromatin accessibility was determined using MACS2 v2.2.6 callpeak with q value < 0.05 in BAMPE mode (80).
    MACS2
    suggested: (MACS, RRID:SCR_013291)
    We merged overlapping narrow peaks detected in at least two samples with bedtools v2.26 (81).
    bedtools
    suggested: (BEDTools, RRID:SCR_006646)
    For each library, featureCounts v1.6.3 was used to count the number of unique fragments overlapping the targeted regions and determine the fraction of reads in peaks (FRIP) (82).
    featureCounts
    suggested: (featureCounts, RRID:SCR_012919)
    We used filterByExpr from edgeR v3.30.3 to select 13,398 peaks with more than ten fragments in 80% of the 15 libraries.
    edgeR
    suggested: (edgeR, RRID:SCR_012802)
    We estimated Storey’s false discovery rate (FDR) with qvalue v2.20.0 package (86).
    qvalue
    suggested: (Qvalue, RRID:SCR_001073)
    ATAC-seq motif enrichment, gene ontology and pathway analysis: To test for enrichment of TF binding motifs in DAC regions, we used the findMotifsGenome function implemented in HOMER v4.11 (88).
    HOMER
    suggested: (HOMER, RRID:SCR_010881)
    By combining the TSS and GeneHancer approaches 76.7% of the 13,398 tested peaks were assigned to at least one gene.
    GeneHancer
    suggested: None
    Next, SparK v2.6.2 was used to calculate the average base-pair accessibility and standard deviation per group (93).
    SparK
    suggested: (Spark, RRID:SCR_006207)
    We used removeBatchEffect from limma to calculate the residual log2 CPM after regressing out covariates included in the linear models (84).
    limma
    suggested: (LIMMA, RRID:SCR_010943)
    Raincloud, ggpubr and ggplot2 were used to produce the remaining ATAC-seq plots (95, 96).
    ggplot2
    suggested: (ggplot2, RRID:SCR_014601)
    In summary, raw reads were trimmed for quality (quality score ≥30) using trimmomatic (v0.36) and to remove sequencing adaptors (73), and then aligned to the human genome hg38 reference genome using Dynamic Read Analysis for Genomics (DRAGENTM) epigenome pipeline (v3.4.5) which utilizes Illumina’s DRAGENTM Bio-It Platform with the improved and highly optimized mapping algorithms.
    trimmomatic
    suggested: (Trimmomatic, RRID:SCR_011848)
    Duplicate reads were then removed using Picard (v2.9.0) [https://broadinstitute.github.io/picard/] and methylation calls were extracted using the methylation caller included as part of Bismark (v. 0.18.1) (98).
    Picard
    suggested: (Picard, RRID:SCR_006525)
    GenPipes (v3.1.6-beta) “DNA Seq’’ pipeline was used to analyze WGS data (97).
    GenPipes
    suggested: (GenPipes, RRID:SCR_016376)
    Seq’’
    suggested: None
    WGS
    suggested: None
    Alignment refinement was done using GATK (v3.8) and include marking duplicate, indel realignment and base recalibration process.
    GATK
    suggested: (GATK, RRID:SCR_001876)
    Methylation raw data were smoothed using a 500-bp smoothing window as a parameter of the BSmooth function implemented in the bsseq Bioconductor package (v. 1.20.0) (100).
    bsseq
    suggested: (bsseq, RRID:SCR_001072)
    Differential methylation was detected using the R Bioconductor package DSS (v2.32.0) adjusting for patient age and sex (101).
    Bioconductor
    suggested: (Bioconductor, RRID:SCR_006442)
    Finally, Distal DMRs were selected by overlapping Promoter DMRs with the GeneHancer v5 (2019) dataset using mergeByOverlaps from IRanges v2.22.2 package, and annotated based on the associated gene (90).
    IRanges
    suggested: (IRanges, RRID:SCR_006420)
    GO Biological Process and KEGG enrichments were detected using the clusterProfiler v3.16.0 package and Reactome enrichment using the reactomePA v1.32.0 package from R.
    GO Biological
    suggested: None
    KEGG
    suggested: (KEGG, RRID:SCR_012773)
    Pathways and GO terms detected by more than one assay with at least 10 genes and FDR < 5% (scRNA and ATAC) or 5 genes and FDR < 10% (DNA methylation) were considered for downstream drug discovery.
    ATAC
    suggested: (Atac, RRID:SCR_015980)
    Due to the hierarchical nature of GO terms and Reactome datasets we selected one term or pathway per hierarchical branch to reduce redundancies for plotting only (full datasets are provided in Suppl.
    Reactome
    suggested: (Reactome, RRID:SCR_003485)
    To identify drugs interacting with genes and hubs in Reactome pathways we used ReactomeFIViz to access the DrugCentral database in Cytoscape (102–104).
    DrugCentral
    suggested: (DrugCentral, RRID:SCR_015663)
    Cytoscape
    suggested: (Cytoscape, RRID:SCR_003032)
    We started by extracting the proteins targeted by the candidate drugs in DrugBank.
    DrugBank
    suggested: (DrugBank, RRID:SCR_002700)

    Results from OddPub: Thank you for sharing your code and data.

    Results from LimitationRecognizer: An explicit section about the limitations of the techniques employed in this study was not found. We encourage authors to address study limitations.

    Results from TrialIdentifier: We found the following clinical trial numbers in your paper:

    IdentifierStatusTitle
    NCT04632381RecruitingIntravenous Zotatifin in Adults With Mild or Moderate COVID-…
    NCT04330690RecruitingTreatments for COVID-19: Canadian Arm of the SOLIDARITY Tria…

    Results from Barzooka: We did not find any issues relating to the usage of bar graphs.

    Results from JetFighter: We did not find any issues relating to colormaps.

    Results from rtransparent:
    • Thank you for including a conflict of interest statement. Authors are encouraged to include this statement when submitting to a journal.
    • Thank you for including a funding statement. Authors are encouraged to include this statement when submitting to a journal.
    • No protocol registration statement was detected.

    Results from scite Reference Check: We found no unreliable references.

    About SciScore

    SciScore is an automated tool that is designed to assist expert reviewers by finding and presenting formulaic information scattered throughout a paper in a standard, easy to digest format. SciScore checks for the presence and correctness of RRIDs (research resource identifiers), and for rigor criteria such as sex and investigator blinding. For details on the theoretical underpinning of rigor criteria and the tools shown here, including references cited, please follow this link.

  3. Version published to 10.1101/2021.09.14.21262309v1 on medRxiv