{'fraction_non_alphanumeric': 0.026646037476362384, 'fraction_numerical': 0.008079766202509884, 'mean_word_length': 5.020006898930665, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Transcriptogram analysis reveals relationship between viral titer and gene sets responses during Corona-virus infection.

To understand the difference between benign and severe outcomes after Coronavirus infection, we urgently need ways to clarify and quantify the time course of tissue and immune responses. Here we re-analyze 72-hour time-series microarrays generated in 2013 by Sims and collaborators for SARS-CoV-1 in vitro infection of a human lung epithelial cell line. Transcriptograms, a Bioinformatics tool to analyze genome-wide gene expression data, allow us to define an appropriate context-dependent threshold for mechanistic relevance of gene differential expression. Without knowing in advance which genes are relevant, classical analyses detect every gene with statistically-significant differential expression, leaving us with too many genes and hypotheses to be useful. Using a Transcriptogram-based top-down approach, we identified three major, differentially-expressed gene sets comprising 219 mainly immune-response-related genes. We identified timescales for alterations in mitochondrial activity, signaling and transcription regulation of the innate and adaptive immune systems and their relationship to viral titer. At the individual-gene level, EGR3 was significantly upregulated in infected cells. Similar activation in T-cells and fibroblasts in infected lung could explain the T-cell anergy and eventual fibrosis seen in SARS-CoV-1 infection. The methods can be applied to RNA data sets for SARS-CoV-2 to investigate the origin of differential responses in different tissue types, or due to immune or preexisting conditions or to compare cell culture, organoid culture, animal models, and human-derived samples.

They harvested total RNA in triplicate at the time points described above, and analyzed global mRNA expression. They detected little cellular gene differential expression in comparison with Mock samples during the first 24 ℎ of icSARS-CoV infection. Figure S 1. Upper panel: Relative Transcriptograms (Radius 30) using Mock 0 ℎ as a control for both WT at 0 ℎ and DORF6 (mutant strain that does not express ORF6 protein) at 0 ℎ for expression data downloaded from the Gene Expression Omnibus site. The offset in the raw data is clear. Bottom panel: Transcriptograms produced from expression data re-normalized so that the average expression for each sample is one, eliminating the offset. All analyses in this paper used re-normalized data.

Volcano plots
We produced Volcano plots comparing time-matched Wild Type samples against Mock samples for all 11 time points. As expected, the locations in the plots of the genes in clusters A, B, and C change with the time after first RNA harvest. Before 24 ℎ, most genes in these clusters are not differentially expressed. At and after 24 ℎ, most genes in cluster A and most genes in cluster B localize to the down-regulated branch (left branch), while genes in cluster C mostly localize to the up-regulated branch (right branch). Genes in cluster B follow the trend of cluster C before 48 ℎ, after which they localize in regions with more modest values of fold change. The localization of the genes from different Transcriptogram-identified clusters to different compact regions in the Volcano plots supports the validity of our Transcriptogram clustering. However, without the Transcriptogram analysis, we could not have identified these genes as functionally significant, because they would be hidden among the many other genes that the Volcano plots find are differentially expressed and assign to the same regions.

Figure S 2. Volcano plots (P-value versus Fold-Change) for all genes with probes in the microarrays, comparing timematched WT samples to Mock samples. We highlight the genes belonging to the clusters found with the Transcriptogram method in pink (Cluster A), blue (Cluster B) and green (Cluster C). The Transcriptogram clusters
populate the regions with low P-values and large fold changes in the Volcano plots, supporting the statistical significance of their differential expression. However, the Volcano plots also find many other genes which overlap these clusters, which the Transcriptogram analysis does not find to be biologically relevant.

Correlation between gene expression and virus titer
For all genes with probes in the microarray, we calculated the time correlation between the differential expression of the -th gene and the virus titer, that is: where ( ) is the differential expression defined in Eq.(1) in the main text (the ratio between the average expression of WT samples and Mock samples at time t), � is the time-average of differential expression for the -th gene, and is the corresponding standard deviation. ( ) is the virus titer at time , ̅ is the virus titer time average and the corresponding standard deviation. Figure S 3 shows for all genes, in random order. We have also highlighted the genes belonging to clusters A (pink), B (blue), and C (green). Not surprisingly, genes in cluster B have the highest correlation with virus titer: their differential expression increases and decrease in parallel with changes in the virus titer. Differential expression of genes in cluster C usually has positive correlation with changes in viral titer, but the correlation is weaker than for genes in cluster B, since their differential expression keeps increasing after the virus titer starts to decrease. Finally, differential expression of genes belonging to cluster A correlates negatively with changes in virus titer. These results agree with those found using the Transcriptogram method, but without the Transcriptogram analysis, we could not have identified these functionally relevant genes among the host of other genes whose differential expression has similar net time correlations with changing virus titer.  Positively-and negatively-correlated genes segregate after 24 ℎ, with maximal separation at 30 ℎ and 36 ℎ, with positively-correlated (red) genes lying in the up-regulated branch of the Volcano plots and negatively-correlated (blue) genes lying in the down-regulated branch. At later times, some positively-correlated genes migrate to the down-regulated region of the Volcano plots though the overall separation is maintained. Figure S 4 again shows that classical analysis, by itself, cannot identify biologically meaningful genes.

Figure S 4. Volcano plots (P-value versus Fold-Change) for all genes with probes in the microarrays, comparing timematched WT samples to Mock samples. We have colored each points by the gene's time correlation between its differential expression and the virus titer, following the color bar in the right lower corner of the panel. Before 12 ℎ, genes with positive (red) and negative (blue) temporal correlation with the viral titer occur in both the up-regulated and down-regulated branches of the Volcano plots.
At 30 ℎ and 36 ℎ the positively-correlated genes mostly lie in the up-regulated branch, while the negatively-correlated genes lie in the down-regulated branch. Overall segregation persists but gradually weakens at later times.

t-test for differential gene expression in the clusters
Here we present t-tests for the mean expression of the clusters, comparing Wild Type samples to Mock samples. Although the Transcriptograms suggest that mean differential expression (as defined in Eq. (1) in the main text) of these clusters before 12 ℎ are not relevant, and genes with positive and negative temporal correlation with virus titer are present in both branches of the Volcano plots, t-tests indicate that differential expression of genes in cluster B, calculated between the Wild Type and Mock samples is statistically significant, presenting a low P-value (see Table S1, below). However, while statistically significant, this differential expression is probably not functionally relevant in this experiment. The t-test is overly sensitive and is detecting small variations that are not that significant to the cells' immune response to viral infection. After 24 ℎ, the differential expression between the clusters is both statistically significant and biologically relevant, as indicated by our Transcriptogram analyses and by the segregation of genes positively and negatively correlated with the virus titer in the Volcano plots.

GSEA and Transcriptogram results
The classical analyses confirm the results obtained using Transcriptograms, but would not, on their own be sufficient to produce gene sets A, B and C. Here we present results using a bioinformatic tool specially designed for gene-set analysis, Gene Set Enrichment Analysis (GSEA) [5].
Gene Set Enrichment Analyses (GSEA) investigate differential expression in genomewide microarray or RNA-Seq transcription data. After quality control and normalization, GSEA can test for differential mean expression of gene sets between two classes of samples. Gene sets may be custom-built or use public data bases, such as Gene Ontology terms or KEGG pathways. GSEA first produces a gene list, ranking each gene by one of several possible metrics used to express the differential expression between genes of two classes (WT versus Mock, for example). The default metric is Signal-to-Noise ratio (calculated as the difference between each gene expression averaged over the replicates of each class divided by the sum of the corresponding standard deviations) but "Ratio between Classes" (calculated as the ratio between the gene expression average over replicates for each class) is also an option. The second step locates all genes from a given gene set (for example, the sets of genes belonging to a KEGG metabolic pathway or a term in the Gene Ontology) within the ranked list. If genes from a gene set are predominantly located in high-(low-) scoring regions of the ranked gene list, the gene set is identified as differentially up(down)-expressed. The significance of the identification may be estimated by calculating the False Discovery Rate (FDR) calculated under permutation of either the samples or the genes. For sample permutation, identification is accepted as significant when the FDR < 25%. For gene-list permutation, identification is accepted as significant when the FDR < 5%. When all classes have few replicates (less than 7), the acceptance criterion is based on gene-set permutation.
We ran GSEA analyses for the 11 time points comparing WT and Mock classes, with the following selection of metrics: a. Expression data, Signal-to-Noise metric, gene-set permutation, testing all GO:Biological Function terms. Expression data, Ratio-between-Classes metric, gene-set permutation, testing all GO:Biological Function terms.
Like GSEA, Transcriptograms help identify biologically-meaningful genes in the presence of noisy signals. We applied GSEA to the gene sets identified by the Transcriptograms of radius 30 (as in the original analysis). This workflow considers only genes that are in the Transcriptogram ordered gene list. The Transcriptogram's running-window average reduces signal noise and assign a value to each gene in the ordering. We then ran GSEA for the Transcriptogram values as follows: b. Transcriptogram values (after averaging over neighbors) with radius 30, Ratio-between-Classes metric, gene permutation, for all GO:Biological Function terms.  GSEA results depend dramatically on the choice of metric and data pre-processing. Comparing the two metrics applied to the normalized data, at 24 ℎ after first RNA harvest and after, GSEA with the Signal-to-Noise metric finds hundreds of significantly differentiallyexpressed terms, while GSEA with the Ratio-between-Classes metric finds essentially none. This sensitivity to the metric used is a concern for the reliability of the method. Using the Transcriptogram to smooth the data with a running window average over the ordered gene list decreases noise. Consequently, the number of terms found to be significantly differentiallyexpressed increases. In fact, the ES (enrichment score) for each gene produced by GSEA applying Ratio-between-Classes to Transcriptogram-smoothed data is identical to the relative Transcriptogram presented in Figure 3 in the main text or Figure S 6 below, differing only in the order of genes: the assigned value for each gene is the same. The hundreds of terms from GO:BP which GSEA finds to be differentially expressed are statistically significant in the sense that their difference in expression is much larger than their standard deviation. However, we may expect that many of these statistically-significant terms are not relevant to the biology of the problem (for example, 409 terms were found up-regulated in Mock samples, when applying GSEA to expression data at 3 ℎ after RNA harvest). The Transcriptograms in Figure S 6 show the variation in scale of the differential expression of different gene sets: the alterations at, for example, 7 ℎ after first RNA harvest are small compared to those observed after 36 ℎ. This comparison may be used to set a scale for signal relevance, used to reduce the number of genes to further analyze. Our analysis pipeline sets the scale for relevant alterations to be a 9/5-fold change in the Transcriptogram values, neglecting peaks with small differential Transcriptogram values in comparison to the most altered ones observed after 36 ℎ. Naturally, when increasing the threshold required to select gene sets to be further analyzed, and hence decrease the number of false-positive errors, we will simultaneously increase the incidence of false-negative errors. We stress that the present analysis tries to address problems caused by the abundance of false-positive errors in finding biologically-relevant genes in experiments showing genome-wide transcription alterations.
Finally, we tested the gene sets in Clusters A, B, and C using GSEA using the GSEA default options (Signal-to-Noise, gene-set permutations). GSEA found these sets to be significantly differentially expressed after 30 ℎ after first RNA harvest (see also Table S 3 below), in agreement with our Transcriptogram results. GSEA defines a normalized enrichment score (NES) by normalizing ES (here the ratio between Transcriptogram values for WT and for Mock samples) by the expected value for a gene set with a given number of genes. Figure S 5 shows the evolution for each cluster of the average NES for each cluster follows the same trends as in Figure 5 in the main text.  Table S 2). After 12 ℎ the NES for each cluster follows the trends found using Transcriptograms.

Transcriptogram analyses for DORF6 samples
In the original paper, the authors produced microarrays for Calu3 2B4 cells infected with a mutant strain of SARS-CoV with the ORF6 gene deleted (DORF6). ORF6 binds to karyopherin α2, trapping these proteins on intracellular membranes and sequestering karyopherin β1. This sequestration may prevent nuclear import in general and hence act on transcription factors effect. The authors found that deletion of ORF6 attenuated the effects of viral infection both in vitro and in vivo, although viral titers were not significantly different from those for WT virus. To highlight the differences between the expression changes in WT and DORF6 samples, we produced time-matched relative Transcriptograms for DORF6 using WT samples as controls, as shown in Figure S 9. Overall transcription differs between WT and DORF6 after 48 ℎ, but differences between WT and DORF6 are much weaker than between WT and Mock or DORF6 and Mock. In host cells ORF6 acts to prevent transcription factors from entering the nucleus and hence promote changes in gene expression. Transcription factors regulate many different pathways, affecting gene expression in a genome-wide manner. The differences we observe between WT and DORF6, though smaller, are statistically significant and genome wide, supporting the idea that the action of transcription factors produces the observed whole-cell phenotypic changes. Complete analysis of Figure S 9 is beyond the scope of this paper. Sims and collaborators examined these differences in detail using other methods [1]. We examined expression changes for genes in the A, B and C clusters we identified using our Transcriptogram analysis for the WT virus strain. Figure S 10 reproduces relative expression averages for clusters A, B, and C for DORF6 (bright colors) and WT (pale colors) data series. See Fig. 5 in the main text. Before 24 ℎ the WT and DORF6 curves are similar for A, B and C. Between 24 ℎ and 72 ℎ, cluster B and C DORF6 curves grow steadily. The WT cluster B curve decreases after 30 ℎ and the WT cluster C curve decreases after 54 ℎ. The difference between DORF6 and WT cluster B time series is particularly dramatic, with DORF6 levels significantly and persistentlyrita higher after 5 ℎ. Since clusters B and C mainly comprise genes involved in immune response, we conclude that ORF6 affects the cellular immune response, in agreement with the original paper. The Transcriptogram analysis suggests that the most relevant differences in gene expression occur in genes in cluster B. and WT. Cluster B responds very fast to viral infection, with DORF6 differential-expression below that of WT samples at 0 ℎ, rapidly increasing to become much larger than for WT samples around 25 ℎ. Cluster C differential-expression dynamics are similar for DORF6 and WT samples, with DORF6 levels slightly lower than WT between 20 ℎ and 55 ℎ.