dataesr/scientific-paragraphs-categorization huggingface dataset:
In addition, incorporating biological functionality from different molecular layers, such as RNA, proteome, and metabolome results, can boost the power of genetic mapping. In this study, we conducted an integrative analysis of GWAS and transcriptomic profiling for NSCLC using network building based on an algorithm that searches for direct interactions from a high-quality, manually curated database of genetic and physical interactions. To evaluate the identified networks, we repeated network building from a large set of randomly generated gene lists for distributions of network complexity and gene ratio. We also used hub genes identified from significant networks for targeted sequencing and further validation in the Transdisciplinary Research in Cancer of the Lung-International Lung Cancer Consortium (TRICO-ILCCO) GWAS meta-analysis.
Study Population This study was based on data derived from 1000 NSCLC cases and 1000 cancer-free controls, frequency-matched by age (±5 years), gender, and smoking status (by packyears) as previously described  [19]. All cases were recruited at Massachusetts General Hospital (MGH) from 1992 to 2004, were N 18 years old, and had newly diagnosed, histologically confirmed primary NSCLC. Controls were healthy, nonblood-related family members and friends of patients with cancer or with cardiothoracic conditions undergoing surgery. Histological classification was done by two staff pulmonary pathologists at MGH according to the International Classification of Diseases for Oncology (ICD-O3). For histology analysis, the following codes were used: adenocarcinoma, 8140/3, 8250/3, 8260/3, 8310/3, 8480/3, and 8560/3; large cell carcinoma, 8012/3 and 8031/3; squamous cell carcinoma, 8070/3, 8071/3, 8072/3, and 8074/3; and other non-small cell carcinomas, 8010/3, 8020/3, 8021/3, 8032/3, and 8230/3.
GWAS Dataset DNA was extracted from peripheral white blood cells using standard protocols and was genotyped using the Human610-Quad BeadChip (Illumina, San Diego, CA). Before association tests, we conducted a systematic quality evaluation of raw genotyping data according to a general quality control (QC) procedure described by Anderson et al. [20]  Briefly, unqualified samples were excluded if they fit the following QC criteria: (i) overall genotype completion rates b95%; (ii) gender discrepancies; (iii) unexpected duplicates or probable relatives (based on pairwise identity by state value, PI_HAT in PLINK N0. 185); or (iv) heterozygosity rates N6 standard deviations from the mean. Unqualified SNPs were excluded if they fit the following QC criteria: (i) overall genotype completion rates b95%; (ii) gender discrepancies; (iii) unexpected duplicates or probable relatives (based on pairwise identity by state value, PI_HAT in PLINK N0.
185); (iv) heterozygosity rates N6 standard deviations from the mean; or (v) individuals were non-Caucasians (using the HapMap release 23, including JPT, CEPH, CEU, and YRI populations as reference). Unqualified SNPs were excluded if they fit the following QC criteria: (i) not mapped on autosomes; (ii) call rate b 95% in all GWAS samples; (iii) MAF b 0. 01; or (iv) genotype distributions deviated from those expected by Hardy-Weinberg equilibrium (p b 1. 0 × 10 -6 ). After quality evaluation, we had a dataset of 984 cases and 970 controls with 543,697 autosomal SNPs for epistasis analysis.
External Transcriptomic Data Two transcriptomic datasets of NSCLC were selected and raw data were downloaded from Gene Expression Omnibus (GEO) database (https://www. ncbi. nlm. nih. gov/geo): GSE10072 and GSE18842. GSE10072 focused on lung adenocarcinoma and included 180 paired adenocarcinoma and non-affected tissue samples  [23]. GSE18842 included 91 samples of mainly squamous cell carcinoma  [24]. GSE10072 and GSE18842 were generated from fresh snap-frozen samples from surgical resection and profiled on Affymetrix Human Genome U133 array (Affymetrix, Santa Clara, CA).
Transcriptomic Data Analysis Although FFPE profiles and external data were generated from different platforms, we used DNA-Chip Analyzer 2006 (dChip, http:// www. dchip. org) software, which applied an invariant set of genes for normalization and calculation of expression values across all microarrays, to normalize raw microarray signals. This analysis assumed that a subset of genes had constant expression among all cell subtypes  [25]  Only paired tumor and non-affected tissue samples were used in the analysis, including 18 FFPE adenocarcinomas, 8 FFPE squamous cell carcinomas, 33 snap-frozen adenocarcinomas (GSE10072), and 32 snapfrozen squamous cell carcinomas (GSE18842).
Integrated Analysis by Network Building Gene lists from GWAS and/or transcriptomic profiling were uploaded into MetaCore GeneGo database (https://portal. genego. com, Thomson Reuters, New York, NY) for network building which has N1. 7 million molecular interactions, 1600 pathway maps, and 230,000 gene-disease associations  [26]. MetaCore is an integrated online software suite for functional analysis of omics data that is based on a high-quality, manually-curated database of molecular interactions, molecular pathways, gene-disease associations, chemical metabolism, and toxicity information. We used direct interaction algorithms on the MetaCore platform to build gene networks consisting only of uploaded genes and their direct interactions, without adding other genes/objects from the GeneGo database.
Network Evaluation by Randomly Generated Gene Lists SNPs of GWAS data and probes of transcriptomic profiling mapped to a total of 24,847 genes. From these genes, we randomly generated 6 sets of gene lists, with each set containing 100 gene lists, for a total of 600 random gene lists. Lists from each set had the same number of genes, but lists from different sets had different numbers of genes-either 50, 100, 200, 300, 400, or 500 genes. Each random gene list was individually uploaded into MetaCore GeneGo database for network building, and network parameters, including network complexity and gene ratio, from the largest networks were recorded. Quantile regression at 95th percentile was performed to estimate the 95% confidence interval.
Targeted and Whole Exome Sequencing Targeted and whole exome sequencing was performed at the Center for Inherited Disease Research. Ninety-nine custom regions targeted for a total of 17. 26 Mb of custom content was captured and sequenced. BAM files were created by aligning FASTQ files to GRCh37 and joint sample variant calling and variant site filtering was performed. Genotypes for biallelic SNPs were further refined using CalculateGenotypePosteriors and allele frequency information from 1000 genomes phase 3 data as well as the Exome Aggregation Consortium data. Further details were described in Supplementary Methods.
Statistical Analysis Clinical characteristics were described as mean ± standard deviation (SD) for continuous variables or n (%) for categorical variables. Student's t-test or Fisher's exact test was used for comparison between groups for continuous or categorical variables, respectively. We carried out gene-based analysis on GWAS data and targeted sequencing data using SKAT-O  [27]. SKAT-O aggregates weighted variance-component score statistics for each SNP/SNV within a set using a kernel function and tests for associations between groups of SNPs/ SNVs and a phenotype while adjusting for relevant covariates  [28]. For GWAS data, initially all SNPs that passed QC were mapped to human genes within ±20 kb regions based on information curated in the RefSeq database(NCBI build GRCh37. p13). Separate analyses were conducted individually for all SNPs and rare SNVs with MAF b0.
01. Models were adjusted for age (years), gender, smoking status, and top significant eigenvectors. SKAT-O analyses were carried out using the SKAT package (R v. 2. 13. 0). In meta-analysis of GWAS datasets within TRICL-ILCCO, we combined imputed genotypes from 13,479 lung cancer cases and 43,218 controls undertaken by the previous TRICL-ILCCO GWAS  [22, 29, 30]. We excluded poorly imputed SNPs defined by imputation quality Rsq b0. 3 or Info b0. 4 for each study and conducted fixed effects meta-analysis with inverse variance weighting and random effects meta-analysis from the DerSimonian-Laird method  [31]. We also generated an index of heterogeneity (I2) and p-value of Cochran's Q statistic to assess heterogeneity in meta-analyses. We only considered SNPs with MAF N0. 005 and that showed little evidence for effect heterogeneity between studies (Cochran's Q statistic p N.
Initial Screening of GWAS and Transcriptomic Profiling Data We produced a genomic dataset of germline polymorphisms (GWAS genotyping data on the Illumina 610 Quad platform) and a transcriptomic dataset of tumor and non-affected tissue (genomewide expression profiling on Illumina DASL HT platform) from NSCLC samples  [29]. The GWAS dataset included 543,697 single nucleotide polymorphisms (SNPs) from 984 NSCLC cases and 970 healthy controls after quality assessment. Transcriptomic profiling was carried out on formalin-fixed paraffin-embedded (FFPE) paired samples of tumor and non-affected tissues from 30 NSCLC cases, including 18 adenocarcinoma, 8 squamous cell carcinoma, and 2 unclassified cases. Eighteen NSCLC cases had both GWAS and transcriptomic data. Patient characteristics are described in Table  1.
Instead of analyzing individual SNP, we conducted gene-based analysis of GWAS data using the optimal unified sequence kernel association test (SKAT-O) method and applied gene lists in the subsequent network analysis  [32]. Among 21,981 mapped genes of GWAS data, there were 103 genes with p b. 005 (Supplementary Table  S1 ), 232 genes with p b. 01, and 1007 genes with p b. 05. The top genes associated with risk of developing NSCLC were HYKK (also known as AGPHD1, 15q25, p = 2. 30 × 10 -6 ), CLPTM1L (5p15, p = 3. 54 × 10 -5 ), CHRNA3 (15q25, p = 6. 77 × 10 -5 ), and DNAJC16 (1p36. 1, p = 7. 12 × 10 -5 ), with 3 genes located within the two previously identified risk loci at 5p15 and 15q25  [29]. We also screened transcriptomic data for differentially expressed genes with N2 fold changes (FC) between tumor and non-affected lung tissues of 18 pairs of FFPE samples (lung adenocarcinoma, including bronchioloalveolar carcinoma) obtained by surgical biopsy.
1b ). Moreover, the network was centered on important genes belonging to the hypoxia-inducible factor (HIF) family, including hub genes HIF1A, ARNT (also known as HIF1B), and EPAS1 (also known as HIF2A), which are transcription factors that respond to changes in available oxygen in the cellular environment  [33]. We therefore named this the HIFs-EGFR-HDAC4-TERT network. Similarly, using top-ranked genes from a GWAS dataset containing only lung adenocarcinomas (597 cases and 970 controls), we built a significant network from a combined list of GWAS genes (90 genes with p b. 005), transcriptomic profiling genes (75 genes with p b. 001), and core genes.
Network Evaluation by Random Gene Lists To examine the possibility that significant networks were formed by random chance of the increased number of genes from combined gene lists, we conducted a series of network analyses using randomly generated lists with different numbers of genes (total list: n = 600) and evaluated the largest network built from each random list. As the number of genes for network building increased, we more frequently observed an increased number of genes and connections within networks, resulting in increased gene ratios. However, the complexity of these networks remained relatively unchanged (data not shown). Networks were better evaluated in a two-dimensional space of network complexity and gene ratio than any individual parameter we tested. As shown in Fig. 2 , the large complex network built from combined gene lists was located far above from a 95% upper-tail conference interval (95% CI) curve, whereas most networks from individual gene lists or combinations of any two lists were located either under or around the 95% CI curve.
Network Validation by External Transcriptomic Data An external transcriptomic dataset (GSE10072) was generated on a different platform (Affymetrix Human Genome U133A array) using 33 fresh-frozen pairs of tumor and non-affected tissues from NSCLC adenocarcinomas collected in the Lombardy region of Italy  [23]. We identified 85 top-ranked, differentially expressed genes (FC ≥ 2; p b 10 -16 ). Similar to the FFPE dataset, a significant network could only be built from combined lists of top-ranked GWAS genes, top-ranked transcriptomic profiling genes, and a core list of 29 non-overlapping genes (Supplementary Fig. S2 ). This network had a complexity and gene ratio comparable with that from the FFPE dataset and also had common hub genes (≥5 connections within the network) shared with the FFPE dataset (Table  3 ).
Moreover, we identified 88 common, differentially expressed genes (FC ≥ 2; p b. 05) between FFPE and GSE10072 datasets, with 85 genes (97%) having expression changes in the same direction between tumor and non-affected tissues (Supplementary Table  S3 ). Although only a simple and small network could be built from this common transcriptomic list, a significant network with the same hub genes, including HIF1A, FOS, HDAC4, and EDN1, could be built by direct combination with the GWAS list (Table  3 ).
05/13 genes = 0. 0038). Further, there was no significant aggregation of variants with moderate to high functional impacts (http://useast. ensembl. org/info/genome/ variation/predicted_data. html) in either NSCLC cases or controls in the other 12 hub genes (Supplementary Table  S5 ). In EPAS1, sequencing identified 2061 SNVs, including 1617 rare/low-frequency SNVs [minor allele frequency (MAF) ≤ 0. 01] and 36 SNVs with moderate-high functional impact. One common SNV (SNP: rs12614710) located within the first intron and identified by sequencing reached genome-wide significance (MAF = 0. 45; OR = 1. 50; 95% CI: 1. 31-1. 72; p = 7. 75 × 10 -9 ) (Fig. 3 ). We further validated SNP rs12614710 in a much larger GWAS dataset using meta-analysis. A fixed-effect model was applied to estimate pooled effects of each SNP using the TRICL-ILCCO GWAS dataset, which included 13,479 lung cancer cases and 43,218 controls (Supplementary Table  S6 )  [34].
Discussion In this study, we conducted an integrative analysis of multi-omic datasets of NSCLC to assess associations beyond the genome and transcriptome. A large, complex, and robust gene network containing well-known lung cancer-related genes, including EGFR and TERT, was identified for lung adenocarcinoma from the combined gene lists. However, the framework of this network was built by key members of the HIF gene family. Subsequent sequencing of network hub genes within a subset of consortium samples revealed a SNP (rs12614710) in EPAS1 associated with NSCLC that reached genome-wide significance based on whole exome sequencing data. Although this SNP was not covered in any GWAS dataset, we used imputed data to find that this SNP is borderline significant in the entire TRICL-ILCCO GWAS dataset.
8 × 10 -8 ; rs9679290, p = 5. 75 × 10 -8 ; rs4953346, p = 4. 09 × 10 -14 ) and prostate cancer  [40] [41] [42]. In a small study of 346 NSCLC patients and 247 controls from a Japanese population, SNP rs4953354 was associated with increased risk of lung adenocarcinoma (OR = 1. 80; 95% CI, 1. 16-2. 79; p =. 008) [42]. In our TRICL-ILCCO GWAS dataset of 13,479 cases and 43,218 controls  [34] , all of these SNPs were significantly associated with NSCLC (rs11894252, p =. 043; rs9679290, p =. 0011; rs4953346, p =. 0015; rs4953354, p =. 025). All previously reported SNPs are located in the first intron, except for rs4953354, which is located in the third intron. In addition, some SNPs are not correlated  [40] , including rs12614710 from our sequencing project. Moreover, bioinformatic analyses using Genome Browser (http:// genome.
ucsc. edu) suggest that most of these first intron SNPs are located in histone mark H3K27Ac, which is defined by a ChIP-seq assay related to enhanced gene transcription  [43]. Further, analysis of ChIP-seq datasets from ENCODE identified binding sites and binding activities for C/EBP-β, AP-1, and MYC families of transcription factors in many cancer cell types within the first intron of EPAS1. Further, the A allele of rs13419896 is associated with enhanced EPAS1 expression and poor prognosis of 76 NSCLC patients  [44]. It is likely that genetic polymorphism of EPAS1 may lead to varied gene expression through either changes in binding sites and binding activities for certain transcription factors or modification of histone epigenetic regulation. In a study of chronic obstructive pulmonary disease, hypermethylation of EPAS1 is correlated with decreased EPAS1 expression and is significantly associated with disease severity  [45].
Although GWAS has provided useful insights into the genetic architecture of complex diseases, there is weak evidence for how GWAS findings improve understanding of molecular pathways involved in disease, thus bringing post-GWAS challenges to the characterization of molecular data. Therefore, it is important to assess how diverse omic datasets at different biological levels can be integrated to exploit the full potential of information to identify causal genes and networks, regulatory genes and networks, and predictive markers for complex traits. Using direct interaction algorithms for network building, we successfully conducted an integrated study of multi-omic data for exploration beyond GWAS. This approach implemented a stringent criterion of only searching for direct gene-gene interactions within a manually curated database (MetaCore, https://portal.
genego. com), while using less strict p-value cut-offs to select gene lists from different omic datasets. Thus, we could explore less significant genes, which often do not reach genome-wide significance, in omic datasets. The underlying hypothesis is that, because genes are highly networked and coordinated and do not act alone, polymorphisms of several genes in one biological process might not reach genome-wide significance. Initially, we could not find any meaningful networks (with relatively higher complexities and gene ratios) from gene lists selected from single omic datasets. This result is reasonable in that a list of top genes with small p-values is more likely to contain discrete genes, especially for transcriptomic data, as differentially large expressions were often downstream target genes in a transcription pathway and might not have direct interactions.
Only combined gene lists from both GWAS Fig. 1. Networks built from GWAS and transcriptomic profiling data. a, Significant network built from combining top-ranked GWAS genes (103 genes with p b. 005), topranked transcriptomic profiling genes (75 genes with p b. 001), and a core list of genes (37 non-overlapping genes out of 46 core genes). Blue indicates genes only identified from transcriptomic profiling; green indicates genes only identified from GWAS; and red indicates genes from the core list. b, Network of 13 hub genes (≥5 connections within the network). and transcriptomic datasets plus a core list of common genes could build a large, complex HIFs-EGFR-HDAC4-TERT network. However, with even less stringent p-value cut-offs, we built a large network from genes selected from a single omic dataset.
A supergene network always was large in size but low in complexity. The HIFs-EGFR-HDAC4-TERT network had a moderate size but high complexity without supergenes. From the distribution of network properties of 600 randomly selected gene lists plotted for network size against network complexity, the HIFs-EGFR-HDAC4-TERT network was above the 95% CI curve, suggesting that this network was not randomly built. We further validated the HIFs-EGFR-HDAC4-TERT network generated from multi-omic datasets by using different transcriptomic datasets. The GSE10072 dataset was from a study of gene expression signatures of cigarette smoking and its role in lung adenocarcinoma development and survival, and it contained 33 fresh-frozen pairs of tumor and nonaffected tissues from NSCLC adenocarcinomas  [23]. Combined with our GWAS dataset, which contained a majority of NSCLC adenocarcinomas, and a core gene list, we also built a network with network size and complexity comparable to the HIFs-EGFR-HDAC4-TERT network and several shared hub genes.
Different combinations of the GWAS dataset, transcriptomic dataset, and core gene list provided similar results. We also analyzed a transcriptomic dataset of squamous cell carcinoma, including our 8 pairs of tumor and non-affected tissues from FFPE samples and 32 pairs of fresh-frozen samples from GSE18842  [24] , and found that no complex networks except several supergene Fig. 2. Distribution of networks built by randomly selected genes. From a total of 24,847 genes, we randomly generated 6 sets of gene lists, with each set containing 100 gene lists. Each set contained 50, 100, 200, 300, 400, or 500 genes, denoted by R50, R100, R200, R300, R400, and R500, respectively. Networks were built from each gene list using MetaCore GeneGo database and were used to calculate a 95% CI curve. GWAS denotes network from GWAS data; FFPE denotes network from transcriptomic profiling data; and GWAS+FFPE+core denotes networks built from combining GWAS data, transcriptomic profiling data, and a core gene list.
networks could be built (data not shown). This might be because a majority of our GWAS samples were lung adenocarcinomas (60. 7%), and squamous cell carcinomas only accounted for ~22% of samples. Meanwhile, no identified network from squamous cell GWAS transcriptomic datasets suggests that our integrated approach of multi-omic data was sensitive to tumor histology. Sequencing of hub genes identified a new locus in EPAS1 that reached genome-wide significance and was validated in the largest lung cancer consortium, providing additional evidence that the HIFs-EGFR-HDAC4-TERT network is associated with NSCLC adenocarcinoma. Further, a recent study reported that EPAS1 could specifically bind to tyrosine kinase inhibitor (TKI)-resistant T790 M EGFR in NSCLC cell lines and enhance amplification of MET  [46].
These findings suggest that EPAS1 is a key factor in EGFR-MET crosstalk in conferring TKI resistance in NSCLC cases and provide in vitro support of the HIFs-EGFR-HDAC4-TERT network. At present, FFPE samples, which represent the greatest stock of archived disease entities, are limited mainly to investigations of a small number of genes using quantitative real-time PCR or global micro-RNA profiling, which is much more stable than mRNA  [47, 48]. The main reason for this restriction is that RNA is often altered and degraded within FFPE samples from the impact of collection and storage  [49, 50]. In our study, from the top 808 differentially expressed genes from FFPE samples, we identified 88 common, differentially expressed genes between FFPE and the GSE10072 dataset, with 85 genes (97%) having expression changes in the same direction between tumor and non-affected tissues.
This common gene list could be used as a transcriptomic list to build a network containing HIFs without incorporating a core gene list. These results demonstrate that FFPE samples could generate a transcriptomic profile for integrated analysis, as we found similar networks with shared hub genes as compared to freshfrozen samples. During the analyses, we evaluated different network results by exploring different p-value cut offs for the selection of genes from different omics datasets. We didn't find any significant networks as all networks were under the 95% CI curve in random gene list evaluation. We noticed there was a limitation that we hadn't control this level of multiple comparisons. Currently, we didn't have a proper analytic method to adjust the multiple comparisons. Therefore, we used external data to evaluate the network results.