1. Study search and selection

Publicly available datasets were collected from GEO [2], ArrayExpress [3], and TCGA [4]. A systematic search of studies was conducted between September 2018 and January 2019 following the preferred reporting items for systematic reviews and meta-analyses (PRISMA) guidelines [5]. Several keywords were employed in the search, including lung adenocarcinoma (ADC), non-small-cell lung carcinoma (NSCLC), Homo sapiens and excluding cell lines. 11 variables were taken into account to evaluate each individual study, including clinical characteristics of the patients (e.g., sex and smoking habit) and experimental design (e.g., sample size and sample extraction source). The final inclusion criteria were:

• Sex, disease stage and smoking habit variables were registered
• RNA was extracted directly from human lung biopsies
• Both normal and lung ADC samples were available
• Patients had not undergone treatment prior to biopsy
• Sample size was > 3 for case and control groups in both sexes

Finally, normalized gene expression data of 6 array NSCLC datasets (GSE10072, GSE19188, GSE31210, GSE32863, GSE63459 and GSE75037) and counts matrix of 3 RNA-seq NSCLC datasets (GSE81089, GSE87340 and TCGA-LUAD) were retrieved.

2. Individual transcriptomics analysis

Individual transcriptomics analysis consisted of three steps: pre-processing, differential expression analysis, and functional enrichment analysis (Figure 1a). Data preprocessing included the standardization of the nomenclature of the clinical variables included in each study, normalization of RNA-seq counts matrix and exploratory analysis. RNA-seq counts were pre-processed with the edgeR [6] R package using the trimmed mean of m-values (TMM) method [7]. Annotation from probe set to Entrez identifiers from National Center for Biotechnology Information [8] database and gene symbol was carried out with the biomaRt [9] R package. When dealing with duplicated probe-to-entrez mappings, the median of their expression values was calculated. The exploratory analysis included unsupervised clustering and principal component analysis (PCA), to detect patterns of expression between samples and genes, and the presence of batch effect in each study (Figure 1b). Differential expression analyses were performed using the limma [10] R package. To detect differentially expressed genes in the progression of lung ADC between men and women, the following contrast was applied:

(ADC.W - Control.W) - (ADC.M - Control.M) where ADC.W, Control.W, ADC.M and Control.M correspond to lung adenocarcinoma affected women, control women, lung adenocarcinoma affected men and control men, respectively. Paired samples design was implemented and tobacco consumption was included as a batch effect to reduce its impact on data. P-values were calculated and corrected for False Discovery Rate (FDR) [11]. This comparison allows us to detect which genes and functions have been altered by the disease and have higher or lower activity in women when compared to men. Significant functions and genes were considered when FDR < 0.05. Functional enrichment analyses were performed using Gene Set Enrichment Analysis (GSEA) implemented in the mdgsa [12] R package. P-values were, again, corrected for FDR. For functional annotation, two functional databases were used: Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway database [13], and Gene Ontology (GO) [14]. GO terms were analyzed and propagated independently for each GO ontology: biological processes (BP), molecular functions (MF) and cellular components (CC). Those annotations excessively specific or generic were filtered, keeping functions with blocks of annotations between 10 and 500. Intersections within groups were analyzed with UpSet plots [15] (Figure 1c).

3. Functional meta-analysis

Functional GSEA results were integrated into a single functional meta-analysis [16] implemented with mdgsa and metafor [17] R packages. Meta-analysis was applied under the DerSimonian & Laird random-effects model [18], taking into account individual study heterogeneity. This model takes into account the variability of individual studies by increasing the weights of studies with less variability when computing meta-analysis results. Thus, the most robust functions between studies are highlighted. A total of 6.330 GO BP terms, 745 GO CC terms, 1193 GO MF terms and 212 KEGG pathways were evaluated. P-values, FDR corrected p-values, the logarithm of odds ratios (LOR) and 95% confidence intervals (CI) were calculated for each evaluated function. Functions and pathways with FDR < 0.05 were considered significant, and both funnel and forest plots were computed for each one (Figure 1d). Funnel plots were analyzed to assess for possible biased results, where LOR represents the effect size of a function and the standard error (SE) of the LOR serves as a study precision measure [19].