{'fraction_non_alphanumeric': 0.03219194598504943, 'fraction_numerical': 0.01716099991962061, 'mean_word_length': 5.097280078412154, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Host transcriptional response to SARS-CoV-2 infection in COVID-19 patients

Background One of the most perplexing aspects of infection with the SARS-CoV-2 virus has been the variable response elicited in its human hosts. Investigating the transcriptional changes in individuals affected by COVID-19 can help understand and predict the degree of illness and guide clinical outcomes in diverse backgrounds. Methods Analysis of host transcriptome variations via RNA sequencing from naso/oropharyngeal swabs of COVID-19 patients. Results We report strong upregulation of the innate immune response, especially type I interferon pathway, upon SARS-CoV-2 infection. Upregulated genes were subjected to a comparative meta-analysis using global datasets to identify a common network of interferon stimulated and viral response genes that mediate the host response and resolution of infection. A large proportion of mis-regulated genes showed a reduction in expression level, suggesting an overall decrease in host mRNA production. Significantly downregulated genes included those encoding olfactory, taste and neuro-sensory receptors. Many pro-inflammatory markers and cytokines were also downregulated or remained unchanged in the COVID-19 patients. Finally, a large number of non-coding RNAs were identified as down-regulated, with a few of the lncRNAs associated with functional roles in directing the response to viral infection. Conclusions SARS-CoV-2 infection results in the robust activation of the body’s innate immunity. Reduction of gene expression is well correlated with the clinical manifestations and symptoms of COVID-19 such as the loss of smell and taste, and myocardial and neurological complications. This study provides a critical dataset of genes that will enhance our understanding of the nature and prognosis of COVID-19.


Introduction
SARS-CoV-2, the coronavirus causing COVID-19 has brought enormous economical, health and human cost to the world. SARS-CoV-2, like MERS-CoV and SARS-CoV that caused epidemics in 2003 and 2012, respectively, is the latest coronavirus to infect humans and can cause severe disease with long term and/or fatal consequences in some individuals [1]. COVID-19 progression is unpredictable and ranges from asymptomatic, mild infection to life-threatening acute respiratory distress syndrome (ARDS) [2]. Mild infections are characterized by symptoms such as fever and cough, and are usually resolved spontaneously in a majority of the infected individuals. However, in severe cases, patients develop ARDS and acute lung injury with damage to alveolar lumen leading to inflammation and pneumonia, and resulting in Intensive Care Unit (ICU) admissions or even death [3]. The cause of such differential outcomes is unknown though older age, presence of other disease comorbidities and the host immune response to the virus are believed to be some of the key factors in disease prognosis.
During a typical infection, the immune system activates two distinct antiviral responses in response to pathogen entry. The first defense relies on the body's innate immunity and is mediated by transcriptional induction of type I and III interferon (IFN) and the subsequent activation of IFNstimulated genes (ISGs) [4]. This then triggers the recruitment and coordination of adaptive immune cells via the secretion of chemokines. In some individuals this process can be associated with an excessive increase in pro-inflammatory cytokines. Finally, sometimes the virus may evade the body's innate immunity, leading to an inadequate or delayed response allowing unrestrained viral replication subsequently leading to hyper inflammatory responses [5].
Recent transcriptomic studies investigating the host response using RNA-sequencing (RNA-seq) have tried to identify a molecular signature of the response to SARS-CoV-2 infection. Single cell RNA-seq performed on human cell lines infected with SARS-CoV-2 has demonstrated that canonical ISGs are broadly induced [6]. The failure to launch an appropriate degree of interferonmediated response despite virus replication has been shown to be associated with high levels of chemokine induction in some individuals [7]. The reasons underlying the diversity in the immune response and the correlates for severe infection and death due to COVID-19 are still unclear and need to be investigated across diverse ethnicities, age groups and susceptibilities.
In this work, we have identified the changes induced in the transcriptome of individuals positive for COVID-19 infection. Our analysis of samples from Indian patients reveals a massive downregulation of host genes in response to SARS-CoV-2 infection. The up-regulated genes comprised only ~27% of the mis-expressed genes, most of which were part of a strong immune response and innate immunity network. Down-regulated genes were enriched in pathways that are associated with systemic symptoms of COVID-19 such as cardiac, pancreatic and muscle related issues as well as neuronal and sensory abnormalities. We also observed significant downregulation of smell and taste receptor genes, which correlates well with olfactory dysfunction, now acknowledged as a key symptom of COVID-19. Interestingly, many long noncoding RNA (lncRNAs) were also found mis-regulated, with functional roles in a myriad of functions including mounting a strong antiviral response [8].

Extensive downregulation of host transcriptome in COVID-19 patients
We analyzed total RNA isolated from nasopharyngeal samples of 36 COVID-19 positive patients   from the southern state of Telangana, India, along with 5 COVID-19 negative samples as controls   (see Methods and Supplementary table 1). Principal component analysis ( Figure 1A) shows that the patient samples did not cluster very well together, suggesting high sample variability and possibly a diverse transcriptomic response by individuals against SARS-CoV-2. However, there was a clear distinction between the positive and the control samples and a total of 56,640 genes were analyzed for differential expression (methods).
While most of the up-regulated DEGs were protein coding (220 genes), only 35% (3252 genes) of the down-regulated DEGs were protein coding (Figure 1B and Supplementary table 3). The levels of expression of all the 3472 differentially regulated protein coding genes are summarized across the samples in Figure 1C. A large proportion of the down-regulated DEGs were non-coding RNAs, lncRNAs being the biggest contributors (4018 lncRNAs). The large number of differentially expressed lncRNAs suggests a key role for them during viral infection as described in earlier studies [8,9].
The downregulation of a large proportion of the DEGs is an interesting observation that may be associated with the phenomenon of viral control and switching off host protein production [10].
Besides affecting the immune response, the host cell proteins fade and viral proteins become more prominent with the progress of the infection. The precise mechanisms remain unclear but host shut-off is likely brought about at multiple levels, including down-regulation of host transcription and at the level of translation, where the mis-regulated lncRNAs may also play a role.

Up-regulation of immune related genes and the interferon-mediated immune response
We found that protein coding up-regulated genes were enriched mostly for immune system related GO terms and pathways, especially those associated with Type I interferon signaling and defense response to viral infections ( Figure 2A). Prominent clusters of genes associated with multiple viral infections could also be identified ( Figure 2B We performed a meta-analysis of published datasets (see methods) to identify a set of 19 upregulated genes that showed increased expression in our analysis as well as in other studies on COVID-19 samples (Supplementary Table 6). Protein-protein interaction (PPI) network analysis showed that most of these genes are highly interconnected ( Figure 2C), involved in "innate immune response" (red nodes), the "type I interferon signaling pathway" (green nodes), and/or "defense response to virus" (blue nodes). Interferon (IFN)-induced antiviral proteins play a crucial role in distinguishing self and non-self mRNAs during viral infection by inhibiting viral mRNA expression which lacks 2'-O-methylation of the 5' cap. IFI6, recently proposed to be a part of a novel transcriptional signature for COVID-19 infected cells [13], formed one of the key PPI nodes along with IFIT family gene paralogs that restrict viral translation and modulate viral response pathways [14]. IFIH1 and ISG15 that drive the innate immune response upon sensing viral RNA were also identified as key upregulated genes [15]. The induction of these interferon stimulated genes (ISGs) appears concomitant with the upregulation of STAT1 in these patients. STAT1 binding to ISGs mediates the interferon triggered host response and its dysfunction has been associated with hyperactivation of inflammatory pathways in individuals with acute COVID-19 pathophysiology [16]. IFN mediated activation of the JAK-STAT signaling pathway also has a role in inducing necroptosis, found enriched in our pathway analysis and associated with ARDS development in some patients [17].
We identified a few other genes involved in viral response that were upregulated in patients with active COVID-19 infection. Increased OAS1 levels confer protection against severe COVID-19 [18] and this family of antiviral genes has been shown to be associated with the host response [19]. However, many of the proinflammatory markers (TNF, IL6, IL1B, CCL3, CCL4 and CXCL2) were not significantly changed, suggesting an absence of hyperinflammation in these patients. Overall, our results suggest the activation of a robust innate immune response in COVID-19 patients.

Down-regulated protein coding genes are associated with systemic and cellular host processes
Many of the down-regulated protein coding genes in our analysis were associated with signaling and ion transport related GO terms. We found significant enrichment for neurotransmission, heart and muscular contraction related terms ( Figure 3A), indicating impairment of these cellular processes in the COVID-19 patients. KEGG pathway analysis ( Figure 3B) confirmed the association of the downregulated genes in calcium signaling and the RAS and cAMP signaling pathways, which may help drive viral replication. The RAS signaling pathway may be involved in pathogenesis of COVID-19 and plays an important role in cardiovascular diseases, neurodegenerative diseases, and acute lung injury [26]. We also noticed enrichment of a group of cardiac function related pathways ( Figure 3B). Many of the genes from these KEGG pathways were related to cellular calcium signaling (CACNs) and were down-regulated. The key cardiac proteins, troponin and tropomyosin, were also down-regulated. Calcium ions are required for proper functioning of troponin and tropomyosin dependent cardiac muscle contraction. A disruption of cellular calcium signaling together with down-regulated cardiac proteins indicates possibility of myocardial issues in COVID-19 patients. A few of the down-regulated genes were associated with the pancreatic and insulin secretory systems, in agreement with recent work showing that the insulin requirement for patients with diabetes mellitus is considerably higher at the peak of COVID-19 illness [27].
Interestingly, the KEGG pathways "olfactory transduction" and "taste transduction" were found significantly enriched ( Figure 3B). 105 out of 118 genes from the "olfactory transduction" pathway, were identified as down-regulated, along with several peripheral sensory receptors, especially the olfactory and taste receptors. With better understanding of the disease over the last year, olfactory dysfunction has emerged as a key symptom of COVID-19 [28] and our results suggest sensory impairment, associated with a loss of smell and taste [29], in these patients. There may also be implications for neuronal infectivity via the olfactory and respiratory tracts. The patients of COVID-19 due to anti-NMDAR encephalitis [32]. In this context, many of the genes involved in synapse organization and regulation were also downregulated in our analysis highlighting the potential for neurological complications observed in COVID-19 patients [33].

Systemic downregulation of long non-coding RNAs upon SARS-CoV-2 infection
We identified 4025 lncRNAs mis-regulated in COVID-19, most of them expressed to a lesser extent than the protein coding genes (Supplementary figure S1). A few of these lncRNAs are known to have functional roles during viral infection. For example, ZBTB11-AS1, an anti-sense lncRNA to ZBTB11, that regulates neutrophil development was up-regulated [34]. Anti-sense lncRNA are associated with regulation of transcriptional gene activation triggered by small RNAs [35], and we found that indeed both ZBTB11-AS1 and ZBTB11 were significantly up-regulated in COVID-19 patients (log2FC 2.01 and 1.05, respectively). A few other lncRNAs such as HEIH and

IGF2-AS (associated with recurrence in Hepatitis C virus related hepatocellular carcinoma and
Hep C viral replication) [36,37], PACER (positively regulates COX2 production and NF-kB [38]) and MALAT1 (reduces activation of NF-kB, MAPK and inflammatory factors) were also significantly mis-regulated.
However, the functional role of many lncRNAs in response against viral infections is not well characterized. We therefore asked if any associated protein-coding genes were affected by the lncRNA mis-regulation. Since a substantial proportion of the functionally annotated lncRNAs are cis-acting [39], we searched for the nearest protein-coding genes to the mis-regulated lncRNAs (methods) and identified 720 differentially expressed protein coding genes. We then grouped the lncRNAs based on distance and strand information in context of the nearest protein-coding genes ( Figure 4A). Interestingly, most of the lncRNAs overlap their nearest genes and are on the antisense strand. Further, Gene Ontology analysis showed enrichment of terms only in the groups where the lncRNA was located within the gene or in its upstream promoter ( Figure 4B). These GO terms were associated mostly with brain or neuron development related processes (Supplementary   table 9). Further investigation is needed into the role of the identified lncRNAs in these pathways and their possible association with the complications seen in COVID-19 patients.

Conclusion
This study identifies the key enriched pathways and processes associated with COVID-19 (summarized in Figure 5

Sample metadata
We received 41 samples where 36 samples were confirmed to be COVID-19 positive using RT-PCR method and 5 samples were negative. The age of the patients ranges from 10 to 80 years old.
There were 11 females and 30 males in our dataset. We also had severity of COVID-19 for 30 positive patients. 7 patients required intensive care unit (ICU) intervention while 23 were discharged from COVID-19 ward (W). The detailed metadata information about the samples is provided in Supplementary table 1.

Sample preparation and sequencing
This study was approved by the Institutional Ethics Committee to use patient samples for sequencing and all guidelines were followed. The samples were collected from nasopharyngeal or oropharyngeal swabs as previously described [40]. Total RNA was extracted from 1ml of VTM

Data Processing and Analysis
Raw sequencing reads were checked for quality using FASTQC v0.11.9 [41] and illumina universal adapters and low-quality reads were removed using cutadapt 2.8 [42]. We removed reads with quality scores less than 20 and discarded reads smaller than 36 bp. The processed reads were then again checked for quality using FASTQC and then mapped to the human genome GRCh38 using STAR 2.7.3a with default parameters [43]. The STAR index was generated using the GRCh38 genome and annotation gtf file downloaded from ENSEMBL [44]. BAM files generated were then sorted according to position using Samtools 1.10 [45]. The quality of BAM was assessed using QualiMap v.2.2.2-dev [46]. All the quality reports were compiled using multiqc 1.9 [47].
Uniquely aligned reads were counted using HTSeq [48]. There were 60683 genes in the gtf file for which we had count information. Genes with read count 0 across all the samples were removed resulting in 56640 genes for further analysis. Differential gene expression analysis was performed using DESeq2_1.24.0 [49]. Genes with adjusted p-value < 0.05 and absolute log2 Fold change > 1 were considered differentially expressed resulting in 9319 genes, remaining 47321 genes were removed from further analysis. For PCA plot and heat map, the raw read counts were rlog normalized, available with the DESeq2 package.

Functional enrichment analysis
For functional enrichment analysis, clusterProfiler_3.12.0 [50]was used for GO term enrichment and ShinyGO [51] was used for Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. We only used the Biological process for GO term enrichment analysis. Similar enriched terms were further merged using the 'simplify' function of clusterProfiler with similarity cutoff set to 0.7. 'p.adjust' was used as a feature to select representative terms and 'min' was used to select features. 'Wang' was used as a method to measure similarity. For ShinyGO the web tool was run with the default parameter. The complete list of enriched KEGG pathways were downloaded. For better clarity, only the top 20 significant pathways based on FDR were shown in the figures. Edge cut-off was set to 0.3, so edges were created only if two nodes shared at least 30% of the genes. Darker nodes represent more significantly enriched pathways and size of the node represents the number of input genes in that pathway. Thickness of edges was proportional to percentage of overlapping genes.

Protein-protein interaction network
Protein-protein interaction network was created using STRING-db [52]. Both physical and functional protein associations were considered for edges. All the active interaction sources were considered for generating the network. The medium confidence required for interaction score was set to 0.7. Only input proteins were used for generating the network. Few chosen enriched biological processes or KEGG pathways as reported by STRING-db were shown by color code of the nodes.

Identifying nearest genes for the lncRNA
We used bedtools (v2.26.0) [53] to find nearest genes for the lncRNA. Bedtools closest searches for overlapping features in two bed files. If there is no overlap, the nearest feature is reported with genomic distance between the two features. Closest features were reported irrespective to strand.
The bed file containing differentially expressed lncRNA genomic locations was used as a query bed file (option -a) and a bed file containing all protein coding RNA genomic locations was used as reference bed file (option -b). The lncRNA were classified based on the closest gene. If both the gene and lncRNA were on the same strand, lncRNA was classified as "sense", otherwise it was classified as "anti-sense". Next, the lncRNA was classified based on the distance from the closest gene. If the gene and lncRNA were overlapping, lncRNA was classified as "Gene". If the lncRNA is within 1kb upstream from the closest gene, it is classified as "upstream promoter". Finally, lncRNA upstream of the closest gene by more than 1 kb was classified as "upstream' and lncRNA downstream of the gene was classified as "downstream". This resulted in 8 sets of classification for the lncRNA.

Meta-analysis of existing transcriptomic data for COVID-19
We performed data mining and compiled a list of 9 published works on transcriptomic host response to COVID-19 [54][55][56][57][58][59][60][61][62]. We then overlapped our up-regulated/down-regulated genes with a list of genes that were reported to be up-regulated/ down-regulated in the published datasets. For down-regulated genes the overlap was not good between the published dataset. For up-regulated genes, we generated a list of genes that were found to be up-regulated in at least 3 publications.
These 19 genes are given in Supplementary table 6.

Data availability
All the raw RNA-seq data and their respective count data has been submitted to Gene Expression Omnibus (GEO) database and can be accessed under the accession number GSE166530.

Conflict of Interest
The authors declare no conflict of interest.  genes are connected by edges whose thickness represents percentage of common genes. Size of the node represents the number of genes in that pathway (ranging from 17 to 4) and shade indicates the significance, darker shades being more significant. (C) Protein-protein interactions of 19 genes showing a tightly connected network involved in "innate immune response" (red nodes), "defense response to virus" (blue nodes) and "type I interferon signalling pathway" (green nodes). These genes were up-regulated during COVID-19 infection in our analysis as well as in other published datasets. Edges depict both functional and physical associations. Edge thickness indicates the confidence in the interaction. All active interaction sources in the STRING database are considered. The minimum interaction score for an edge is set at a high confidence level of 0.7. Nodes sharing 30% or more genes are connected by edges whose thickness represents percentage of common genes. Size of the node represents the number of genes in that pathway (ranging from 136 to 7) and shade indicates the significance, darker shades being more significant. (C) Protein-protein interactions of 55 genes showing three tightly connected clusters of genes involved in "cognition" (red nodes), "learning and memory" (blue nodes) and "sensory perception" (green nodes). These genes were involved in 4 addiction related pathways. Edges depict both functional and physical associations. Edge thickness indicates the confidence in the interaction. All active interaction sources in the STRING database are considered. The minimum interaction score for an edge is set at a high confidence level of 0.7. The lncRNAs are divided into 8 groups based on their location -genomic strand and distance to their closest differentially expressed protein coding gene. Number of lncRNA and their closest protein coding genes for each category is shown in the table. Most of the lncRNA were on the anti-sense strand and overlapping with genic region or promoter of the closest protein coding gene. (B) GO enrichment analysis of the nearest protein coding genes associated with the differentially expressed lncRNAs. Enriched terms were identified only when the associated lncRNA was present in the genic or promoter contexts, "Gene antisense" and "Upstream promoter antisense". Only top 10 enriched terms, based on adjusted p-value (p.adj) values, are plotted. Size of dots represents the number of genes and color signifies p.adj value.

Figure 5 : Enriched GO processes and KEGG pathways associated with significant DEGs
Volcano plot showing rlog transformed expression values of all the genes. Each gene is plotted based on the log2 fold change value (X-axis) and -log10 adjusted p-value (Yaxis). Vertical dashed line (absolute (log2 fold change) = 1) and horizontal dashed line (padj = 0.05) shows the criteria set for defining significant DEGs. Genes not changing significantly are colored grey. Up-regulated genes (colored red) are on the right side of the plot while down-regulated genes (colored blue) are on the left side of the plot. Genes from significantly enriched GO terms and KEGG pathways are highlighted. For example, genes labelled as "Interferon related" are associated with the GO terms "type I interferon signaling pathway" and "cellular response to type I interferon". Genes labelled as "Response to virus" are associated with the GO terms "defense response to virus" and "response to virus". Genes marked as "Addiction related" are from 4 enriched KEGG pathways namely "Nicotine addiction", "Morphine addiction", "Cocaine addiction" and "Amphetamine addiction". The other labels correspond to genes from the KEGG pathway "Insulin secretion" and "Olfactory transduction".
Supplementary figure S1: Boxplot showing distribution of expression data for proteincoding genes, lncRNAs and pseudogenes. Distribution of expression value (rlog transformed) for pseudogenes, lncRNAs and protein-coding genes displayed in the form of boxplot . All the pseudogenes, lncRNAs and protein-coding genes are shown on the left side of the plot while differentially expressed are shown on the right side of the plot. Median of the distribution is represented by a horizontal line and whiskers are extended to 1.5 times interquartile range.