key: cord-0786067-hb8wnfim
authors: Szachnowski, Ugo; Bhargava, Anvita; Chazal, Maxime; Foretek, Dominika; Aicher, Sophie-Marie; da Fonseca, Juliana Pipoli; Jeannin, Patricia; Beauclair, Guillaume; Monot, Marc; Morillon, Antonin; Jouvenet, Nolwenn
title: Transcriptomic landscapes of SARS-CoV-2-infected and bystander lung cells reveal a selective upregulation of NF-κB-dependent coding and non-coding proviral transcripts
date: 2022-02-26
journal: bioRxiv
DOI: 10.1101/2022.02.25.481978
sha: cb83b6aaebd2f16bb3837393970c20e6fa5ed59e
doc_id: 786067
cord_uid: hb8wnfim

Detailed knowledge of cellular networks that are modulated by Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is needed to understand viral replication and host response. So far, transcriptomic analyses of interactions between SARS-CoV-2 and cells were performed on mixed populations of infected and uninfected cells or using single-cell RNA sequencing, both leading to inaccurate or low-resolution gene expression interpretations. Moreover, they generally focused on annotated messenger RNAs (mRNAs), ignoring other transcripts, such as long non-coding RNAs (lncRNAs) and unannotated RNAs. Here, we performed deep polyA+ transcriptome analyses of lung epithelial A549 cells infected with SARS-CoV-2, which were sorted based on the expression of the viral protein spike (S). To increase the sequencing depth and improve the robustness of the analysis, the samples were depleted of viral transcripts. Infection caused a massive reduction in mRNAs and lncRNAs, including transcripts coding for antiviral innate immune proteins, such as interferons (IFNs). This absence of IFN response probably explains the poor transcriptomic response of bystander cells co-cultured with spike positive (S+) ones. NF-κB and inflammatory response were among the pathways that escaped the global shutoff in S+ cells. In agreement with the RNA-seq analysis, inflammatory cytokines, but not IFNs, were produced and secreted by infected cells. Functional investigations revealed the proviral function of the NF-kB subunit p105/p50 and some of its known target genes, including IL32 and IL8, as well as the lncRNA ADIRF-AS1, which we identified as a novel NF-kB target gene. Thus, analyzing the polyA+ transcriptome of sorted populations of infected lung cells allowed unprecedented identification of cellular functions that are directly affected by infection and the recovery of coding and non-coding genes that contribute to SARS-CoV-2 replication.

to dissect virus-cell interactions.

Here, we investigated the coding and non-coding transcriptional landscape of lung cells infected 116 with SARS-CoV-2 and sorted according to the expression of the viral protein spike (S). Our deep 117 transcriptome analysis using annotated RNA genes and reference-based RNA profiler uncovered 118 pathways that are directly affected by infection and identified coding and non-coding genes contributing 119 to an optimal SARS-CoV-2 replication. 

PolyA+ RNAs were isolated from mock-infected, S + and Scells. Around 85% of the total reads mapped 133 to the viral genome in S + cells, while less than 5% of the total reads aligned with the viral genome in S -134 cells (Fig. 1C) , validating our sorting approach. The large dominance of viral reads over cellular reads 135 illustrates the ability of the virus to hijack the cellular machinery for its replication. A similar proportion 136 of SARS-CoV-2 reads in the RNA pool was previously reported in lung epithelial carcinoma Calu-3 137 cells infected for 8 hours [20] . These differences in the representation of viral RNA between S + and S -138 5 cells altered the robustness of the statistical analysis used to identify DEGs. To overcome this limitation, 139 the samples were depleted of viral RNAs (vRNAs) using a set of oligonucleotide probes covering the 140 entire viral genome (Fig. 1A) . Following depletion, viral reads represented between 0,01 and 2,8% of 141 the total reads both in S + and Scells (Fig. 1C ).

Coding and long non-coding genes were identified using gencode annotation (v32), while 143 unannotated RNAs were recovered with Scallop assembler [36] (Fig. 1D-F, Fig. S1 and tables S1-S3).

Principal component analysis (PCA) of polyA+ transcriptomes segregated S + cells from Sand mock-145 infected ones (Fig. 1D ). This segregation based on S expression represented around 92% of the 146 transcriptomic differences between the samples (Fig. 1D ). Only subtle differences (2,5%) distinguished 147 bystander and mock-infected cells (Fig. 1D) , suggesting that the transcriptional landscapes of these 2 148 cell populations were very similar. An absence of response of Scells was unexpected since cytokines, 149 which are commonly secreted by virally infected cells, activate an antiviral state in bystander cells 150 through surface receptors.

Analysis of gene expression allowed identification of thousands of annotated coding and non-152 coding genes that were differentially expressed (absolute fold change ≥2, p-value < 0.05) in S + cells as 153 compared to Sor mock-infected ones (Fig. 1E -F and tables S1-S3). We identified around 13 times more 154 downregulated coding genes than upregulated ones in S + cells ( Fig. 1E-F) , suggesting that infection 155 triggers a massive, but incomplete, shutoff of gene expression. Among the top upregulated coding genes 156 in S + cells, we confirmed candidates revealed by previous analyses performed in non-sorted non-vRNA-157 depleted A549-ACE2 cells, such as CXCL8, CCL20, IL6 and NFKB1 [19, 38, 39] , but also novel highly 158 significant candidates, including IL32 and ITGAM (table S1). The genes encoding IFN type I and type 159 III were not significantly upregulated in S + cells, as compared to mock-infected cells. Accordingly, ISGs 

In agreement with the PCA (Fig. 1D) , volcano plots and heat maps revealed that Sbystander cells 175 and mock-infected control cells exhibited very similar transcriptomic profiles (Fig. 1E-1F and Fig. S1A -176 S1B). Only around 170 polyA+ transcripts were differentially expressed in Scells as compared to 177 mock-infected ones ( Fig. 1E and tables S1-S3). As a comparison, over 13000 DEGs were identified in 178 S + as compared to mock-infected cells ( Fig. 1E and tables S1-S3). These analyses further suggest that 179 S+ cells present none or very little paracrine signaling response. Among the 69 coding genes that were 180 upregulated in Scells as compared to mock-infected, 29 were also upregulated in S + cells (Fig. S1B 181 and table S3). Some of these common genes were inflammatory genes, such as IL32, IL6 and CCL20.

Among the 39 upregulated coding genes that were unique to Scells, 16 were ISGs (examples include 183 MX1, APOL1 and IFI6). The expression of these inflammatory genes and ISGs in bystander cells could 184 be induced early in infection, prior to the production of S proteins.

Our approach reveals that SARS-CoV-2 infection triggers a major shutoff of gene expression in 186 A549-ACE2 cells. It also shows that Scells do not exhibit a strong transcriptional signature despite 187 being cultured with S + cells, suggesting the absence of an efficient paracrine communication.

To compare our differential deep analyses with known datasets, we analyze publicly available 191 polyA+ RNA-seq raw data of unsorted A549-ACE2 infected with SARS-CoV-2 at a MOI of 0.2 [19] .

Viral reads represented around 50% of the total number of reads in these unsorted bulk population of 193 cells [19] , which was expectedly less than in A549-ACE2 cells positive for S (Fig. 1C) . The 2 analyses 194 shared 150 upregulated protein-coding genes and 238 downregulated ones ( Fig. 2A, 

S + ) supports this hypothesis (Fig. S1C ). About 41% of the upregulated protein-coding genes and 16% 203 of downregulated ones that we identified did not appear in conventional RNA-seq analysis of mixed 204 populations [19] . This comparison highlights the accuracy and the depth of our analysis.

To validate the sorting approach combined with vRNA-depletion, we compared mRNA 206 abundances of a few DEGs in a bulk population of cells infected with SARS-CoV-2 for 24 hours, as 207 well as in sorted S + and bystander Scells infected in the same condition. As expected, S+ cells produced 208 approximately 200-fold more intracellular viral RNAs than did Scells (Fig. 2B) . These qPCR analyses 209 confirm that some Scells are at an early stage of viral replication, prior to viral protein expression (Fig. 

1B). We included in the analysis three coding transcripts (IL32, ITGAM and TRAF1), two lncRNAs 7 (WAKMAR2 and AL132990.1) and one unannotated transcript (XLOC_007519) that were identified 212 amongst the most upregulated RNAs in S + A549-ACE2 cells (Fig. 2C , S2A and tables S1-S3). The 213 abundance of IL32 mRNA did not increase significantly in the infected bulk population, as compared 214 to mock-infected cells (Fig. 2C ). By contrast, IL32 mRNA levels increased around 30-fold in S + cells, 215 compared to those in mock-infected cells (Fig. 2C) . This difference explains why IL32 was not identified 216 as an up-regulated gene in previous RNA-seq analysis performed in mixed population of infected A549-217 ACE2 cells [19, 38] . Similarly, the expression of ITGAM, TRAF1, WAKMAR2, AL132990.1 and 218 XLOC_007519 showed a modest increase of mRNA abundances in the bulk population and a significant 219 increase in S + cells, as compared to mock-infected cells ( Fig. 2C and S2B ). The decreased expression 220 of transcripts identified as top downregulated hits in the RNA-seq analysis of S+ cells, such as the coding 221 transcripts FEN1 and SNRPF, the lncRNAs AC016747.1, DANCR and TP53TG1, as well as the 222 unannotated RNA XLOC_049236 (Fig. S2A ), was significantly more pronounced in S + cells than in the 223 mixed population of cells, when compared to mock-infected cells ( Fig. 2D and S2C ). Analysis of RNA 224 abundances in sorted cells thus highlighted the increased accuracy of our approach, compared to 225 classical methods, in detecting up-and down-regulated genes.

To identify pathways affected by infection in A549-ACE2 cells, we performed Gene Ontology 227 (GO) terms and KEGG pathway enrichment analysis on the upregulated coding genes in S + cells, as 228 compared to mock-infected cells (Fig. 2E ). We observed a significant enrichment in several 229 inflammatory signaling pathways, including TNF and NF-κB signatures, which were previously 

To investigate these possibilities, we selected 5 inflammatory chemokines (IL-6, CXCL1, CCL2, 247 CXCL8/IL-8 and CCL20) whose expression was upregulated in S + cells upon infection (table S1) and 8 quantified their intracellular and secreted levels in lysates and supernatants of A549-ACE2 cells infected 249 for 24 hours (Fig. 3) . As a comparison, A549-ACE2 cells transfected with the immuno-stimulant 250 poly(I:C) were included in the analysis. In a mixed population of S + and Scells, mRNAs of these 5 251 cytokines were significantly more abundant than in mock-infected cells (Fig. S3A ), in agreement with 252 the increased levels of mRNAs detected in S + cells by RNA-seq (table S1). Their expression was also 253 induced by poly(I:C) (Fig. S3A ). All five cytokines were expressed at detectable levels in cells as little as 50 pg/ml of IFNβ, which was similar to the quantity secreted by mock-infected cells and 283 lipofectamine-exposed cells, likely representing baseline levels (Fig. 3D ). MeV infected cells secreted 284 around 1000 pg/ml of IFN-λ1 and 5000 pg/ml of IFN-λ2/3 while no IFN-λ was detected in the 9 supernatant of SARS-CoV-2 infected cells (Fig. 3D ). This baseline level of IFN type-I secretion and 286 absence of IFN type-III release by SARS-CoV-2-infected cells is likely to be responsible for the lack of 287 paracrine signaling revealed by the RNA-seq analysis (Fig. 1 ).

Upregulated NF-κB target genes contribute to an optimal SARS-CoV-2 replication 290 Numerous genes associated with the NF-κB signaling pathway fall into the category of genes that 291 escaped the virus-induced cellular shutoff (Fig. 2E ). To determine which of these genes were directly 292 controlled by NF-κB, we cross-compared the upregulated genes with known NF-κB target genes.

Among the 68 upregulated NF-κB-targets in S + cells, we identified cytokines such as CXCL8/IL8 and 294 IL32 (Fig. 4A , S4A and table S5). NFKB1, which codes for the p105/p50 subunit of the transcription 295 factor, and is itself a NF-κB-target gene [49, 50] , also showed a significant transcriptional induction in 

These results confirmed the pro-SARS-CoV-2 activity of NFKB1 in A459-ACE2 cells [38] and revealed 321 that CXCL8, IL32 and ADIRF-AS1 also exhibited significant proviral functions.

10 Thus, our sorting approaches identified coding and non-coding genes that contribute to an optimal 323 SARS-CoV-2 replication. 

Although our RNA-seq analysis identified over 12000 host transcripts that were significantly 370 reduced during SARS-CoV-2 infection as compared to control cells, it also recovered around 1500 371 transcripts whose levels were significantly elevated and 2800 transcripts whose levels were unchanged 372 upon infection. Among top upregulated genes in S + cells, we identified numerous proinflammatory 373 cytokines, such as IL6, CXCL1, CCL2, IL8/CXCL8 and CCL20. ELISA analysis confirmed that 374 infected cells were producing these inflammatory cytokines. They were previously identified in bulk or 375 sc-RNA analysis of A549-ACE2 cells as upregulated [19, 38, 39] , while others, such as IL32, were 

Our data suggests that the genes that are refractory to the viral-induced shutoff are proviral genes.

Understanding the molecular mechanisms underlying the selectivity of the shut-off would be interesting.

Since coronavirus Nsp1 induces the cleavage of the 5'UTR of capped transcripts bound to 40S 441 ribosomes, the 5'UTR length and/or structure may affect Nsp1 binding and subsequent degradation.

Alternatively, the extent of transcript reduction may be linked to their GC content and/or their lengths, 443 which could affect the specificity of the host RNase that is presumably recruited by Nsp1. Discovering and bystander vs infected), and genes were retained as differential if adjusted p-value was < 0.05 and 520 log fold-change > 1 or < -1. All plots were made using custom script, except for heatmaps that were 521 done using pheatmap package (RRID:SCR_016418).

Fastq files produced in the study of Blanco-Melo et al (2020) 

Differential analysis was performed as described above.

GO enrichment analysis. The GO enrichment and KEGG pathway analysis were performed using 536 DAVID online tool (updated version 2021) [101, 102] . Upregulated protein-coding genes from each 537 comparison were taken for the analysis with default background for Homo sapiens.

GOTERM_BP_DIRECT and KEGG pathway were retained and top 10 results based on adjusted p-539 value (Benjamini) were plotted using ggplot2 R package (v 3.3.0).

Identification of NF-κB target genes. A list of coding genes that are known targets of NF-kB is 541 available on Gilmore's laboratory website (https://www.bu.edu/nf-kb/gene-resources/target-genes). We 542 selected genes from this list that were shown to be direct targets of NF-κB, and for which the gene 543 symbol could be retrieved in gencode annotation (354 genes). For identifying lncRNAs and 544 16 unreferenced RNAs that possess NF-κB binding site in their promoter, we used p65 ChIP-seq data from 545 GEO dataset GSE34329 [51] -one input file and 2 ChIP replicates, 38nt long reads, single-end. Reads 546 were mapped using bowtie2 v2.4.1 using hg38 as reference, and SAMtools was used to retained the one 547 flagged as primary alignment, with mapping quality > 30, and to remove PCR duplicates (markdup, 548 with -r option). NF-κB binding sites were then detected using macs2 v2. 

Step Kit (Thermo Fisher Scientific). Data were analyzed using the 2-ΔΔCT method, with all samples 563 normalized to endogenous BPTF, whose gene expression was confirmed as homogenous across samples Figure S4 . (A) Visualization of read coverage (tag/nucleotide) from polyA+ RNA-seq and ChIP-seq (IP -Input) at NFKB1, CXCL8, IL32 and ADIRF-AS1 loci. RNA-seq and ChIP-seq data were normalized independently, on ERCC reads for RNA-seq and on library size for ChIP-seq. (B) RT-qPCR quantification of NFKB1, CXCL8 and ADIRF-AS1 transcripts induction in A549-ACE2 cells, 24 hours post infection with SARS-CoV-2 (MOI of 1) (normalized fold change over mock-infected, n=3 independent experiments, ratio-paired t test, line at mean ± SEM).

COVID-19): A 610 Review

On the 612 whereabouts of SARS-CoV-2 in the human body: A systematic review

CoV-2: a storm is raging

Clinical and immunological 617 features of severe and moderate coronavirus disease 2019

Longitudinal 620 analyses reveal immunological misfiring in severe COVID-19

Complex Immune Dysregulation in COVID-19 Patients with Severe 624 Respiratory Failure

Inflammatory responses and 627 inflammation-associated diseases in organs

Clash of the titans: interferons and SARS-CoV-2

Stimulation of Innate Immunity by Host and Viral RNAs. 632

Interferon-Stimulated Genes: A Complex 634 Web of Host Defenses

Recent advances in antiviral interferon-stimulated gene biology

Viruses and interferon: a fight for supremacy

Evasion of Type I 641 Interferon by SARS-CoV-2

Mechanisms of Antiviral Immune Evasion of SARS-CoV-2

Interplay between SARS-CoV-2 and 645 the type I interferon response

A virus-derived 648 microRNA targets immune response genes during SARS-CoV-2 infection

SARS-CoV-2 651 expresses a microRNA-like small RNA able to selectively repress host genes

Impaired type 654 I interferon activity and inflammatory responses in severe COVID-19 patients

Imbalanced Host Response to SARS-CoV-2 Drives Development of COVID-19

SARS-CoV-2 660 uses a multipronged strategy to impede host protein synthesis

Transcriptomic 663 profiling of SARS-CoV-2 infected human cell lines identifies HSP90 as target for COVID-19 664 therapy. iScience

Transcriptomic 666 characteristics of bronchoalveolar lavage fluid and peripheral blood mononuclear cells in 667 COVID-19 patients

Deconvolution of pro-and antiviral genomic responses in Zika virus-infected and bystander 671 macrophages

Single-cell 674 analyses reveal SARS-CoV-2 interference with intrinsic immune response in the human gut

Single cell resolution of SARS-CoV-2 tropism, antiviral responses, and susceptibility to 678 therapies in primary human airway epithelium

Single-cell longitudinal analysis of SARS-CoV-2 infection in human airway epithelium 682 identifies target cells, alterations in gene expression, and cell state changes

Single-cell landscape of 685 bronchoalveolar immune cells in patients with COVID-19

Single-Cell RNA-Seq Technologies and Related 688 Computational Data Analysis

Determining sequencing depth in a single-cell RNA-seq 690 experiment

Functional Classification and Experimental Dissection of Long 692 Noncoding RNAs

Encode Project Consortium. An integrated encyclopedia of DNA elements in the 697 human genome

Chromatin 699 signature reveals over a thousand highly conserved large non-coding RNAs in mammals. 700

RNAs: Emerging and Versatile Regulators in Host-Virus Interactions

Interferon Regulates the Expression of Long Non-Coding RNAs

Accurate assembly of transcripts through phase-preserving 708 20 graph decomposition

Bridging the gap between reference and real transcriptomes

The NF-712 κB Transcriptional Footprint Is Essential for SARS-CoV-2 Replication

SARS-CoV-2 infection induces a pro-inflammatory cytokine response through cGAS-STING 716 and NF-κB

Noncoding RNA Downregulated in Human Chronic Wounds

Motility and Production of Inflammatory Chemokines

LncRNA HOXA-AS2 represses 722 endothelium inflammation by regulating the activity of NF-κB signaling

Long Noncoding RNA Blocks IκB Phosphorylation and Suppresses Breast Cancer Metastasis

NF-kappaB and virus infection: who controls whom

Signalling pathways of the TNF superfamily: a double-edged sword

CoV-2 Disrupts Splicing, Translation, and Protein Trafficking to Suppress Host Defenses

Dynamic 735 competition between SARS-CoV-2 NSP1 and mRNA on the human ribosome inhibits 736 translation initiation

Nonstructural Protein 1 739 of SARS-CoV-2 Is a Potent Pathogenicity Factor Redirecting Host Protein Synthesis 740 Machinery toward Viral RNA

Measles virus activates 743 NF-kappa B and STAT transcription factors and production of IFN-alpha/beta and IL-6 in the 744 human lung epithelial cell line A549

Modulation of NF-κB-dependent gene transcription using programmable DNA minor groove 753 binders

p50-associated COX-2 extragenic RNA (PACER) 755 activates COX-2 gene expression by occluding repressive NF-κB complexes

SARS-CoV-2 tropism, entry, replication, and propagation: Considerations for drug discovery 759 and development

CoV-2 disrupts the mRNA export machinery to inhibit host gene expression

SARS coronavirus nsp1 protein induces template-dependent endonucleolytic cleavage of 765 mRNAs: viral mRNAs are resistant to nsp1-induced RNA cleavage

A two-pronged 768 strategy to suppress host protein synthesis by SARS coronavirus Nsp1 protein

Severe 771 acute respiratory syndrome coronavirus nsp1 protein suppresses host gene expression by 772 promoting host mRNA degradation

The N-terminal domain of SARS-CoV-2 nsp1 plays key roles in suppression of cellular gene 776 expression and preservation of viral gene expression

Comparative 779 transcriptomic analysis of SARS-CoV-2 infected cell model systems reveals differential 780 innate immune responses

SARS-CoV-2 N protein 782 antagonizes type I interferon signaling by suppressing phosphorylation and nuclear 783 translocation of STAT1 and STAT2

CoV-2 Orf6 hijacks Nup98 to block STAT nuclear import and antagonize interferon 787 signaling

SARS-CoV-2 triggers an MDA-5-dependent interferon response which is unable to control 790 replication in lung epithelial cells

Experimental and natural evidence of SARS-CoV-2-infection-induced activation of type I 793 interferon responses. iScience

MDA5 795 Governs the Innate Immune Response to SARS-CoV-2 in Lung Epithelial Cells

Activation and evasion of type 798 I interferon responses by SARS-CoV-2

Critical 801 Role of Type III Interferon in Controlling SARS-CoV-2 Infection in Human Intestinal 802

An organoid-derived bronchioalveolar model for SARS-CoV-2 infection of human 805 alveolar type II-like cells

Critical 807 Role of Type III Interferon in Controlling SARS-CoV-2 Infection in Human Intestinal

Suppression of NF-κB 810 Activity: A Viral Immune Evasion Mechanism. Viruses

SARS-CoV-2 infection induces a pro-inflammatory cytokine response through cGAS-STING 814 and NF-κB

Activation of NF-κB and induction of proinflammatory 816 cytokine expressions mediated by ORF7a protein of SARS-CoV-2

SARS-CoV-2 Nsp14 819 activates NF-κB signaling and induces IL-8 upregulation

SARS-CoV-2 Nsp5 Activates NF-822 κB Pathway by Upregulating SUMOylation of MAVS

Interleukin-32 in Pathogenesis of Atopic Diseases: 825 Proinflammatory or Anti-Inflammatory Role

Extensive Antiviral Function via Selective Stimulation of Interferon λ1 (IFN-λ1) *

Endogenous IL-32 controls cytokine and HIV-1 production

IL-32: A Host Proinflammatory 834 Factor against Influenza Viral Replication Is Upregulated by Aberrant Epigenetic 835 Modifications during Influenza A Virus Infection

Interleukin-32: a cytokine and 838 inducer of TNFalpha

The α Chemokine, Interleukin 8, Inhibits the Antiviral Action of Interferon α

Hepatitis B virus (HBV) induces the expression of interleukin-8 that in turn reduces HBV 844 sensitivity to interferon-alpha

Negative 846 regulation of the interferon response by an interferon-induced long non-coding RNA

Unique 849 signatures of long noncoding RNA expression in response to virus infection and altered innate 850 immune signaling

NRAV, a long noncoding 852 RNA, modulates antiviral responses through suppression of interferon-stimulated gene 853 transcription

The LPS-inducible lncRNA 855 Mirt2 is a negative regulator of inflammation

The Neat Dance of COVID-19: NEAT1, DANCR, 858 23 and Co-Modulated Cholinergic RNAs Link to Inflammation

Protein Coding and Long Noncoding RNA 864 (lncRNA) Transcriptional Landscape in SARS-CoV-2 Infected Bronchial Epithelial Cells 865 Highlight a Role for Interferon and Inflammatory Response

Frenkel-Morgenstern M. mRNA-lncRNA Co-868

Expression Network Analysis Reveals the Role of lncRNAs in Immune Dysfunction during 869 Severe SARS-CoV-2 Infection

A 871 Molecularly Cloned Schwarz Strain of Measles Virus Vaccine Induces Strong Immune 872 Responses in Macaques and Transgenic Mice

Cutadapt removes adapter sequences from high-throughput sequencing 880 reads

STAR: 882 ultrafast universal RNA-seq aligner

888 deepTools2: a next generation web server for deep-sequencing data analysis

891 Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and 892 isoform switching during cell differentiation

BEDTools: a flexible suite of utilities for comparing genomic 895 features

featureCounts: an efficient general purpose program for 897 assigning sequence reads to genomic features

R: A Language and Environment for Statistical Computing

Moderated estimation of fold change and dispersion for 902 RNA-seq data with DESeq2

Systematic and integrative analysis of large 904 gene lists using DAVID bioinformatics resources

Bioinformatics enrichment tools: paths 907 toward the comprehensive functional analysis of large gene lists