Irati Garmendia

Aditi Varthaman

Solenne Marmier

Mahmud Angrini

Ingrid Matchoua

Aurelie Darbois-Delahousse

Nathalie Josseaume

Pierre-Emmanuel Foy

Lubka T Roumenina

Naïra Naouar

Maxime Meylan

Sophie Siberil

Jules Russick

Pierre-Emmanuel Joubert

Karen Leroy

Diane Damotte

Audrey Mansuet-Lupo

Marie Wislez

Marco Alifano

Laurie Menger

Ignacio Garcia-Verdugo

Jean-Michel Sallenave

Olivier Lantz

Florent Petitprez

Dr Isabelle Cremer
email: isabelle.cremer@sorbonne-universite.fr

Acute influenza infection promotes lung tumor growth by reprogramming the tumor microenvironment

Keywords: non-small cell lung cancers, influenza virus, immune exhaustion, immunosuppression, drug metabolism

The authors show in lung-cancer models that acute influenza virus infection has a long-term effect on the tumor microenvironment impacting antitumor immune responses, tumor progression and survival, with implications for prophylactic and therapeutic management of lung cancer patients.

.

Introduction

The tumor microenvironment (TME) is a complex ecosystem comprising tumor cells, fibroblasts, endothelial and immune cells which strongly influences clinical outcome and response to treatments [START_REF] Gentles | The prognostic landscape of genes and infiltrating immune cells across human cancers[END_REF][START_REF] Damotte | The tumor inflammation signature (TIS) is associated with anti-PD-1 treatment benefit in the CERTIM pancancer cohort[END_REF][START_REF] Petitprez | The Tumor Microenvironment in the Response to Immune Checkpoint Blockade Therapies[END_REF]. In the past decade, the central role of the lung TME in tumor control or progression has been underlined [START_REF] Altorki | The lung microenvironment: an important regulator of tumour growth and metastasis[END_REF][START_REF] Germain | Presence of B cells in tertiary lymphoid structures is associated with a protective immunity in patients with lung cancer[END_REF]. In non-small cell lung cancer (NSCLC), the most common subtype of lung cancer, the study of the immune landscape shows robust infiltration of T lymphocytes co-expressing immune checkpoints, B lymphocytes, NK cells, macrophages and neutrophils [START_REF] Germain | Presence of B cells in tertiary lymphoid structures is associated with a protective immunity in patients with lung cancer[END_REF][START_REF] De Chaisemartin | Characterization of chemokines and adhesion molecules associated with T cell presence in tertiary lymphoid structures in human lung cancer[END_REF][START_REF] Platonova | Profound coordinated alterations of intratumoral NK cell phenotype and function in lung carcinoma[END_REF][START_REF] Lavin | Innate Immune Landscape in Early Lung Adenocarcinoma by Paired Single-Cell Analyses[END_REF](9). In addition, the TME is highly heterogeneous between patients, but how the individual TME is shaped and evolves during tumor progression is not completely understood. Several factors influence the TME, such as specific tumor mutations and tumor mutational burden (10,11), epigenetic modifications [START_REF] Hua | Genetic and epigenetic intratumor heterogeneity impacts prognosis of lung adenocarcinoma[END_REF], tumor cell metabolism [START_REF] Vanhove | The Metabolic Landscape of Lung Cancer: New Insights in a Disturbed Glucose Metabolism[END_REF], chronic inflammation, such as in the context of chronic obstructive pulmonary disease (COPD), and the lung microbiome [START_REF] Ramírez-Labrada | The Influence of Lung Microbiota on Lung Carcinogenesis, Immunity, and Immunotherapy[END_REF][START_REF] Jin | Commensal Microbiota Promote Lung Cancer Development via γδ T Cells[END_REF][START_REF] Tsay | Airway Microbiota Is Associated with Upregulation of the PI3K Pathway in Lung Cancer[END_REF][START_REF] Tsay | Lower Airway Dysbiosis Affects Lung Cancer Progression[END_REF]. However, the effect of acute respiratory viral infections has barely been investigated so far. NSCLC patients may experience influenza A virus (IAV) infection during the course of their cancer, but the impact of acute IAV on lung tumor progression and antitumor immune responses is largely unexplored.

IAV infection occurs worldwide in a seasonal cycle. It causes acute respiratory disease and results in 3 to 5 million cases of severe illness [START_REF] Iuliano | Estimates of global seasonal influenza-associated respiratory mortality: a modelling study[END_REF]. Elderly and immunocompromised individuals, such as those with lung cancer, are more likely to be infected and to develop severe respiratory disorders [START_REF] Angrini | To Vaccinate or not: Influenza Virus and Lung Cancer Progression[END_REF][START_REF] Kinslow | Influenza and mortality for non-small cell lung cancer[END_REF][START_REF] Kunisaki | Influenza in immunosuppressed populations: a review of infection frequency, morbidity, mortality, and vaccine responses[END_REF].

Among the potential mechanisms, IAV is a single-stranded RNA (ssRNA) virus that infects airways epithelial cells and is intracellularly recognized by RNA sensors, notably toll-like receptor 7 (TLR7) [START_REF] Pang | Efficient influenza A virus replication in the respiratory tract requires signals from TLR7 and RIG-I[END_REF].

We have previously shown that TLR7 is highly expressed by malignant cells in patients with NSCLC and confers a poor clinical outcome [START_REF] Chatterjee | TLR7 promotes tumor progression, chemotherapy resistance, and poor clinical outcomes in non-small cell lung cancer[END_REF]. Stimulation of tumor cells with TLR7 synthetic agonists induces NF-κB activation and Bcl-2 expression, increased tumor progression and resistance to chemotherapy in murine tumor models [START_REF] Cherfils-Vicini | Triggering of TLR7 and TLR8 expressed by human lung cancer cells induces cell survival and chemoresistance[END_REF][START_REF] Dajon | Toll like receptor 7 expressed by malignant cells promotes tumor progression and metastasis through the recruitment of myeloid derived suppressor cells[END_REF]. In this study, we sought to analyze the impact of IAV infection on the long-term progression of NSCLC using subcutaneous, orthotopic and genetically induced tumor models. We found that acute IAV infection of lung tumor-bearing mice increased tumor progression, induced tumor nodules of a higher grade and favored immunosuppression. RNA sequencing of the tumor-bearing lungs from non-infected and IAV-infected mice revealed a distinct transcriptomic profile suggesting a specific TME induced by viral infection. Finally, we show that in patients, expression of this signature in epithelial cells predicted poor overall survival in the lung adenocarcinoma (LUAD) subtype of NSCLC.

Material and Methods

Cell lines

KP cells (lung adenocarcinoma subtype expressing luciferase), derived from tumor-bearing lungs of Kras LSL-G12D/+ Trp53 flox/flox mice, were kindly provided by Dr. Laurie Menger (Institut Curie, Paris). KP cells were cultured in RPMI 1640 medium (Gibco, cat# 61870) supplemented with 10% heat-inactivated FCS (Eurobio, cat# CVFSVF00-01) and 1% penicilin/streptomycin (Gibco, cat# 15140) (the cells were not cultured more than passage 15). LLC cells (Lewis Lung Carcinoma, ATCC CRL-1642; gift from Vincenzo di Bartolo in 2015, Institut Pasteur, Paris) were stably transfected with firefly luciferase. LLC cells were cultured in DMEM + Glutamax (Gibco, cat# 61965) with 10% heat-inactivated FCS, 1% penicillin/streptomycin and 1% hygromycin B (Thermo Scientific Chemicals, cat# J60681) (the cells were not cultured more than passage 15). All cell lines were routinely tested for Mycoplasma with a kit (Ozyme, cat# LT07-710) and used only if negative. No specific authentication of the cell lines was performed.

Tumor induction in mice

Subcutaneous tumor models: B6(C)/Rj-Tyr c/c mice (Janvier Labs) were subcutaneously injected with 10 5 KP cells or 2x10 4 LLC cells in their right flanks. Tumor progression was measured every 2-3 days by bioluminescence using IVIS Lumina II In vivo Imaging System (Perkin Elmer) for KP cells and by caliper for LLC cells. To measure tumors from KP cells, mice were injected intraperitoneally with D-Luciferin (XenoLight D-Luciferin-K + Salt, PerkinElmer, cat# 122799; 150 mg Luciferin/kg body weight). Ten minutes after D-Luciferin injection, the bioluminescent signal was acquired during 5 min using the IVIS Lumina II system (Perkin Elmer, Waltham, MA, USA). During the acquisition procedure, mice were anesthetized with 2.5% isoflurane (XGI-8 Gas anesthesia system, Perkin Elmer) for 10 min. Data were analyzed with Living image software by defining a region of interest on the thoracic area of each mouse and extracting the corresponding total flux. For LLC cells, tumor growth was monitored via repeated measurements of the tumor size using a digital caliper. Tumor volume was calculated using the following formula: tumor size (mm 3 ) = (length × width × height)/8 × 4/3 × π.

Orthotopic tumor model: B6(C)/Rj-Tyr c/c mice were intravenously injected with 5x10 4 KP cells in their tail vein. Tumor progression was measured weekly by measuring bioluminescence. We choose to use B6(C)/Rj-Tyr c/c mice because, being albinos, they are more suitable for IVIS imaging.

Transgenic mice: Lung adenocarcinoma were induced in 8-12 weeks-old female Kras LSL- G12D/+ Trp53 flox/flox mice (referred to as KP mice, kindly provided by Dr. Tyler Jacks, Massachusetts Institute of Technology, Cambridge, MA) by intratracheal administration of 2x10 4 lentiviral (LTV) particles as previously described [START_REF] Dupage | Conditional mouse lung cancer models using adenoviral or lentiviral delivery of Cre recombinase[END_REF]. Lentivirus was generated from Cre recombinase-expressing plasmids in Dr. Olivier Lantz's laboratory to express the MHC class I OVA epitopes fused to the Cterminal end of the luciferase protein [START_REF] Alonso | Induction of anergic or regulatory tumor-specific CD4+ T cells in the tumor-draining lymph node[END_REF]. In vivo tumor progression was measured weekly by measuring bioluminescence using IVIS Lumina II In vivo Imaging System (Perkin Elmer).

All mice were housed in specific pathogen free (SPF) conditions in full accordance with FELASA recommendations. All procedures were performed in compliance with European Union Directive 63/2010 and received approval from the Charles Darwin Ethics Committee for animal experimentation.

(APAFIS 18589 -201901211717868).

IAV infection of tumor-bearing mice

For the subcutaneous lung tumor models, when tumors were appreciable (10-13 days post-cells injection), mice were randomly divided into the different groups: the groups of IAV-infected mice were intratumorally infected either with the Influenza A/Scotland/20/74 (H3N2 strain) virus or the A/PR/8/34 (H1N1 strain) (150 or 1000 pfu in 100 μL PBS). Non-infected control mice were injected with PBS.

For the orthotopic lung tumor model, when tumors were detectable (total flux around 1x10 5 photon/s), mice were randomly divided into two groups, one intranasally infected with H3N2 IAV (150 pfu in 40 μL PBS) and the other one with PBS.

For the genetic lung tumor model, eleven weeks after lentiviral inoculation, KP mice were infected with H3N2 IAV by intranasal inoculation (150 pfu in 40 μL PBS). Control mice received PBS alone. IAV infection was confirmed by measuring the body weight loss of the mice daily after viral infection and quantifying anti-IAV immunoglobulins. Mice were euthanized 16 weeks after lentiviral delivery.

Both IAV strains, originally obtained from Professor S. Van der Werf (Institut Pasteur, Paris) were purified as previously described [START_REF] Sanchez-Guzman | Silver nanoparticle-adjuvanted vaccine protects against lethal influenza infection through inducing BALT and IgA-mediated mucosal immunity[END_REF]. Tumor growth was followed up by IVIS Lumina II In vivo Imaging System (Perkin Elmer).

Immunohistochemistry

Mice were perfused intracardially with PBS followed by 4% paraformaldehyde (Chem Cruz, cat# sc-281692) at room temperature before dissection. Lung tissue was embedded in paraffin using a Leica ASP300/Histocore Arcadia H instrument. Histological 3-5 µm sections were processed by cutting 20 slides to perform hematoxylin (Dako, cat# 53309) and eosin (Abcam, cat# ab246824) (HE) staining.

Three slides for each were used for the quantification (using the first slide, the 10 th and the 20 th ). All slides were scanned with the NanoZoomer 2.0HT (Hamamatsu manufacturer) and analyzed by NDP.view2 viewing software (Hamamatsu manufacturer). We counted the number of tumor nodes and identified the grade of the tumors in HE sections, based on the tissue histology, as previously described [START_REF] Dupage | Conditional mouse lung cancer models using adenoviral or lentiviral delivery of Cre recombinase[END_REF] and with the help of a pathologist, who validated our quantification. Finally, the surface area of the tumor and the tissue were quantified with NDP.view2 viewing software (Hamamatsu manufacturer).

For immunofluorescence staining of murine and human tumors (see below), antigen retrieval was done at 97 °C for 30 min at pH 9 using a PT-Link instrument (Dako). Primary antibody staining was carried out for 15 min at RT, followed by the use of polymers for 10 min at RT. Antibody detection was performed with Opal technology. A list of the antibodies and other reagents used in the immunofluorescence staining can be found in Supplementary Table S1. All sections were stained with DAPI (1 μg/mL, Invitrogen, cat# D21490) for 5 min and mounted in Prolong Diamond antifade mounting (Invitrogen, cat# P36961). Slides were scanned using Axioscan (Zeiss) and the analysis and quantification were done with Halo software version 3.5 (Indica Labs).

Anti-influenza antibody quantification by ELISA

Blood collected from anesthetized mice was centrifuged at 10,000 rpm for 10 min at 4 °C. Purified IAV was coated at 1 x 6 pfu/mL overnight at 4 °C in 96 microtiter plates (Greiner). After washing, blocking was performed with 2% BSA and 0.05% Tween-20 in PBS for 2 hours at RT. Serial dilutions of plasma were prepared in 0.1% BSA, 0.05% Tween and PBS buffer, and incubated for 1.5 hours at 37°C. Mouse HRP-conjugated IgG (Southern Biotech, cat# 1030-05) or mouse HRP-conjugated IgM (Millipore, cat# AP500P) diluted at 1:10,000 were used as detection antibodies for 30 min at 37°C. Peroxidase activity was revealed using sigmaFAST OPD (Sigma, cat# P9187). The colorimetric reaction was stopped using H2SO4 and plates were read at 492 nm on the SPARK 10M reader (Tecan). IgG and IgM concentrations were determined using a standard curve of mouse purified total IgG (Sigma, cat# PP54) or IgM (Sigma, cat# PP50).

Flow cytometry of in vitro IAV-infected KP cells

In vitro IAV-treated (KP cells: 0.01 or 0.05 MOI for 24, 48 and 72 hours; LLC cells: 1 or 2 MOI for 24, 48 and 72 hours) and non-treated cells were detached using trypsin (Gibco, cat# 25300-054) and washed in RPMI 1640 10% heat-inactivated FCS (complete medium). Fc receptors were blocked using purified anti-mouse CD16/CD32 (1:50, BD Pharmingen). Cells were stained with LIVE/DEAD™ Fixable Yellow dead cell staining kit (1:100, Invitrogen, cat# L32250), anti-PD-L1 and anti-MHC class-I (H-2Kb) (information regarding the antibodies can be found in Supplementary Table S2). Cells were washed, fixed and permeabilized using the intracellular fixation and permeabilization buffer set (eBioscience, cat# 88-8823-88), according to the manufacturer's guidelines. Intracellular staining was performed using an anti-influenza A virus nucleoprotein (NP) (Supplementary Table S2). Cells were analyzed using a BD Fortessa X20 flow cytometer (BD Biosciences). Data was studied using Kaluza Analysis software, version 1.3 (Beckman Coulter).

Flow cytometry analyses of lungs from KP mice

Mouse lungs were harvested in RPMI 1640 medium (Gibco) into gentleMacs C tubes (Miltenyi) and enzymatically digested with the Miltenyi tumor dissociation kit (cat# 130-096-730) by using the gentleMACS™ octo dissociator at 37°C, according to manufacturer's instructions. Tumor adjacent normal tissue was excluded at maximum so that only tumor tissue was included in these analyses.

Tumor-bearing lung pieces were manually filtered through a 70 µm cell strainer. Erythrocytes were lysed in ACK buffer for 5 min and the reaction was halted with RPMI 1640 medium. Approximately 2 million cells were transferred into a 96-well plate for flow cytometry staining. Fc receptors were blocked with anti-mouse CD16/CD32 (1:100, BD Pharmingen) along with Zombie™ Aqua Fixable Viability dye (BioLegend) for 10 min at RT. Cells were washed with PBS followed with specific antibody staining diluted in PBS for 15 min at RT. All antibodies used are listed in Supplementary Table S2.

Mouse MHC class I tetramers

Taqman assay for viral sensors and Interferon-Stimulated Genes

RNA was extracted using the RNeasy Mini Kit (Qiagen, cat# 74136) from KP cells infected with 0.01 or 0.05 MOI of IAV for 24, 48 and 72 hours. Non-infected cells were used as control. The quality and quantity of extracted RNA were tested by using Agilent RNA 6000 Nano Kit (Agilent Technologies cat# 5067-1511) and Agilent Bioanalyzer according to the manufacturers' guidelines. Reverse transcription of RNA to cDNA was performed using 1 μg of total RNA using high-capacity cDNA Reverse transcription kit (Applied Biosystems, cat# 10400745) and Thermocycler Eppendorf nexus Mastercycler. For individual gene expression analysis, qPCR reactions were carried out by using FastStart Universal Probe Master Mix (Roche, cat# 4913949001) and inventoried TaqMan probes (FAM-MGB dye) specific for murine Tlr7 (Thermofisher Scientific, cat# Mm00446590_m1), Tlr3 (Thermofisher Scientific, cat# Mm01207404_m1), Ddx58 (Thermofisher Scientific, cat# Mm01216852_m1), Ifih1 (Thermofisher Scientific, cat# Mm00459183_m1) and Gapdh (Thermofisher Scientific, cat# Mm99999915_g1), according to manufacturer guidelines. These experiments were performed in triplicates. Gene expression was determined by fold-change (2-∆∆Ct) comparing infected and non-infected samples.

TaqMan® gene expression array was performed using the TaqMan Array Mouse Immune Panel (Applied Biosystems, cat# 4367786) following recommendations of the manufacturer, for the study of interferon-stimulated genes' (ISGs) expression. The plates were read on a 7900HT Fast Real-Time PCR System (Thermofisher).

Bulk RNA sequencing

Primary human lung tumors were obtained after resection surgery at Cochin Hospital (Paris). Written informed consent was obtained from the eight patients before participation. The patients included in this study signed an informed consent form prior to inclusion. The research was approved (CPP Ile-de-France II n°2012-06-12) by the medical ethics board at Cochin Hospital and conducted according to the recommendations in the Helsinki Declaration. Tumors were cut in two pieces of approximately 5 mm 3 : one part was infected with IAV (closed to MOI 1) for 24 hours, while the other part was not infected and served as control.

Total RNA was extracted from human tumors and tumor-bearing lungs from mice infected or not with IAV. RNA was extracted using the Maxwell 16 LEV simply RNA tissue kit (Promega, cat# AS120) and was quantified using the fluorometric qubit RNA assay (Life Technologies, cat# Q32852). RNA integrity (RIN) was determined with the Agilent 2100 bioanalyzer (Agilent Technologies). RNA sequencing was performed at the GENOM'IC platform (Cochin Institute, Paris). To construct the libraries, 1 µg of highquality total RNA (RIN ≥ 8 for murine tumors, RIN ≥ 5 for human tumors) was processed using TruSeq Stranded mRNA kit (Illumina, cat# 20020594). Briefly, after purification of poly-A containing mRNA molecules, these were fragmented and reverse-transcribed using random primers. Replacement of dTTP by dUTP during the second strand synthesis permitted high strand specificity. Addition of a single A base to the cDNA was followed by ligation of Illumina adapters. Libraries were quantified by qPCR using the KAPA library quantification kit for Illumina libraries (Kapa Biosystems, cat# 5067-4626) and library profiles were assessed with the DNA high sensitivity labchip kit on an Agilent bioanalyzer.

Libraries were sequenced on an Illumina NextSeq 500 system using 75 base-lengths read V2 chemistry in a paired-end mode. After sequencing, a primary analysis based on AOZAN software (ENS, Paris) was applied to demultiplex and control the quality of the raw data (based on FastQC modules, version 0.11.5).

Statistical and bioinformatics analyses

Comparisons of tumor growth were analyzed by using the two-way ANOVA with Sidak correction for Differential gene expression analysis from RNA-Seq data (mouse and human tumors) was performed using Galaxy Software (ARTBio platform, IBPS, Paris). Fastq files obtained were aligned using the HiSAT2 algorithm (on MM10 murine and on GRCH38 human reference genomes) and quality control of the alignment was performed with the RSeQC package. Reads were counted using htseq-count and the analysis of differential expression was performed with the DESeq2 package and Biomart (mouse GRC38p6). KEGG and GO Biological processes were interrogated and GSEA analysis was done using the R package GSEABase. These analyses were performed using R version 4.1.3.

TCGA LUAD analysis

For TCGA analysis, gene expression and clinical data from the TCGA LUAD cohort [START_REF]Comprehensive molecular profiling of lung adenocarcinoma[END_REF] were obtained from the Genomics Data Portal (https://portal.gdc.cancer.gov) for n=473 tumors. FPKM data were normalized and set to logarithmic scale. For the REVISE signature definition, human orthologs of genes overexpressed in infected compared to non-infected mice were obtained using the R package biomaRt.

The Pearson correlation matrix of the expression of these genes in the TCGA LUAD cohort was examined and clustered using hierarchical clustering with Euclidian distance and Ward criterion. Five groups were identified from the clustering, each with strong intra-group correlation and low correlation between groups. These groups were labeled 1 to 5. The impact on survival was estimated using Kaplan-Meier curves and logrank test.

In all statistical studies, *: P < 0.05, **: P < 0.01, ***: P < 0.001 and ****: P < 0.0001.

OncoSg analysis

cBioPortal (https://www.cbioportal.org) for n=169 samples. The impact on survival was estimated using Kaplan-Meier curves and logrank test with the R package survival.

GSE131907 single-cell RNA-seq analysis

For Single-cell RNA-sequencing analysis, gene expression (log2 TPM) and metadata from cohort GSE131907 [START_REF] Kim | Single-cell RNA sequencing demonstrates the molecular and cellular reprogramming of metastatic lung adenocarcinoma[END_REF] were downloaded from the Gene Expression Omnibus portal. Only cells from lung tumor samples (Sample_Origin "tLung", n = 11 samples) were kept. Data visualization and analysis was performed using the R package Seurat v4 [START_REF] Hao | Integrated analysis of multimodal single-cell data[END_REF].

COPD cohort analysis

Data was obtained from cohort GSE124180 (34) Gene counts and metadata were downloaded from the Gene Expression Omnibus portal and transformed into TPM values. Data visualization and analysis was performed using R.

Results

IAV infection enhances tumor growth in subcutaneous, orthotopic and genetically-induced autochthonous lung adenocarcinoma

In an attempt to characterize the molecular modifications induced after IAV infection, KP and LLC cells were in vitro infected (H3N2 strain). We observed productive infection as revealed by viral nucleoprotein (NP) expression, which increased between 24 and 72h after infection (Fig. 1A, Supplementary Figure 1A,B), concomitant with upregulation of the viral sensors Interferon Induced With Helicase C Domain 1 (Ifih1), TLR3, TLR7 and DExD/H-box helicase 58 (Ddx58), in addition to CD274 and MHC-I (Fig. 1A). We also found that many genes encoding molecules involved in immune response regulation including the cytokines IL-1α, TNF-α, IL-6 and IL-13, the chemokines CXCL10 and CXCL11, and Nos2 or HMOX1 were highly increased by IAV in tumor cells, suggesting that IAV may influence the TME composition (Fig. 1B). Additionally, IAV infection did not alter cell viability of LLC cells and reduced KP cell viability by 40% 72h after infection (Supplementary Fig. S1A,B).

We analyzed the effect of IAV infection in subcutaneous and orthotopic murine models of lung carcinoma cells injection. KP and LLC cells were first subcutaneously injected in B6(C)/Rj-Tyrc/c mice and intratumoral IAV infection (150 and 1000 pfu/mice) was performed. Increased tumor growth was induced by IAV infection as compared to non-infected subcutaneous tumors both in the KP (Fig. 1C-D) and in the LLC tumor model, whatever the strain of IAV and the dose of virus injected (Supplementary Fig. S2). For the orthotopic model, KP cells were inoculated in mice by intravenous injection. As soon as lung tumors were detected around 13 days post-cells injection, mice were intranasally infected with IAV (150 pfu/mice). Again, IAV infection significantly increased the growth of already established tumors (Fig. 1E-F).

We then studied the effect of intranasal IAV infection on lung tumor progression in a mouse model of genetically-induced autochthonous lung adenocarcinoma (LUAD, the most common histologic subtype of NSCLC) driven by the expression of oncogenic Kras G12D and loss of p53 (Kras LSL-G12D/+ Trp53 flox/flox , KP mice) [START_REF] Dupage | Conditional mouse lung cancer models using adenoviral or lentiviral delivery of Cre recombinase[END_REF]. Tumor progression was monitored weekly by bioluminescence imaging and when 50% of mice developed detectable tumors (total flux <2x10 5 photon/s, week 11 post-LTV inoculation), mice were intranasally administered with a sublethal dose of IAV or PBS (experimental design in Fig. 2A).

Non-LTV-inoculated mice, all of which did not develop tumors, served as controls. IAV infection in tumor-bearing mice resulted in significantly larger tumors at 5 weeks post-infection (Fig. 2B-C) and reduced survival (Fig. 2D) as compared to non-infected tumor-bearing mice. IAV infection resulted in a comparable reduced body weight in all infected mice with or without tumors (loss of 15-20%, 10 days after infection) (Fig. 2E). Furthermore, both control and tumor-bearing mice had detectable anti-IAV IgM and IgG titers in the serum following viral infection when sacrificed (Fig. 2F), which further confirmed that all mice were successfully infected. Histological characterization of the lung tumors (performed as in ( 26)) revealed higher tumor grade in IAV-infected lungs compared to the non-infected ones (Fig. 2G). Viral hemagglutinin (HA) protein was detected in both epithelial pan cytokeratin (panCK + ) tumor cells (median of 4.6% HA + cells) and CD45 + immune cells (median of 2.2% HA + cells) (Fig. 2H-I), showing that IAV infects both tumor and immune cells and persists long term in the lungs (up to 5 weeks after infection). In addition, viral infection induced high PD-L1 expression in panCK + malignant cells (from 2.8% to 32%, P=0.009) (Fig. 2J-K), consistent with the in vitro upregulation of PD-L1 after IAV infection in KP cells (Fig. 1A).

These results demonstrate that acute IAV infection of tumor and immune cells in the TME accelerates lung tumor growth.

Influenza A virus infection affects tumor-associated CD8 + T cells

To understand the underlying causes of the accelerated tumor growth upon IAV infection, we performed a comprehensive immune cell phenotyping analysis of the lungs of tumor-bearing and nontumor-bearing KP mice by flow cytometry, at sacrifice (Supplementary Fig. S3A-D describes the gating strategies). Although the distributions of CD4 + and CD8 + T cells, NK cells, B cells and myeloid cells were similar between the groups (Supplementary Fig. S4A), the T-cell compartment was markedly affected by IAV infection. In particular, IAV infection led to an increase in effector and effector memory CD8 + Tcell populations in lungs both with and without tumors, concomitant with a decrease of naïve cells (Fig. 3A, Supplementary Fig. S4B). CD8 + T cells harboring a tissue-resident memory phenotype (TRM CD8 + T) were enriched in infected lungs as compared to the non-infected ones (Fig. 3B, Supplementary Fig. S4C). In contrast, no notable changes were observed in the tumor-infiltrating effector CD4 + , TRM CD4 + and regulatory CD4 + T cell (Tregs) compartments (Supplementary Fig. S4D-F).

We then evaluated the exhaustion state of CD8 + T cells by characterizing the expression of negative immune checkpoints. Unsupervised analysis of the flow cytometry data from tumor-bearing mice with a T-cell marker panel using the Excyted pipeline (https://github.com/maximemeylan/Excyted) [START_REF] Meylan | Excyted pipeline[END_REF], which performs automated cell gating and clustering, showed sixteen T-cell clusters (Fig. 3C), based on the expression of the markers included in the panel (Fig. 3D). In-depth analysis of each cluster showed that five CD8 + T-cells clusters (clusters 2, 3, 6, 12 and 13) were differentially represented in infected and non-infected tumors (Fig. 3E). Among them, cluster 12, which was significantly overrepresented in IAV-infected tumors, corresponded to exhausted TRM CD8 + T cells expressing PD-1, TIGIT, LAG-3, CD39, CD69 and CD103 (Fig. 3E). This result was further confirmed using a classical gating strategy. The expression of PD-1, CTLA-4, TIGIT, LAG-3 or TIM-3 was significantly increased in TRM CD8 + T cells in tumor-bearing lungs following IAV infection, PD-1 being the highest expressed in IAV-infected tumor-bearing and control mice (Fig. 3F, Supplementary Fig. S4G). However, IAV did not induce expression of the five immune checkpoints in either the overall CD8 + T-cell compartment (Supplementary Fig. S4H) or the non-TRM CD8 + T cells (Supplementary Fig. S4I). antigens. IAV-specific CD8 + T cells were strongly increased in number in IAV-infected tumors (Fig. 3G).

Since malignant cells in the KP

In contrast, OVA-specific total CD8 + T cells were diminished in number in infected as compared to noninfected lung tumors (Fig. 3H). Furthermore, IAV infection also led to increased PD-1 expression in these OVA tumor-specific CD8 + T cells (from 17% to 32%, P=0.048) (Fig. 3I). Importantly, PD-1 expression was significantly higher in tumor-specific than in IAV-specific CD8 + T cells (Fig. 3J).

These data show that IAV infection promoted a decrease in the number of tumor-specific CD8 + T cells in the lungs, and that these cells are strongly exhausted, which suggests an impaired anti-tumor response in IAV-infected mice.

IAV induces minor changes in the NK-cell and B-cell compartments and no changes in myeloid cells

Unsupervised analysis of T, NK and B cells from the lungs of KP tumor-bearing mice revealed 18 clusters (Supplementary Fig. S5A), which were separated in four groups of cells (B cells, CD11b + memory B cells, T cells and NK cells), as evidenced by the expression of the markers used in the panel (Supplementary Fig. S5B). Four clusters were significantly enriched (T-cell clusters 1 and 2, B-cell cluster 3 and CD11b + memory B-cell cluster 4) and three excluded (NK-cell cluster 6, B-cell cluster 11, and CD11b + memory B-cell cluster 13) in tumors of IAV-infected mice as compared to the non-infected ones (Supplementary Fig. S5C). Since unsupervised analysis of NK cells and B cells revealed distinct profiles between the different groups of mice, we further examined B-cell subsets by classical flow cytometry evaluation. Viral infection did not modify B-cell compartments based on the expression of IgM and IgD (Supplementary Fig. S5D), nor B-cell activation (Supplementary Fig. S5E). In addition, no significant impact of viral infection could be observed on myeloid-cell subsets (Supplementary Fig. S5F). Altogether, these results show modifications in B-cell and NK-cell clusters.

IAV infection induces a distinct transcriptional profile in tumor-bearing lungs

To decipher the mechanisms involved in the accelerated lung tumor growth we observed upon IAV infection, we analyzed the transcriptional profile of KP tumor-bearing and non-tumor-bearing lungs from IAV-infected and non-infected mice, by RNA-Seq analysis at sacrifice point. The four groups of mice were delineated in four well-defined clusters in a principal component analysis (PCA) (Fig. 4A).

IAV infection led to the upregulation of 273 genes and downregulation of 12 genes, as represented in the volcano plot (Fig. 4B).

Enrichment analyses on the 273 upregulated genes were carried out using Gene Ontology (GO) and KEGG pathways. They revealed enrichment of pathways related to immune response (Fig. 4C), antiviral response, carcinogenesis, lipid metabolism and drug metabolism (Fig. 4D). To further annotate the biological functions enriched in the 273 upregulated genes, we performed gene set enrichment analyses (GSEA). Here again, we observed upregulation of pathways involved in response to viruses (Fig. 4E), chemical carcinogenesis (Fig. 4F), drug metabolism (Fig. 4G-J) and cAMP signaling pathway (Fig. 4K). We performed an additional differential gene expression analysis between the four groups of mice by pairwise comparison. The Venn-diagram showed 52 genes (Supplementary Table S3) that were specifically overexpressed only in IAV-infected tumors as compared to the other groups (Fig. 4L).

Altogether, our data show that IAV infection induces a distinct transcriptional signature that distinguishes infected from non-infected tumors. The overexpressed genes are associated with antiviral response, lipid metabolism and drug metabolism, as well as an immune response-related profile.

IAV infection induces a distinct transcriptional profile in ex vivo cultures of primary human lung tumors

To characterize the effect of IAV infection in human lung tumors, we infected ex-vivo eight primary tumor pieces (7 adenocarcinoma, 1 squamous cell carcinoma) obtained immediately after surgery at the Cochin Hospital (Paris, France) with IAV (MOI around 1). Non-infected tumor pieces served as control. We detected viral infection of tumor cells by HA staining of cytokeratin expressing tumor cells (Supplementary Fig. S6A). We then examined the transcriptional profile by RNA-Seq analysis 24h after infection. We found 27 overexpressed and 17 downregulated genes (log2 Fold Change>1 and P<0.05) (Supplementary Fig. S6B, Supplementary Table S4) associated with immune response processes (Supplementary Fig. S6C) and the immune antiviral response (Supplementary Fig. S6D). cGMP-PKG, PPAR and Ras signaling pathways, influenza A pathway and TLR signaling pathways were also enriched in KEGG analysis (Supplementary Fig. S6E). These results demonstrate that IAV infection modulates the human lung TME, affecting both tumor-associated features (e.g., lipid metabolism and oncogenesis) and the immune response, at early stages of the infection.

The REVISE signature identifies patients with LUAD with poor survival

Because IAV infection induced a specific gene expression signature in lung tumor-bearing mice, we analyzed the cohort of human LUAD patients from TCGA [START_REF]Comprehensive molecular profiling of lung adenocarcinoma[END_REF] in order to identify such a signature and study its impact on patients' clinical outcome. Of the 273 genes upregulated in the mouse model, we identified 217 human orthologs with expression values in the TCGA dataset. A correlation heatmap of these genes revealed 5 modules, labeled 1 to 5 (Fig. 5A, Table 1). These subsignatures were identified by a strong internal correlation and by their low correlation between each other. To estimate which cell populations expressed these modules, we analyzed a single-cell RNA-seq dataset of lung adenocarcinoma from GSE131907 [START_REF] Kim | Single-cell RNA sequencing demonstrates the molecular and cellular reprogramming of metastatic lung adenocarcinoma[END_REF]. We selected only lung tumor samples. Our cell type annotation based on the data provided by the original authors is shown in Fig. 5B. We observed the expression pattern of the full signature of 217 genes and of each module (Fig. 5C-D). Only module 1 was specific of epithelial cells. We named this gene signature module REVISE (Respiratory Virus Infection Signature in Epithelium, Table 1).

We further investigated the clinical relevance of the REVISE signature in human lung adenocarcinoma.

Given the distribution of the REVISE expression in the TCGA LUAD cohort, we chose the third quartile as cut-off between high and low expression (Fig. 5E). A high expression of the REVISE signature was found to be associated with a poor prognosis in the TCGA LUAD cohort (Fig. 5F). To validate this, we obtained data from an independent cohort of human lung adenocarcinoma, OncoSg [START_REF] Chen | Genomic landscape of lung adenocarcinoma in East Asians[END_REF]. We similarly separated the cohort into high and low expression of the REVISE signature (Fig. 5G) and found that there was a trend towards a significantly worse prognosis for patients with tumors having a high REVISE signature expression (Fig. 5H). Finally, to estimate whether the expression of the REVISE signature stems from previous infections of the lung, we analyzed RNA-seq data from large airways epithelium in patients with or without chronic obstructive pulmonary disease (COPD) (Supplementary Fig. S7). No significantly different expression of the REVISE signature was observed between patients experiencing COPD and the control cases, which could indicate that the expression of the REVISE signature in the tumor is related to tumor-specific mechanisms and is not due to pre-existing inflammation like COPD.

Discussion

Patients with lung cancer frequently exhibit a compromised immune response that makes them prone to microbial infection. Reciprocally, the consequences of viral infections and the ensuing response on lung cancer progression remains ambiguous. Several studies and clinical case reports have suggested that vaccination against IAV, yellow fever virus and even, SARS-CoV-2, leads to tumor regression [START_REF] Aznar | Repurposing the yellow fever vaccine for intratumoral immunotherapy[END_REF][START_REF] Sousa Lg De | Spontaneous tumor regression following COVID-19 vaccination[END_REF][START_REF] Newman | Intratumoral injection of the seasonal flu shot converts immunologically cold tumors to hot and serves as an immunotherapy for cancer[END_REF]. It has also been shown in murine models, that IAV infection prior to lung tumor initiation enables better control of tumor progression because of the expression of abnormal self-antigens before tumor challenge, which induce immune memory and cancer immunosurveillance [START_REF] Masemann | Oncolytic influenza virus infection restores immunocompetence of lung tumor-associated alveolar macrophages[END_REF][START_REF] Iheagwara | Influenza virus infection elicits protective antibodies and T cells specific for host cell antigens also expressed as tumor-associated antigens: a new view of cancer immunosurveillance[END_REF]. On the other hand, recent studies have shown that SARS-CoV-2 infection leads to a worse outcome in cancer patients [START_REF] Winter | Blood Cytokine Analysis Suggests That SARS-CoV-2 Infection Results in a Sustained Tumour Promoting Environment in Cancer Patients[END_REF], and epidemiological studies evidenced an increased risk of developing lung cancer due to cumulative exposure to IAV [START_REF] Inaida | Previous Infection Positively Correlates to the Tumor Incidence Rate of Patients with Cancer[END_REF][START_REF] Weng | Association between the risk of lung cancer and influenza: A population-based nested case-control study[END_REF]. In line with our observations herein that acute IAV infection leads to tumor progression, Kohlhapp and colleagues demonstrated that lung IAV infection of mice harboring distal melanoma enhances tumor growth and reduces survival [START_REF] Kohlhapp | Non-oncogenic acute viral infections disrupt anti-cancer responses and lead to accelerated cancer-specific host death[END_REF]. We found that tumor cells are infected with IAV and respond with specific transcriptional programs. Indeed, responses to live virus infection and vaccination can be distinct with vaccinations bringing additional adjuvant-induced responses. T cells specific to viral epitopes have also been found in tumors and, upon specific reactivation can lead to tumor growth control by creating an inflammatory environment within the tumor [START_REF] Rosato | Virus-specific memory T cells populate tumors and can be repurposed for tumor immunotherapy[END_REF]. We found that, upon acute IAV infection, the expansion of virus-specific T cells was accompanied by an exhaustion of tumor antigen-specific T cells. Together, the response of tumor cells to virus infection and the expansion of virus-specific T cells to acute IAV infection led to a lower anti-tumor immune response and an increase in tumor growth.

In conclusion, differences observed between studies could be argued due to the different mouse tumor models and viral concentrations used in each of them. In this respect, the KP mouse model used in the present study, harboring Kras G12D mutation with p53 deletion, finely resembles human lung cancer, and IAV infection may mimic its impact on human patients. Additionally, this model is suitable to study the impact of IAV infection, since the Kras mutation itself is unlikely to directly affect the immune cells in the TME. Indeed, in a previous study, we found no difference in immune cell densities (mature DC and CD8 + T cells) between mutated and wild-type KRAS (10).

Mechanistically, IAV HA protein persisted in malignant cells and in immune cells five weeks postinfection suggesting sustained responses in the TME. IAV considerably impaired tumor-specific T-cell responses. TRM CD8 + T cells were notably increased in the lung tumors following IAV infection, but, crucially, they were exhausted and unable to support an antitumor immune response. IAV infection induced the expression of checkpoint molecules, including PD-1, on tumor-specific TRMs, as well as a strong expression of PD-L1 in the lung epithelial tumor cells. PD-1 expression was significantly higher in tumor-specific TRMs as compared to IAV-specific TRM CD8 + T cells, suggesting that the antigen specificity of the TRMs drove their exhaustion. Echoing our results, Kohlhapp and colleagues showed that upon IAV lung infection, distal melanoma-specific CD8 + T cells migrated to the lungs and expressed high levels of PD-1. Although the main target in their study were not TRM but circulating T cells, PD-1 therapeutic blockade inhibited tumor growth in IAV-infected mice suggesting that PD-1 expression might not be related to TRM cells alone [START_REF] Kohlhapp | Non-oncogenic acute viral infections disrupt anti-cancer responses and lead to accelerated cancer-specific host death[END_REF]. In addition, although to a lesser extent as compared to PD-1, the expression of CTLA-4, TIGIT, LAG-3 and TIM-3 was also stimulated by IAV in our study and their involvement may be interesting to study further in our mouse model. The effect of IAV infection on T-cell exhaustion or T-cell dysfunction has already been described in humans [START_REF] Liu | Cycling and activated CD8+ T lymphocytes and their association with disease severity in influenza patients[END_REF] and mice [START_REF] Erickson | Viral acute lower respiratory infections impair CD8+ T cells through PD-1[END_REF] in responses to different viral strains [START_REF] Ou | Enhanced Programmed Cell Death Protein 1/Programmed Cell Death Ligand 1 Expression Induced by Severe Influenza A Virus Infection Impairs Host's Antiviral Response[END_REF] and may be due to the infection of antigen-presenting cells, which impairs their ability to cross-present antigens to CD8 + T cells [START_REF] Smed-Sörensen | Influenza A virus infection of human primary dendritic cells impairs their ability to cross-present antigen to CD8 T cells[END_REF].

Regarding the NK-cell and B-cell compartments, our unsupervised analysis revealed significant differences in seven clusters in the tumors of IAV-infected and non-infected mice. Among them, B-cell and CD11b + memory B-cell clusters were indistinctly over-and under-represented. The markers included in the panel were not sufficient to more deeply analyze the specific characteristics of each of these clusters. The current study has implications for the evaluation of the B-cell compartment in future studies.

Using RNA sequencing, we identified genes overexpressed in IAV-infected tumor-bearing mice. Strikingly, some of these genes were also up-regulated when primary human lung tumors were infected ex vivo with IAV. Despite the different kinetics of the experiments (RNA-Seq experiments were performed 5 weeks and 24h post-infection for in vivo mouse models and ex vivo primary human tumors, respectively), the induction of IFIT1, OASL and DUOXA2 was found in both cases, highlighting their importance in the process from an early time point. These three genes are well-known for their antiviral role against respiratory viruses [START_REF] Zhao | A Long Non-coding RNA IVRPIE Promotes Host Antiviral Immune Responses Through Regulating Interferon β1 and ISG Expression[END_REF][START_REF] Zhu | Antiviral activity of human OASL protein is mediated by enhancing signaling of the RIG-I RNA sensor[END_REF][START_REF] Grasberger | Identification of the maturation factor for dual oxidase. Evolution of an eukaryotic operon equivalent[END_REF], although their implication in lung cancer has not been deeply studied yet. It will be interesting to assess the implication of these genes in lung cancer.

The lack of data on IAV infection experienced by patients from the public cohorts LUAD TCGA and OncoSg cohort prevented us from confirming its direct implication in the accelerated lung tumor progression. However, based on the murine signature, we could translate our data from murine experiments to human. The analysis of differentially expressed genes in the cohort of 273 lung adenocarcinoma patients from the TCGA database revealed a novel signature (REVISE), expressed by epithelial cells, that is likely predictive of an inflammatory gene signature that supports tumor growth, and that could stem from IAV or other previous infections, which may also affect the TME and tumor progression. Of clinical interest, the high expression of the REVISE signature was significantly associated with a poorer overall survival. We propose that the REVISE signature may be used to identify patients with poor clinical outcome. In summary, our findings indicate a pro-tumoral role of IAV infection in established lung tumors and highlight the impact of the virus on the TME and may explain the epidemiological data showing increased cancer-specific deaths following IAV infection [START_REF] Li | Influenza in hospitalised patients with malignancy: a propensity score matching analysis[END_REF].

  (H-2K(b) chicken OVA 257-264 SIINFEKL-PE tetramer and H-2D(b) Influenza A PA (polymerase acidic protein) 224-233 SSLENFRAYV-BV421 tetramer) were kindly provided by the NIH Tetramer Core Facility. Tetramer staining was performed in PBS at RT for 15 min prior to staining of cell surface proteins. The anti-CD8 clone KT15 (at 1:200 dilution) was used in concert with tetramer staining to avoid antibody competition with the CD8 molecule for MHC class I molecule binding. After staining, cells were resuspended in PBS 4% paraformaldehyde for 10 min, washed twice in PBS and analyzed using a BD Fortessa X20 flow cytometer (BD Biosciences). Absolute cell counts were determined by using CountBright Absolute Counting Beads (Invitrogen, cat# C36950) following the manufacturer's recommendations. Results were plotted and analyzed using the FlowJo 10 software. Unsupervised analysis of the flow cytometry data was performed using RStudio software 1.3.959 and the Excyted pipeline (https://github.com/maximemeylan/Excyted) (29). T cells and NK and B cells panels were analyzed by the unsupervised strategy. Intensity values of events gated from live cells were normalized with the Logicle transformation. Uniform Manifold Approximation and Projection (Umap) calculation and unsupervised clustering was performed with 10,000 events for each sample using k=30 and Rphenograph (Leland McInnes arXiv:1802.03426, https://github.com/icyto/Rphenograph). Further re-clustering in the T cells panel was carried out by re-running the Excyted pipeline on selected CD3-positive clusters using k=100.

  multiple comparisons. Correlation between IAV infection and the expression of viral sensors, PD-L1 and MHC-I was calculated by Spearman correlation. Differences in cytokine expression between IAVinfected and non-infected KP cells were studied by linear regression with Benjamini-Hochberg correction for multiple testing. The percentage of tumor nodules with different grades per total lung lobule were compared between groups with the Chi-Squared test. Data obtained from immunofluorescent staining were compared with the Mann-Whitney U test. Differences in IgM and IgG anti-IAV were studied with the ordinary one-way ANOVA with Sidak correction for multiple comparisons. Data obtained from flow cytometry and studied by the supervised gating strategy were analyzed by multiple Mann-Whitney tests with Benjamini-Hochberg correction, except for data concerning Tregs, B cells activation and myeloid subsets (1-way ANOVA) as well as PD-1 expression in OVA-tetramer + as compared to IAV-tetramer + CD8 + T cells (Mann-Whitney U test). Differences in clusters found by the unsupervised analysis of flow cytometry data were studied by the Mann-Whitney U test. All these statistical analyses were performed using GraphPad Prism software (version 9.2.0).

  mice express the SIINFEKL OVA peptide (27), we used specific MHC class I tetramers to analyze the specific T-cell responses against tumor (OVA tetramer) or viral (IAV tetramer)

FiguresFig. 1 .

 1 Figures

Fig. 2 .Fig. 4 .Fig. 5 .

 245 Fig. 2. Intranasal infection by IAV in tumor-bearing lungs significantly enhances tumor growth with

Table 1 . Genes included in REVISE signature defined in the cohort of LUAD patients from TCGA.

 1

	REVISE signature
	MSLN
	ERN2
	FXYD3
	PLA2G10
	AGR2
	FAM83E
	ENTPD8
	SPINK4
	KCNE3
	MISP
	SDCBP2
	TFF1
	GCNT3
	FA2H
	BCAS1
	CDHR5
	TM4SF5
	ANKS4B
	HNF4A
	EPS8L3
	ONECUT2
	KRT20
	HNF4G
	C2orf72
	VIL1
	CDH17
	ANXA13
	BAIAP2L2
	SMIM24
	HSD17B2
	LGALS4
	CYP3A5
	TM4SF4
	DUOX2
	DUOXA2
	SMPD3
	IHH
	ALDOB
	NR1I2
	ANXA10
	MYO1A
	FER1L6
	UGT2A3
	TFF2
	REG4
	PSCA
	ARL14
	SULT1B1
	PTPRR
	CYP3A4
	REP15

Acknowledgements

We thank the Cochin hospital for contributing to the tissue collection. We also thank the CHIC (Center of Histology Imaging and Cytometry) and the CEF (Center for animal experiment) facilities of the Centre de Recherche des Cordeliers. This work was supported by the "Institut National de la Sante et de la Recherche Medicale" (INSERM), Sorbonne Universite, Universite de Paris, the LabEx Immuno-Oncology, the Institut National du Cancer (2016-PLBIO) and foundation ARC.

This work was supported by the "Institut National de la Sante et de la Recherche Medicale" (INSERM), Sorbonne Universite, Universite Paris Cite, the LabEx Immuno-Oncology, the Institut National du Cancer (2016-PLBIO) and foundation ARC.

Data availability statement

The data generated in this study are available at the European Nucleotide Archive (ENA) study, accession number: PRJEB59490 or from the corresponding author upon reasonable request.


Asmaa Rassil
email: asmaarassil@gmail.com

Hiba Chougrad
email: chougrad.hiba@gmail.com

Hamid Zouaki
email: hamidzouaki@yahoo.fr

Graph Neural Networks to Advance Anticancer Drug Design

Keywords: Graph Neural Networks, Transductive learning, Inductive learning, Representation vectors, Anticancer activity prediction

Predicting the activity of chemical compounds against cancer is a crucial task. Active chemical compounds against cancer help pharmaceutical drugs producers in the conception of anticancer medicines. Still the innate way of representing chemical compounds is by graphs, the machine learning algorithms can not handle directly the anticancer activity prediction problems. Dealing with data defined on a non-Euclidean domain gave rise to a new field of research on graphs. There has been many proposals over the years, that tried to tackle the problem of representation learning on graphs. In this work, we investigate the representation power of Node2vec for embedding learning over graphs, by comparing it to the theoretical framework Graph Isomorphism Network (GIN). We prove that GIN is a deep generalization of Node2vec. We then exert the two models Node2vec and GIN to extract regular representations from chemical compounds and make predictions about their activity against lung and ovarian cancer.

Introduction

Application domains of Artificial Intelligence (AI) became concerned in the last few years with treating complex problems defined on irregular data structure as graphs. The conception of pharmaceutical drugs [START_REF] Li | Learning deep generative models of graphs[END_REF][START_REF] You | Graph convolutional policy network for goal-directed molecular graph generation[END_REF][START_REF] De Cao | Molgan: An implicit generative model for small molecular graphs[END_REF] through the prediction of chemical compounds activities against cancer [START_REF] Wale | Comparison of descriptor spaces for chemical compound retrieval and classification[END_REF][START_REF] Shervashidze | Weisfeiler-lehman graph kernels[END_REF] need more than tabular representations to be solved. The natural way of modeling molecular data for an anticancer activity prediction problem is by graphs. Representing chemical compounds in the shape of graphs keep the structure of the data unchanged, and so no information from the initial collected data will get lost. Graphs are irregular types of data which are defined on non-Euclidean domains, therefore, all the classic machine learning algorithms cannot be directly applied on this type of structure [START_REF] Bronstein | Geometric deep learning: going beyond euclidean data[END_REF]. To deal with such complex data, a number of works on representation learning on graphs have been proposed and proved to be effective in resolving these tasks. These representation approaches aim to map nodes or the entirely given graph to an embedding space of low-dimension (figure 1). The embedding vectors are supposed to summarize the information on top of nodes and their surrounding neighbors [START_REF] Hamilton | Inductive representation learning on large graphs[END_REF]. To do representation learning on graphs two types of methods are proposed transductive and inductive learning. Inductive approaches [START_REF] Atwood | Diffusion-convolutional neural networks[END_REF][START_REF] Hamilton | Inductive representation learning on large graphs[END_REF][START_REF] Gilmer | Neural message passing for quantum chemistry[END_REF][START_REF] Chen | Fastgcn: fast learning with graph convolutional networks via importance sampling[END_REF][START_REF] Ying | Hierarchical graph representation learning with differentiable pooling[END_REF] for representation learning on graphs do not learn the embedding vectors directly from the graph but use instead the training graphs to learn a parametric encoding function to generate embeddings for unseen graphs. Transductive approaches [START_REF] Belkin | Laplacian eigenmaps and spectral techniques for embedding and clustering[END_REF][START_REF] Ahmed | Distributed large-scale natural graph factorization[END_REF][START_REF] Cao | Grarep: Learning graph representations with global structural information[END_REF][START_REF] Ou | Asymmetric transitivity preserving graph embedding[END_REF][START_REF] Shen | Discrete network embedding[END_REF][START_REF] Yang | Binarized attributed network embedding[END_REF] associate input graphs to embedding matrices of N rows and d columns. The parameter N refers to the number of nodes in the given graph and d refers to the dimensionality of the embedding space. Therefore, by transforming the irregular structure of graphs into Nxd matrices, we can straightforwardly apply Machine Learning algorithms to the data. In this work, we investigate the representation power of Node2vec for embedding learning over graphs, by comparing it to the inductive theoretical framework Graph Isomorphism Network (GIN) that proved to be equivalent to the Weisfeiler-Lehman (WL) test of isomorphism. We demonstrate theoretically that GIN is a deep generalization of Node2vec. We employ Node2vec and GIN models to define regular representations of chemical compounds represented by graphs and use these regular representations to predict the activities of the corresponding chemical compounds against lung and ovarian cancer in NCI1 and NCI109 [START_REF] Wale | Comparison of descriptor spaces for chemical compound retrieval and classification[END_REF][START_REF] Shervashidze | Weisfeiler-lehman graph kernels[END_REF] datasets respectively.

Problem Formulation

A graph G = (V, E) is defined by the vertex set V = {v i : i = 1, 2, ..., N } and the edges e ij ∈ E connecting pairs of nodes v i and v j in the graph. For each node v in the graph , we note its neighbors set by N (v) = {u : (u, v) ∈ E}. The adjacency matrix A = (a ij ) 1≤i≤N,1≤j≤N outlines the graph structure by setting a ij = 1 for connected nodes e ij ∈ E and a ij = 0 if e ij / ∈ E.

Graph Embedding

Graph embedding frameworks [START_REF] Belkin | Laplacian eigenmaps and spectral techniques for embedding and clustering[END_REF][START_REF] Ahmed | Distributed large-scale natural graph factorization[END_REF][START_REF] Cao | Grarep: Learning graph representations with global structural information[END_REF][START_REF] Ou | Asymmetric transitivity preserving graph embedding[END_REF][START_REF] Shen | Discrete network embedding[END_REF][START_REF] Yang | Binarized attributed network embedding[END_REF] are unsupervised approaches for learning on graph structured data. The main purpose of network embedding is to map input graphs to a regular low dimensional space where the use of standard machine learning algorithms become feasible [START_REF] Hamilton | Representation learning on graphs: Methods and applications[END_REF][START_REF] Cui | A survey on network embedding[END_REF][START_REF] Cai | A comprehensive survey of graph embedding: Problems, techniques, and applications[END_REF]. To preserve the graph topology and the structural information on graphs in the embedding space, the network embedding variants propose pairwise similarity measures upon the graph structure S G (u, v) and learn the node embeddings by preserving the defined similarity measure in the embedding space S E (e u , e v ) through the following loss function:

L[S G (u, v), S E (e u , e v )]
where e v ∈ R d is the corresponding embedding vector of the node v ∈ V in the embedding space E = R d and L(.) is a loss function. Two principal branches of graph embedding are proposed in the literature [START_REF] Hamilton | Representation learning on graphs: Methods and applications[END_REF] matrix factorization [START_REF] Belkin | Laplacian eigenmaps and spectral techniques for embedding and clustering[END_REF][START_REF] Ahmed | Distributed large-scale natural graph factorization[END_REF][START_REF] Cao | Grarep: Learning graph representations with global structural information[END_REF][START_REF] Ou | Asymmetric transitivity preserving graph embedding[END_REF] and random walk [START_REF] Perozzi | Deepwalk: Online learning of social representations[END_REF][START_REF] Grover | node2vec: Scalable feature learning for networks[END_REF]. Matrix factorization approaches define the proximity measure over the graph S G (u, v) either by an adjacency-based pairwise proximity measure (S G (u, v) = A uv ) [START_REF] Ahmed | Distributed large-scale natural graph factorization[END_REF], or by multihop pairwise proximity measure captured by the k-th power of the adjacency matrix of the graph (S G (u, v) = A k uv ) [START_REF] Cao | Grarep: Learning graph representations with global structural information[END_REF] (see figure 2). The adjacency-based approaches push adjacent nodes (u, v) ∈ E in the graph S G (u, v) = A uv = 1 to be equal in the embedding space S E (u, v) = e T v e u = 1 through the following loss function:

A uv -e T v e u 2 2 (1)
Multi-hop pairwise proximity measures (S(u, v) = A k uv ) [START_REF] Cao | Grarep: Learning graph representations with global structural information[END_REF] push adjacent nodes captured by A k to be equal in the embedding space S E (u, v) = e T v e u through the following loss function:

A k uv -e T v e u 2 2 (2)
Thus, nodes located in k-hop distance from each other, will have similar embedding vectors. Random walk [START_REF] Perozzi | Deepwalk: Online learning of social representations[END_REF][START_REF] Grover | node2vec: Scalable feature learning for networks[END_REF] approaches for graph embedding are heuristic methods for embedding learning. To learn embeddings, authors in [START_REF] Perozzi | Deepwalk: Online learning of social representations[END_REF][START_REF] Grover | node2vec: Scalable feature learning for networks[END_REF] optimize the probability of visiting nodes v that tend to co-occur in random walks on the graph starting from a node u Matrix factorization gives a deterministic definition of similarity upon the graph. Otherwise, random walk methods adopt a stochastic proximity measure based on randomly generated walks on the graph. Through two searching parameters p and q, Node2vec [START_REF] Grover | node2vec: Scalable feature learning for networks[END_REF] controls the propagation of the generated walk. For p > max(q, 1) we reduce the probability of sampling an already visited node in the next two steps, therefore distant nodes from a source node u are pushed to have similar embeddings [START_REF] Grover | node2vec: Scalable feature learning for networks[END_REF]. However, for p < min(q, 1), the walk stays local, then closest nodes to the source are pushed to have similar embeddings (figure 3). Fig. 3. The parameters p and q control the propagation of the generated walk on the graph. For small values of p, local neighbors are likely to be visited. However, when setting the parameter q to a small value, distant nodes tend to co-occur in the random walk

p(v | u) = exp(e T v e u ) vi∈V exp(e T v e vi ) (3)

Graph Neural Networks

Contrary to graph embedding, graph neural networks (GNNs) [START_REF] Atwood | Diffusion-convolutional neural networks[END_REF][START_REF] Hamilton | Inductive representation learning on large graphs[END_REF][START_REF] Gilmer | Neural message passing for quantum chemistry[END_REF][START_REF] Chen | Fastgcn: fast learning with graph convolutional networks via importance sampling[END_REF][START_REF] Ying | Hierarchical graph representation learning with differentiable pooling[END_REF] are deep and inductive approaches for representation learning on graphs. Through an end-to-end network, GNNs learn jointly the embeddings or representation vectors of the nodes and solve the defined problem on the graph structure. Thus, the representation vectors are learned based on the intrinsic topology of the graph and the task we aim to solve. The representation learning process involves two steps (i)aggregating neighboring nodes then (ii)combining them with the node being processed through the following recursive schema:

a (k) v = AGGREGAT E (k) ({h (k-1) u : u ∈ N (v)}) (4)
h (k) v = COM BIN E (k) (h (k-1) v , a (k) v ) (5)
where h

(k) v
is the learned representation vector of node v in the kth layer of the network. Variants of aggregation and combination functions give the GNN variants. Graph Isomorphism Network ( GIN ) [START_REF] Xu | How powerful are graph neural networks? International Conference on Learning Representations[END_REF] is a theoretical variant of GNNs which combines the aggregation and combination step in the following formula:

h (k) v = M LP (k) [(1 + (k) )h (k-1) v + u∈N (v) h (k-1) u ] (6)
where (k) is a model parameter. Graph Isomorphism Network (GIN) [START_REF] Xu | How powerful are graph neural networks? International Conference on Learning Representations[END_REF] is as powerful as the Weisfeiler-Lehman (WL) test of isomorphism [START_REF] Weisfeiler | A reduction of a graph to a canonical form and an algebra arising during this reduction[END_REF]. Thus, the learned graph representation vectors are as powerful as WL in capturing differences of graph topologies.

Deep Generalization of Node2vec

Learning embedding vectors by maximizing the likelihood of visiting the neighbors generated through a combination of Breadth-first Sampling (BFS) and Depth-first Sampling (DFS) strategies bestow Node2vec [START_REF] Grover | node2vec: Scalable feature learning for networks[END_REF] with an impressive representation power. When setting the parameter p to a small value (figure 3), nodes belonging to similar communities [START_REF] Yang | Overlapping communities explain core-periphery organization of networks[END_REF] in a graph are forced to learn approximately the same embeddings. However, when setting the parameter q to a small value (figure 3), nodes representing the same structural role in different communities [START_REF] Henderson | Rolx: structural role extraction & mining in large graphs[END_REF] are pushed to learn similar embeddings. Therefore, according to the given task Node2vec [START_REF] Grover | node2vec: Scalable feature learning for networks[END_REF] is able to learn suitable embeddings. In this work, we demonstrate the analogy between Node2vec [START_REF] Grover | node2vec: Scalable feature learning for networks[END_REF] and the powerful theoretical framework GIN [START_REF] Xu | How powerful are graph neural networks? International Conference on Learning Representations[END_REF]. And, we prove that Node2vec is closely equivalent to the Weisfeiler-Lehman(WL) test of isomorphism [START_REF] Weisfeiler | A reduction of a graph to a canonical form and an algebra arising during this reduction[END_REF].

We consider a shallow neural network with one hidden layer h as described in figure 4. Given a node u, we maximize the probability of visiting a node v ∈ N (v) by generating a random walk starting from u. We initialize the embedding of u by its corresponding column a u ∈ R N in the adjacency matrix  = A + I N , then we feed it as input to the network (figure 4) . We denote the N × d matrices corresponding to the model parameters in the input-to-hidden and the hiddento-output sections of the network by W and W respectively. We assume that W and W are embedding matrices of the input graph. Each row of W and W is a d-dimensional embedding vector e v of a node v ∈ V. In the forward pass,

W T a u = v∈N (u)∪u W T .,k = e u + v∈N (u) e T v = h u = h
the updated embedding vector h u ∈ R d of the input node u is a summation of the embedding vectors e u of u and its surrounding neighbors {e v : v ∈ N (u)} given by the embedding matrix W . Node2vec adopts a sum-aggregator to aggregate neighboring nodes and utilizes the collected information to update the initial node embedding a u . Thus, we may proclaim that node2vec is a shallow approximation of the deep and inductive Graph Isomorphism Network GIN (Eq.6).

The output of the network before applying activation is,

z i = W T .,i h = W T .,i h u = h T vi h u ∀i = 1, ..., N
using a log-softmax activation function Φ(.), we obtain (Eq.3) corresponding to Node2vec:

Φ(h T vi h u ) = log p(v i | u) = log exp(h T vi h u ) v∈V exp(h T v h v ) (7)
We can conclude that, Node2vec learns embedding vectors of the graph nodes using approximately the same process of one layer GIN. Now, considering that GIN can capture and distinguish between different graph topologies as would the WL test of isomorphism, we may assert that embeddings learned by Node2vec will also be powerful in capturing differences of graph topologies.

Experiments and Results

To predict anticancer activity of chemical compounds against lung and ovarian cancer, we investigate the representation power of Node2vec for embedding learning on graphs and two variants of the inductive Graph Isomorphism Network (GIN) using two benchmark datasets NCI-1 and NCI-109 for lung and ovarian cancer respectively. In the following paragraphs we provide a description of the datasets, details about the models architecture and configurations, and the obtained results achieved by each model for predicting anticancer activity.

Datasets

The NCI-1 and NCI-109 graph datasets [START_REF] Wale | Comparison of descriptor spaces for chemical compound retrieval and classification[END_REF][START_REF] Shervashidze | Weisfeiler-lehman graph kernels[END_REF] are provided by the National Cancer Institute. Each NCI dataset contains thousands of chemical compounds modelled as graphs. The atoms are represented as nodes and bonds as edges in the graph representation of the chemical compounds (figure 5). Statistics about NCI-1 and NCI-109 graph datasets are reported in Table 1. NCI-1 dataset contains two categories of chemical compounds active and nonactive against the non-small cell lung cancer. Active compounds are labeled by +1 and non-active by -1. NCI-109 dataset contains two categories of chemical compounds active and non-active against ovarian cancer. Active compounds are labeled by +1 and non-active by -1.

Implementation Details and Results

To predict the anticancer activity of chemical compounds against lung and ovarian cancer, we train the embedding model Node2vec [START_REF] Grover | node2vec: Scalable feature learning for networks[END_REF] on NCI-1 and NCI-109 datasets with the following hyper-parameters: embedding dimensionality ∈ {7, 12}, random walk length ∈ {3, 4}, size of neighbors set (i.e the context size) ∈ {5, 7} and for each node we generate 10 random walks.

To train the deep Graph Isomorphism Network (GIN) [START_REF] Xu | How powerful are graph neural networks? International Conference on Learning Representations[END_REF] on NCI-1 and NCI-109 datasets, we use two variants of the model using (1) multilayer perceptron and (2) single-layer perceptron as combination functions of neighboring node representations. To build the anticancer activity prediction model, we use five GIN layers each followed by batch normalization [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF]. We learn the model parameters by optimizing the cross-entropy loss between the predicted labels and the true labels using the Adam optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF]. We tune the following hyper-parameters: batch size ∈ {32, 42}, hidden units ∈ {16, 32}, learning rate ∈ {0.01, 0.001} and dropout coefficient ∈ {0, 0.2}. All implementations were carried out using Py-Torch Geometric [START_REF] Fey | Fast graph representation learning with PyTorch Geometric[END_REF] and the Cuda enabled GPU 12GB NVIDIA Tesla K80. Results reported in table 2 and figure 6 show that GIN models perfectly fit the data and report approximately the same results for predicting the anticancer activity of chemical compounds on both NCI-1 (GIN-1-layer gives an average accuracy of 83.0 ± 3.5% and GIN-MLP gives 83.4 ± 2.5%) and NCI-109 ( GIN-1-layer gives an average accuracy of 86.8 ± 6.0% and GIN-MLP gives 85.0 ± 6.7% accuracy) datasets, while Node2vec has a low and unstable performance on the two datasets ( average accuracy of 54.9 ± 1.2% on NCI-1 for lung cancer and 51.4 ± 2.3% on NCI-109 for ovarian cancer). Those results are maybe due to the fact that GIN models are deep while Node2vec is a shallow model. Subsequently, we validate the representation power of Node2vec by comparing it theoretically to the Graph Isomorphism Network (GIN) [START_REF] Xu | How powerful are graph neural networks? International Conference on Learning Representations[END_REF] which is proved to be as powerful as the Weisfeiler-Lehman(WL) test of isomorphism [START_REF] Weisfeiler | A reduction of a graph to a canonical form and an algebra arising during this reduction[END_REF] in capturing differences on graph structures and topologies. We use Node2vec and two variants of the GIN model to auspiciously predict the activity of chemical compounds against lung and ovarian cancer over NCI-1 and NCI-109 respectively. Chemical compounds predicted to be active against lung and ovarian cancer can then be decisive in the conception of pharmaceutical anticancer drugs.

Conclusion

In this work, we investigate the representation learning power of Graph Neural Networks on graph data structure. We explore two powerful models that operate on graphs (1) Node2vec and (2) Graph Isomorphism Network and derive that the latter is a deep and inductive extension of Node2vec. Thus, based on the aforementioned results we are able to apprehend the representation power of Node2vec. Predicting effectively the activity of chemical compounds against cancer helps pharmaceutical drugs producers in the fabrication of anticancer medicines. We exploit Node2vec and GIN representation power over graph data to predict the activity of chemical compounds against lung and ovarian cancer and achieve remarkable results using GIN variants. As a future work, we would like to define a new specially-designed deep learning framework for representation learning on chemical compounds modelled as graphs. We also plan to go beyond neighbor aggregation and adopt attention mechanisms to learn the connective strength between nodes.

Fig. 1 .

 1 Fig. 1. Embedding learning. Each node v in the graph is mapped to a d-dimensional embedding vector ev which summarizes the local information of the nodes around v.

Fig. 2 .

 2 Fig. 2. For every node (i.e source node ), matrix factorization approaches learn embedding vectors of source nodes based on information from local neighbors (i.e 1-hope neighbors) or distant neighbors.

Fig. 4 .

 4 Fig. 4. Unsupervised representation Learning on graphs. starting from a source node u defined by its corresponding column in the adjacency matrix  = A + IN , we update the embedding of u by maximizing the probability of visiting the neighbors v ∈ N (u) generated by a random walk. The embedding we are looking for are defined by the parameters W and W of the model.

Fig. 5 .

 5 Fig. 5. Graph representation of a chemical compound The left figure represents a chemical compound and the right figure represents its corresponding graph representation where C refers to carbon atom, H to Hydrogen atom and N to Nitrogen atom.

Fig. 6 .

 6 Fig. 6. Training set performance of GIN-MLP compared to Node2vec on NCI-1 and NCI-109 datasets over 100 epochs

Table 1 .

 1 Datasets statisticsThe table summarizes the number of graphs, the average number of nodes, the average number of edges and the number of node labels in NCI-1 and NCI-109 datasets.

	Datasets		Dataset Statistics
		Graphs Nodes(avg.) Edges(avg.) Node Labels
	NCI-1	4110	29.8	32.3	37
	NCI-109	4127	29.6	32.1	38

Table 2 .

 2 Binary-classification results(%) The table reports the average and standard deviations of accuracies on binary-graph classification

	Datasets	Representation Learning Models
		Node2vec[12] GIN-1-Layer [25] GIN-MLP [25]
	NCI-1	54.9 ±1.2	83.0 ±3.5	83.4 ±2.5
	NCI-109	51.4 ±2.3	86.8 ±6.0	85.0 ±6.7