nf-core/rna-seq

I used the nf-core/rna-seq pipeline (v3.14.0) to analyse RNA-sequencing data.

params yaml

I used GRCh38 as the reference genome. The QC is reasonable.

The pipeline parameters were relatively vanilla. In a larger analysis, there would be a little more optimisation to do.

Transcript abundance to count matrices

I used the Bioconductor package tximport in R version 4.3.1 to summarise transcript-level estimates (salmon files) into a count matrix. Counts are length scaled transcripts per million (nTPM).

Script - 01_tximport_DESeq2.R

Outlier detection and filtering

I filtered all genes with ≤10 normalised TPM.

At this step, I also filtered some genes which are related to U2 and are massively over expressed across the board. These were:

• ENSG00000275616
• ENSG00000274062
• ENSG00000275219
• ENSG00000274862

I performed Generalised PCA on the count matrix, with the package glmpca. I identified samples which were more than 2.5 SD from the mean in either (or both) of the first 2 PCs and removed them.

• GCTU_PN0047
• GCTU_PN0052

Script - 02_outliers.R

I then re-imported the salmon files.

Script - 03_tximport_DESeq2_filtered.R

Differential expression

I used DESeq2 to perform differential expression analysis, based on subphenotype. I used the ‘il6_stnfr1_hco3_poc_alloc’ metadata column to assign samples to a subphenotype. Where this value was missing, I imputed it based on the ‘il6_stnfr1_hco3_elisa_alloc’ column. For sample GCTU_PN0090, I assumed the patient to be hypoinflammatory.

Script - 04_DE.R

At an adjusted P value threshold <0.1, for hyperinflammatory vs. hypoinflammatory, there were:

• 130 (0.47%) up-regulated genes (LFC >0)
• 438 (1.6%) down-regulated genes (LFC <0)

Applying a log2 fold change threshold of > 1.2 or < -1.2, there were:

• 37 up-regulated genes
• 104 down-regulated genes

Up-regulated genes (adj P <0.1 and abs(Log2FC) > 1.2)

Down-regulated genes (adj P <0.1 and abs(Log2FC) > 1.2)

Functional enrichment

I performed functional enrichment of all significantly differentially expressed genes (adj P <0.1 and abs(Log2FC) > 1.2), against the background of all genes with > 10 nTPM.

First, over representation analysis using GO, KEGG, Reactome, WikiPathways, and the Human Protein Atlas. I used the grprofiler2 package, with ordered queires and a FDR <0.01.

Script - 05_functional_enrichment.R

Protein-protein interaction

I created a PPI network using STRINGdb and significantly differentially regulated genes (adj P <0.1 and abs(Log2FC) > 1.2). This mapped 59 proteins. There were 22 edges, giving an enrichment P value of 3.35e-05.

Script - 06_protein_protein_interaction.R

I then MCL clustered this network, identifying 8 clusters.

The following enrich against Reactome:

• Cluster 2 - Classical antibody-mediated complement activation
• Cluster 3 - RHOB/RHOC/RAC2 GTPase cycle
• Cluster 7 - EGFR downregulation
• Cluster 8 - Interaction between L1 and Ankyrins

Overlap with existing inflammatory subphenotype RNA-seq studies

Next, I considered the overlap between DEGs in PHIND with DEGs in two other RNA-seq studies:

• Sarma, A. et al. Hyperinflammatory ARDS is characterized by interferon-stimulated gene expression, t-cell activation, and an altered metatranscriptome in tracheal aspirates. bioRxiv (2022).
• Sinha, P et al. Molecular phenotypes of acute respiratory distress syndrome in the ROSE trial have differential outcomes and gene expression patterns that differ at baseline and longitudinally over time. Am. J. Respir. Crit. Care Med. 209(7):816-828 (2024).

Script - 07_gene_overlap.R

I retrieved hyper vs. hypoinflammatory DEGs from these studies based on the same criteria - adj P value <0.1 and abs(Log2FC) > 1.2.

Sarma

Sinha day 0

Sinha day 2

Clustering signatures based on gene embedding

Given the limited overlap between studies based on gene identity. I embedded the DEGs (gene signatures) from each study using pre-derived functional embeddings:

• Chehen, H. et al. Drug target prediction through deep learning functional representation of gene signatures. Nat. Commun. 15(1):1853 (2024). https://github.com/chenhcs/FRoGS

Embeddings - embeddings.csv

To create a background, I randomly partitioned the PHIND data using the same proportions of hyper and hypoinflammatory patients, and iterated DE analysis 1000 times. This returned 208 random gene signatures with ≥ 10 genes.

Script - compute_embeddings.py

I then formally assessed whether the study derived signatures were more or less similar when compared to the background. I computed the pairwise Euclidean distance between all pairs of study signatures and from each study signature to each random signature.

I then performed paired t tests and Wilcoxon signed-rank tests - comparing the means of the within-study distances to the study-to-random distances.

• Paired t-test results: statistic = -4.52758, p-value = 0.00624
• Wilcoxon signed-rank test results: statistic = 0.0, p-value = 0.03125

The mean within-study distance was 0.28876 vs. 0.45015 from study to random signatures.

Script - embeddings_similarity.py

I also calculated an estimate of the effect size using Cohen’s d which was -2.89138, suggesting that the within study distances are smaller than study to random background and that the magnitude of this difference is substantial.

To defend against assumptions in the t test and Wilcoxon ranked-sign test, I also conducted a permutation test:

1. Observed Difference: First, I calculated the actual difference between the mean of within-study distances and the mean of study-to-random distances. This is the observed difference.
2. Combining Data: I then combined all the distances (both within-study and study-to-random) into a single pool of data.
3. Permutation Process: I randomly shuffled this combined pool of data. After shuffling, we split the data back into two groups of the same sizes as our original groups. I calculated the difference in means between these new, randomly assigned groups. I repeated this 10,000 times.
4. Comparison: For each permutation, I compared the randomly generated difference to the observed difference. I counted how many times the permuted difference is as extreme as or more extreme than our observed difference.
5. P-value Calculation: The p-value is calculated as the proportion of permutations where the permuted difference was as extreme as or more extreme than the observed difference.

The P value was < 0.0001

Script - embeddings_cohensd_and_perm.py