dataesr/scientific-paragraphs-categorization huggingface dataset:
The data have been continuously collected using massive financial and scientific efforts. Now, the question arises if this resource can also be used to support clinicians in making rapid and accurate assessments leading to tailored treatments for individual cancer patients. Therefore, methods that can incorporate newly generated "omics" data (further addressed as "investigation dataset") in large public reference datasets are required for improved diagnostics of new patients. Here we propose an ICA-based method to combine reference and investigation datasets for patient diagnostics and detailed inspection of biological processes in cutaneous melanoma  (Hayward, Wilmott et al. , 2017). Melanoma arises through the malignant transformation of melanocytes and presents a very aggressive form of skin cancer with increasing global case numbers (https://www.
, 2013 , Zhang, Dutton-Regester et al. , 2016). Nevertheless, according to mutational hotspots, melanoma patients have recently been allocated to four groups: BRAF, NF1, RAS or Triple WT (wild type), while the analysis of gene expression data resulted in three functional groups: "immune subclass", "keratin subclass" and "MITF-low subclass", which according to the authors have implications for patient survival  (Cancer Genome Atlas, 2015). Interestingly, the majority of primary melanomas belonged to the "keratin subclass" having a worse prognosis than the other two subclasses. In this study, we used the SKCM TCGA cohort with over 470 patients diagnosed with cutaneous melanoma as the reference dataset. Two layers of sequencing data were considered and integrated: mRNA and microRNA (miRNA). The investigation dataset included a small cohort of three primary melanoma tumours and two controls: matched cancer patient-derived normal skin and NHEM cells (normal human epidermal melanocytes).
First, for the reference cohort, we demonstrated that ICA deconvolution was successfully applied to classify patients based on their tumour subtypes and to build the hazard score that predicts patient survival. The hazard score was then tested using an independent validation cohort of 44 patients, obtained by microarray gene expression technology. The strong difference between reference RNA-seq data and microarray-derived validation datasets was overcome by our method. Next, the investigated dataset was studied in depth and key processes involved in cancer aetiology were detected: immune response and inflammation, angiogenesis, selfsufficient cell proliferation and other. We show here that consensus ICA can integrate data from different biological sources, time points and platforms to predict clinically important characteristics of cancer in a bias-free, unsupervised and potentially automatable fashion, suggesting consensus ICA as a useful module of future clinical "oneclick" support systems.
ICA of combined data sets can remedy technical biases Reference, validation and investigation datasets. Using publicly available data sets for diagnostic/prognostic purposes on new patients is often hampered by the fact that different technical platforms, data formats and analysis pipelines are used to generate the reference data and the patient samples to be investigated. Thus, reference and newly generated data do not come from the same distribution, which often introduces strong technical biases despite different normalisation schemes. It has recently been shown that ICA, when applied to heterogeneous datasets affected by technical factors, can identify such biases and single out one or several components to account for them  (Biton et al. , 2014 , Taroni & Greene, 2017). As soon as the technical bias is captured by an independent component, it is isolated from the core data, which can then be analysed in a more accurate way.
In this study, graphically outlined in Figure  1 , we used public TCGA data as the reference dataset, published microarray data  (Bogunovic, O'Neill et al. , 2009)  as a validation set and an investigation set based on clinical samples -three primary melanoma tumours (denoted P2PM, P4PM, P6PM) and two control samples (matched normal skin P4NS and NHEM cell line) as described in Table  1. We developed a method to assign a prognostic hazard score based on the reference set and validated our approach on an independent public microarray data set, followed by testing the method on clinical samples. The underlying ICA was done on single or combined datasets (reference only, reference/validation and reference/investigation) by the developed consensus ICA method.
MICs were reoriented based on the negative sign of correlation with RICs. ICA identifies technical biases. The combined reference/investigation dataset profiled by RNA-seq is presented in the space of two first principal components (Figure  2A ) and of two selected independent components (Figure  2B ). The two principal components included 33% of total variability and mainly reflected technical effects: PC1 was linked to the RNA-seq library size (data not shown) and PC2 segregates reference and investigation data. Among all RICs, the components that reflected data clustering according to gender (RIC3) and sample type (primary or metastatic, RIC5) were chosen. Both components were detected with high reproducibility (R 2 RIC3 = 0. 996, R 2 RIC5 =0. 993). The investigation data were clearly integrated within the reference data and showed reasonable clustering in Figure  2B.
P2PM and P4PM samples originated from primary tumours of female patients, P4NS normal skin corresponding to P4PM samples and NHEM was extracted from eyelid skin of a female. P6PM was a primary tumour of a male patient. Functional analysis showed that genes involved in RIC5 participate in keratinocyte-specific functions and thus weights of RIC5 could be used as a marker of keratinocyte presence. Indeed, the vast majority of metastatic samples had low values of RIC5 weights, while primary tumours showed high values. NHEM (pure melanocytes) are devoid of keratinocytes and therefore clustered with metastatic tissues. An even stronger normalisation effect of ICA was observed for miRNA data. As mentioned above, the techniques used for miRNA detection and quantification were different: miRNA-seq for the reference dataset and whole miRNome qPCR arrays for the new samples.
Combined HS showed high significance (Cox log rank p-value = 2. 2e-13) for the TCGA dataset. In order to validate the proposed hazard scoring approach on an independent cohort of patients, we repeated the ICA on the combined TCGA RNA-seq data (reference set) and microarray data E-GEOD-19234 (additional validation set). Notably, the metagenes of the new decomposition were strongly correlated with the metagenes from the initial reference-set; only ICA and 44 of them showed R 2 >0. 5 (95% prediction interval for R 2 values between metagenes of the same decomposition was: 0<R 2 <0. 004). The components that showed a significant (adj. p-value<0. 05) link to survival on the reference set were then used to compose HS for the validation data and also showed significant prognostic properties (log rank p-value of 0.
0013); Kaplan-Meyer plot shown in Figure  3. The developed HS separated patients with low hazard (only one death among 7 patients, blue line in validation set, Figure  3 ) from the group of patients with a high hazard score. For the three primary melanoma samples from the investigation set, the calculated HS was the lowest for the P2PM sample (HS=1. 16) and the highest for P6PM (HS=1. 92). This was in agreement with clinical observations, as patient P6 suffered from a very aggressive form of melanoma and deceased shortly after sample collection. From the quantitative results obtained from the validation dataset and qualitative differences observed for the investigation dataset, we concluded that weights of independent components can be combined into a hazard score, suitable to predict patient survival.
63, "B cells memory" r=0. 58, "B cells naïve" r=0. 57, Supplementary Results, and also showed the highest correlation with B cell profiles measured in single cells (Supplementary Figure  S3 )). RIC27 showed a very similar collection of enriched gene sets. These two components most probably represent naïve and activated B-cells. Functionally, RIC28 was linked to inflammatory responses to wounding (adj. p-value=6. 3e-22), neutrophil degranulation (adj. p-value=1. 3e-7), TNF (adj. p-value=4. 7e-8) and IL1-mediated signalling pathways (adj. p-value=2. 2e-9); RIC37 was connected to interferon signalling (adj. p-value = 5. 1e-22) whose metagenes were also reciprocally correlated with M5_INTERFERON of the Biton dataset  (Biton et al. , 2014). Components RIC74 and RIC79 contained a very limited number of topcontributing genes, but both were significantly linked to survival (p-values of 1.
p-value < 2e-15 for all mentioned components) and keratinisation (adj. p-value < 1. 4e-10). Enrichr suggested that the signals of these components are specific to skin (adj. p-value < 1e-50). The dataset contained 48 keratins and many of them were observed among the top-contributing genes: 20 for RIC5, 28 (RIC7), 30 (RIC19) and 13 (RIC31). RIC5 and RIC7 were negatively correlated with survival, which is in concordance with previous observations (Cancer Genome Atlas, 2015). MIC9 with the skin-specific miR-203  (Yi, Poy et al. , 2008) , was linked to RIC5, RIC7 and RIC31. Furthermore, several components (RIC4, RIC16, MIC11 and MIC14) were connected to the activity of melanocytes. Top-contributing genes of RIC4 were enriched in the melanin biosynthesis process (adj. p-value 1. 2e-5) and Enrichr linked these genes to melanocytes (adj.
Increased cell proliferation was linked to survival (Cox p-value=3. 0e-2). In the investigated samples, the highest weight was observed for the most aggressive tumour P6PM and the lowest value for normal skin P4NS. Several RICs showed linkage to neural tissue. As an example, both positive and negative topcontributing genes of RIC6 were linked to brain in the ARCHS4 tissue sets of Enrichr (both adj. pvalues < 1e-33). This component was as well associated with patient survival (Cox p-value=5. 5e-3). The component indicates the ability of melanoma cells to show expression patterns specific for cells of the neural crest of human embryos and can be linked to motility of malignant melanocytes. MiRNA component MIC1 showed an interesting bi-modal distribution in the reference dataset (see two clusters in Figure  2D ) and was strongly linked to patient survival (Cox p-value=9.
4e-4), suggesting two subgroups of melanoma patients with different prognosis. This component most probably was linked to regulation of EMT, as many miRNA positively or negatively influencing the component are known to be EMT regulators or linked to metastasis formation: miR-551, miR-206, miR-34a, miR-1269, miR-205, miR-876, miR-301b, and miR-365a. Based on our analysis of the reference TCGA dataset, these miRNA listed in Supplementary Results can be further investigated as potential survival markers for melanoma patients.
The method directly works with the data from bulk samples without any preliminary assumption about the transcriptomes of the purified cell types. Among the components, one can expect to see not only those defined by "pure" tumours or stromal cells, but also those originating from tumour/stroma interactions including tumour-induced stromal cell reprogramming. One example of such interactions is angiogenesis, further discussed below. We implemented a robust consensus ICA method and applied it to several datasets from patients with cutaneous melanoma (SKCM). These included (a) a large cohort of SKCM patients from TCGA, used as reference set; (b) an independent cohort of 44 patients with publicly available microarray mRNA data and (c) 5 in-house clinical investigation samples: 3 primary melanomas, a normal skin sample and a normal melanocyte cell line (NHEM).
Both mRNA and miRNA datasets were obtained for the reference and investigation samples. Despite the fact that different techniques were used for data acquisition, ICA was able to identify common signals in the datasets and properly allocate the new samples within the reference set (Figure  2 ). This was particularly evident for miRNA data where the reference set was obtained by small RNA-seq and the new samples by qPCR arrays with PCA showing a strong difference between these two datasets. With ICA, technical biases in the data were isolated within several components and thus separated from biologically relevant signals leading to a better and more correct characterisation of the samples. The fact that ICA should be rerun for every series of new samples could be considered as a drawback of our approach.
However, recalculation of the components does not require supervision and could be done automatically. In the case when investigation and reference datasets come from the same distribution, one can use the matrix S obtained from the reference dataset in order to define the weights (M) for the samples forming the investigation dataset (Eq. 1). However, in reality, the variability in the data requires recalculation of the components for the new investigated samples. We demonstrate here that when analysing data from melanoma patients, the weights of independent components can be used as predictive features of patient subgroups and can also be linked to patient survival. While the ICA-based feature extraction method has been previously discussed (e. g. (Aziz, Verma et al. , 2016 , Teschendorff et al.
Preparation of the reference and validation datasets Expression data. As reference dataset, we used the open-access TCGA skin cutaneous melanoma   (Nazarov, Muller et al. , 2017) , we shifted and scaled the microarray data: the 5 th percentile of expression was used as the lowest meaningful signal and was subtracted from microarray gene expression. All negative values were set to 0. The data were then scaled to unify the 75 th percentile between reference RNA-seq and validation microarray data.
Clinical data. To explore the possibility of assigning clinical traits to TCGA samples, we compiled a small dataset based on public TCGA data with "fail-safe" items covering gender and sample type (primary tumour and metastatic). Additionally we added publication-based data for RNA-seq clustering (immune / keratin / MITF-low) (Cancer Genome Atlas, 2015) as this information has been claimed to be relevant for disease prognosis. Survival data were extracted by parsing related information out of publicly available individual clinical data files (XML files provided by GDC). We extracted and processed the information assigned to the tags: bcr_patient_barcode, days_to_death, vital_status, year_of_initial_pathologic_diagnosis, age_at_initial_pathologic_diagnosis, days_to_last_followup and person_neoplasm_cancer_status.
Preparation of the investigation dataset: clinical samples, data acquisition and analysis The investigation dataset, represented by RNA-seq and miRNA qPCR array data, is composed of primary tumour samples of three melanoma patients and two control samples (one matched normal skin and a healthy melanocyte cell line, NHEM). Sample annotation is presented in Table  1. Details of sample collection, preparation, transcriptome and miRNome analyses are described in Supplementary Methods. RNA-seq data for these samples are available by GEO accession number GSE116111 and Ct-values for all quantified miRNAs are available in Supplementary Table  S5. To harmonise miRNA annotation of qPCR arrays and TCGA-derived miRNA isoform read count data, we first re-annotated our qPCR arrays to miRNA version 21. To have comparable data between qPCR arrays and TCGA, we worked with miRNA isoform data referring to miRNA IDs, so that mapping of stem loop IDs to mature miRNA IDs was possible.
Data Analysis RNA-seq expression measures. Four metrics of gene expression were considered: raw counts, DESeq2-normalized counts  (Anders & Huber, 2010) , Fragments Per Kilobase of transcript per Million (FPKM) and Transcripts Per kilobase Million (TPM). All expression values were log 2 transformed. Raw gene expression data were represented by 60446 features, of which many were lowly expressed. In order to reduce the number of uninformative features, we applied soft filtering, by the maximum expression level as described in  (Nazarov, Reinsbach et al. , 2013) : only genes that showed over 1000 counts in at least one sample of the reference TCGA SKCM dataset were considered (Supplementary Figure  S6A  shows distribution of maximum gene expression and the threshold in log 2 scale). This resulted in 16579 informative genes (distribution is presented in Supplementary Figure  S6B ).
Unlike nonnegative matrix factorization, both metagenes and weights can be positive or negative and ab initio the selection of the direction is random, depending on the initial estimation. Because of this, ICA may suffer from reduced reproducibility for at least some components. To mitigate this drawback, we ran the analysis multiple times (100 runs during the exploratory step and 1000 for the final analysis). The detailed algorithm of consensus ICA is described in the Supplementary Methods. Multithreading was implemented in R code to speed-up calculations using the foreach package and either doMC (Linux) or doSNOW (MS Windows) packages available in R/Bioconductor. The script of the implemented consensus ICA and other tools used here for investigation of the components are available online: https://gitlab.
The lists of top-contributing genes of each RIC were afterwards used for over-representation (enrichment) analysis. The 16579 informative genes were used as a background gene list and significantly enriched (adj. p-value<0. 01) GO terms were investigated. In order to simplify the interpretation and to increase the robustness for runs on different datasets, we reoriented the components in order to have the most significantly enriched categories associated with positive top-contributing genes (Supplementary Methods). For MICs, the direction could not be identified by enrichment analysis, therefore we reoriented only those MICs that showed strong negative correlation with RICs. As we ran ICA on a combination of reference and investigation datasets, several components captured the platform difference between these datasets.
Prediction of new sample classes. Random forest (RF) classifier, implemented in the randomForest R package  (Liaw & Wiener, 2002) , was used with the default settings to predict classes of patients. Columns of the weight matrix M were used as inputs and clinical variables (e. g. gender, sample type) as outputs. Each variable was analysed independently. First, leave-one-out cross-validation (LOOCV) was performed on the reference set in order to address the ability of predicting sample classes and estimate the accuracy of prediction. Then the RF, trained on all reference data, was used to predict classes for the new clinical samples of the investigation dataset. Integration of mRNA and miRNA expression data for survival prediction. Weights of the components (rows of matrix M) were statistically linked to patient survival using Cox partial hazard regression implemented in the survival package of R  (Therneau & Grambsch, 2000).
Our strategy to investigate the biological relevance of the components is presented in Supplementary Figure  S1. First, we attempted to connect the metagenes of all the components from the mRNA data to biological functions and cell types. We analysed separately the positively and negatively involved genes using several tools. Automatic analysis was done by topGO  (Alexa & Rahnenfuhrer, 2016)  followed by a manual analysis with Enrichr  (Kuleshov et al. , 2016 , http://amp. pharm. mssm. edu/Enrichr/) that checked for enrichment in multiple categories originated from various databases (we used Reactome 2016, GO Biological Processes 2017, Human.
Gene Atlas, ARCHS4 Tissues and Chromosome Location, http://amp. pharm. mssm. edu/Enrichr/#stats). In addition, we compared the metagenes to the ones previously published by Biton et al. (Biton et al. , 2014)  and assigned the component number to the reciprocally corresponding metagene as explained in Cantini et al. (Cantini, Kairov et al. , 2018)  using the DeconICA R package (https://zenodo. org/record/1250070). As enrichment of immune-related processes and functions was observed, we also correlated our metagenes to the immune cell type signature matrix named LM22  (Newman et al. , 2015)  in order to identify components originated from different types of leukocytes; cell-types were associated with components through highest correlation. Finally, for some components we confirmed their biological origin by correlating the metagenes with averaged gene expression profiles of cell types measured at a single-cell level and reported by Tirosh et al.
After automatically extracting miRNA-names and clusters from publication titles and abstracts by an inhouse Python script (available on demand), we compared those with the miRNA-metagenes comprised in the correlation-based networks, to identify a possible enrichment of miRNAs related to the proposed biological function. Additionally, we compared our miRNA-metagenes to miRCancer (http://mircancer. ecu. edu/  (Xie, Ding et al. , 2013) ). In order to explore the link or edges between MICs and RICs we extracted the target genes with a strong support from miRTarBase (http://mirtarbase. mbc. nctu. edu. tw/php/index. php,  (Chou et al. , 2018) ) for those miRNA-metagenes mapping the miRNAs and clusters found by literature mining. Additionally, we filtered these target genes to ensure that they were part of the reference gene set based on top-contributing genes as determined by the mRNA-reference set.
We then overlapped these target genes with the metagenes of the respective linked RIC and applied Enrichr tool through the automated Python-based API using the following reference gene set collections: KEGG_2016, GO_Biological_Process_2017b, GO_Cellular_Component_2017b, Jensen_TISSUES, Jensen_DISEASES (script available on demand). Moreover, we explored the overlapping target-and metagenes by STRING (https://string-db. org/,  (Szklarczyk et al. , 2017) ) to detect significantly enriched protein-protein interaction networks. Both the results of the link-analysis between RICs and MICs as well as the investigation of MIC-metagenes have been finalised by visual inspection and evaluation by biological experts. A detailed schematic illustration is shown in Supplementary Figure  S1. The involvement or the weight of each component in the new samples was not centred and scaled.
Tools: Consensus parallel ICA: https://gitlab. com/biomodlih/consica R/Bioconductor v. 3. 4. 3 with packages fastICA, doMC, doSNOW, topGO, randomForest, survival https://cran. r-project. org/ Enrichr http://amp. pharm. mssm. edu/Enrichr/ STRING https://string-db. org/ miRTarBase http://mirtarbase. mbc. nctu. edu. tw/php/index. php miRBase http://www. mirbase. org/ miRCancer http://mircancer. ecu. edu/ DeconICA https://github. com/sysbio-curie/DeconICA cyto_convert https://www. ncbi. nlm. nih. gov/genome/tools/cyto_convert/.
ACKNOWLEDGEMENT We would like to thank the patients for providing clinical material and the clinical staff of the Dermatology Unit at the University of Freiburg for professionally taking and handling all the samples. Bioinformatics analyses presented in this paper were partly carried out using the high-performance computing facilities of the University of Luxembourg (http://hpc. uni. lu). Furthermore, we thank Demetra Philippidou and Dr. Susanne Reinsbach for processing clinical samples and NHEM cells forming the investigation dataset. We acknowledge Dr. Aurélien Ginolhac and Cristina Maximo for critically reading the manuscript and for valuable input.