dataesr/scientific-paragraphs-categorization huggingface dataset:
Prediction of sample classes Random forest classifier, implemented in therandom-Forest`R-package  [24] , 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 crossvalidation (LOOCV) was performed on the discovery set in order to address the ability of predicting sample classes and estimate the accuracy of prediction. Then the random forest, trained on all discovery data, was used to predict classes for the new clinical samples of the investigation dataset. To ensure accuracy and robustness of our approach to select the number of components, we performed a nested cross-validation, excluding 20% of the data and using the remaining 80% to estimate the optimal number of components and then train the classifier (Additional file 1: Supplementary Methods).
Other dimensionality reduction methods In order to compare the performance of consensus ICA to other available tools, we run benchmarking of 7 approaches, applying them to mRNA expression data. First, we considered standard PCA of the joint dataset and PCA after correction for batch effects between discovery and investigation data sets using ComBat (packagesva`)  [25]  and XPN  [26]. Next, we applied the non-negative matrix factorization (NMF) implemented in theNMFp ackage  [7]  and low-rank approximation based multiomics data clustering (LRAcluster)  [27]. Finally, we investigated several non-linear dimensionality reduction methods, such as locally-linear embedding (LLE) implemented in thelle`package  [28] , Isomap (packageRDRToolbox`)  [29] , as well as t-SNE (packagetsne`)  [30]. To compare these methods to ICA, we performed 20 runs of 5-fold cross validation and estimated the accuracy of patient classification by random forest.
p-value< 0. 05) linked to survival and 0 for other. The applicability of the proposed score was checked using the independent validation dataset. This dataset was extracted from another study and was based on an independent cohort of the patients  [21]. In addition, gene expression for the validation cohort was measured using Affymetrix U133 Plus 2. 0 microarrays, while the discovery dataset was based on RNA-seq.
Biological relevance of the components Our strategy to investigate the biological relevance of the components is presented in Additional file 3: Fig. S1  (green box). 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 contributing genes using several tools. Automatic analysis was done by `topGO`R-package  [32]  followed by a manual analysis with Enrichr  [33]  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). In addition, we compared the metagenes to the ones previously published by Biton et al. [10]  and assigned the component number to the reciprocally corresponding metagene as explained in  [34]  using theDeconICA`R-package (https://zenodo.
ICA of combined data sets can remedy technical biases In this study, graphically outlined in Fig. 1  (see detailed schemes in Additional file 3: Figs. S1  and S2 ), we used public TCGA data as the discovery dataset, published microarray data  [21]  as a validation set. An investigation data set was based on newly obtained clinical samples described in Methods and Additional file 4: Table  S3. ICA was applied to two types of transcriptomic data: mRNA and miRNA expression. The number of components was chosen based on ability of ICA features to classify patients in the discovery set (see Additional file 1): 80 independent components were used for the deconvolution of mRNA data (named RIC1-80) and 40 for miRNA data (denoted as MIC1-40). ICA was run 1000 times in order to achieve robust results.
S4 A, B ). Therefore, the use of consensus ICA may be considered as optional for patient classification, but it is necessary for obtaining reproducible and biologically interpretable components and gene signatures. The combined discovery/investigation dataset profiled by RNA-seq is presented in the space of two first principal components (Fig. 2a ) and weights of two selected independent components (Fig. 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 discovery 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 as an example. The investigation data were clearly integrated within the discovery data and showed reasonable clustering in Fig.
S3B, D ) and only slightly lower for tumour subtypes (Additional file 3: Fig. S3F ), where averaged AUC was reported. The observations were confirmed by Wilcoxon test (p-values are reported in Additional file 3: Fig. S3 ) and by 2-factor ANOVA on log-transformed p-values. Post-hoc analysis confirmed that ICA, on average, gives features that are linked stronger to clinical groups, than PCA (Tukey's HSD p-value = 0. 0175). An even stronger correction effect of ICA was observed for miRNA data, for which discovery data were obtained Fig. 1  Visualization of the approach taken to data analysis. A large discovery dataset and a small investigation dataset from patients (both mRNA) were concatenated and analysed together by ICA. As a result, two matrices were obtained: S (metagenes), containing contribution of the genes to each component, and M (metasamples), presenting the weights of the components in the samples.
Fig. 2  Data overview in the space defined by principal and independent components. Data variability captured by the first components of PCA (a) and two selected components of ICA (b) in gene expression data. Independent components were selected based on the predictive power of their weights for patient gender (RIC3) and sample type (RIC5). MiRNA data showed even higher discrepancy comparing miRNA-seq and qPCR results by PCA (c). However, in the space of independent components (MIC1 and MIC9), the samples studied by miRNA-seq and qPCR overlap (d) Next to cross-validation tests that were run directly on the discovery data, we applied ICA and random forest classification on two independent datasets joint with the discovery data: public validation (E-GEOD-19234, 44 samples) and in-house clinical investigation data (5 samples).
Although these 11 RICs and 3 MICs were statistically linked to survival in our discovery set, the predictive power of any of them was not sufficient to predict survival of new patients. Therefore, we combined the weights of these components into a risk score (RS) as described in Methods. Combined RS showed high significance (p-value = 2. 2e-13) for the TCGA dataset. In order to validate the proposed risk scoring approach on an independent cohort of patients, we applied it on the validation dataset. The components that showed a significant link to survival (adj. p-value< 0. 05) on the discovery set were then used to compose RS for the validation data and also showed significant prognostic properties (LHR = 0. 87, p-value = 0. 0013); Kaplan-Meier plots are shown in Fig. 4. The developed RS separated patients with low hazard (only one death among 7 patients, blue line in the validation cohort, Fig.
4b ) from the group of patients with a high risk score. For the three primary melanoma samples from the investigation set, the calculated RS was the highest for P6PM (RS = 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 risk score, suitable to predict patient survival.
5  f-g). RIC27 showed a very similar collection of enriched gene sets, but was much less correlated to this particular score, suggesting that ICA shows better sensitivity and captures more cell subtypes than ESTIMATE. 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. pvalue = 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  [10]  (Additional file 3: Fig. S5A ). Neither RIC28 nor RIC37 were detected by ESTIMATE scoring. Components RIC74 and RIC79 contained a very limited number of top-contributing genes, but both were significantly linked to survival (p-values of 1.
2 (adj. p-value = 7. 3e-3). The samples from the investigation cohort were characterised by the involvement of the above immune responserelated components (Fig. 5d ). The results are presented in Fig. 6. All components linked to subpopulations of immune cells (RIC2, RIC25, RIC57, MIC20) showed little involvement in the patients of investigation cohort suggesting low overall immune reactions to the tumour except specific interferon responses, which had high weights in the investigation samples (RIC28, RIC37). Similarly, we checked behaviour of these components for the validation dataset (Additional file 3: Fig. S7 ). RIC2 and RIC25 showed tendency to predict better survival (LHR < 0) and their weights are higher for censored patients than for dead. However the p-values from Cox regression on 44 validation samples were not conclusive (RIC2: 0.
Skin-related components RIC5, RIC7, RIC19, RIC31 all showed an enrichment in GO terms related to processes of the skin including epidermis development (adj. 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  [18]. MIC9 with the skin-specific miR-203  [40] , 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-values <1e-33). This Fig. 6  Biologically relevant components and their ranked weights in the investigation dataset. Rank for samples is calculated in comparison to the TCGA discovery set (redweight above median in TCGA samples, bluebelow) component was as well associated with patient survival (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 bimodal distribution in the discovery dataset (see two clusters in Fig. 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 epithelialmesenchymal transition (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.
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 SKCM. These included (a) a large cohort of SKCM patients from TCGA used as discovery 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 discovery 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 discovery set (Fig.
1 ). This was particularly evident for miRNA data where the discovery 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. Such batch correction, of course, could also be performed by other methods. We tested several correction methods together with standard dimensionality reduction methods and showed that overall, ICA performed best across them. We recently applied ICA-based batch correction on single-cell RNAseq data and confirmed its usefulness  [44]. The fact that ICA should be re-run for every series of new samples could be considered as a drawback of our approach.
However, similarly to PCA, recalculation of the components does not require supervision and could be done automatically. In the case when investigation and discovery datasets come from the same distribution, one can use the matrix S obtained from the discovery dataset in order to define the weights (M) for the samples forming the investigation dataset (1). However, in reality, the variability in the data requires recalculation of the components for the new investigated samples. We demonstrate here that the weights of independent components can be used as predictive features of patient subgroups and can be linked to patient survival. We also propose a method to select the number of components, based on the required classification task (Additional file 2 and Additional file 3: Fig. S10 ). While the ICA-based feature extraction method has been previously discussed (e.
Conclusions In conclusion, we used the consensus ICA method to combine transcriptomics data of melanoma patients with large public datasets. Here we showed successfully that the ICA-based decomposition separates true biologically relevant transcriptional signals from technical biases. The obtained ICA-based features were used to predict cancer subtypes and patient survival. We also showed how cellular composition and biological signals can be uncovered within new clinical samples. Transcriptional signals from immune cells, melanocytes, keratinocytes and stromal cells were identified and confirmed by comparison to published signatures. We demonstrated that some of the identified signals, including immune activity and cell proliferation, are linked to the aggressiveness of tumours and could influence patient survival.
Acknowledgements 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. The results shown here are in part based upon data generated by the  TCGA Research Network : https:// www. cancer. gov/tcga.
Availability of data and materials The developed tools are available at https://gitlab. com/biomodlih/consica. RNA-seq data for the samples used to build the investigation dataset are available by GEO accession number GSE116111 (https://www. ncbi. nlm. nih. gov/geo/). Ct-values for all quantified miRNAs are available in Additional file 4: Table  S4. The results of ICA and their biological interpretation are presented in Additional file 2. The validation gene expression dataset used in this study is available from ArrayExpress under the accession number E-GEOD-19234.
This is especially important for several unique exoplanets, like the above-mentioned WASP-12 or WASP-4 demonstrating possible TTV trends. In particular, relatively little attention was paid so far to a yet another explanation of such trends, based on the light-travel effect induced by outer bodies  (Irwin 1952). Finally, this work represents the first big practical test of the EXPANSION project (EXoPlanetary trANsit Search with an International Observational Network), grown on the basis of the ETD (Exoplanet Transit Database) that was used by  Baluev et al. (2015). Now EXPANSION is a standalone international project joining a network of several dozens of relatively small-aperture telescopes, aimed to monitor the exoplanetary transits  (Sokov et al. 2018). This network covers amateur as well as professional observatories spreaded over the world in the both hemispheres.
PHOTOMETRIC AND DOPPLER DATA The EXPANSION project performs a long-term monitoring of exoplanetary transits. Amateur and professional observatories from Russia, Europe, North and South Americas with relatively small telescopes from 25 cm to 2 m are used in the photometric observations  (Sokov et al. 2018). We used data from this network, including all the data from ETD that were used in  (Baluev et al. 2015). Additionally, we used lightcurves published in the literature or kindly provided by the observers, as listed in Table  1. Most of them are available in the VIZIER database. We expanded our targets list by seven exoplanets: Qatar-2, WASP-3, -6, -12, HAT-P-3, -13, and XO-5, thus increasing their number to 17. The total amount of the input data has grown considerably. This time we had ∼ 3×10 5 photometric measurements in ∼ 1000 lightcurves, compared to ∼ 8 × 10 4 measurements in ∼ 300 lightcurves processed by  Baluev et al.
(2015). Whenever necessary, the timestamps in the photometric series were transformed to the BJDTDB system by means of the public IDL software developed by  Eastman et al. (2010). To perform this reduction, we used ICRS coordinates through the SIMBAD database which originate from GAIA DR2  (Brown et al. 2018). We did not apply any correction to these coordinates due to proper motion, since this would imply only a negligible correction to the time (below ∼ 0. 1 sec). Additionally, we used the precision radial velocity (RV) measurements obtained from the archival spectra of the HARPS, HARPS-N, SOPHIE, and HIRES spectrographs. This involves the following targets from our photometry sample: Corot-2, GJ 436, TrES-1, WASP-2, -4, -5, -6, -12, HD 189733, XO-2N. The spectra were processed with the HARPS-TERRA pipeline (Anglada-Escudé & Butler.
2012 ). Some of these data represent reprocessed versions of the RV data available in the literature, e. g. from  (Baluev et al. 2015) , and some are new. Whenever performing a self-consistent transit and radial velocity analysis we transform all the Doppler time stamps to the BJDTDB system consistent with the photometry. However, the RV data that we release here correspond to the UTC rather than TDB system (as traditionally adopted for this type of the data). Since  Wilson et al. (2008) , additional 31 RV measurements have been obtained for WASP-4 with the high resolution spectrograph CORALIE on the Swiss 1. 2 m Euler telescope at La Silla Observatory, Chilie  (Queloz et al. 2001). RVs were re-computed for the new data and the dataset presented in  Wilson et al. (2008) , for 45 measurements in total, by cross-correlating each spectrum with a G2 binary mask, using the standard CORALIE data-reduction pipeline.
By adding two more probe fits starting from τ closer to the low and upper limits of the range, we could robustly detect the red noise in several lightcurves additionally. But even with these very mild detection criteria and multiple trial fits it appeared that only 1/4 to 1/3 of our lightcurves (depending on the target) revealed an individually fittable red noise. This is in agreement with  Baluev et al. (2015) , however such a low fraction of the rednoised lightcurves still appears surprising. The red noise may exist in the rest of lightcurves too, but with ill-fitted individual parameters. Therefore, leaving the noise models in such a partial model-mixed state might make the resulting TTV data less homogeneous. For example, the uncertainties in the white-noise portion of TTV data may appear systematically smaller than in the red-noise one.
Therefore, we need some benchmark of the accuracy and quality of our TTV data. This can be done by comparing them with analogous TTV data from other published works. However, most of the published TTV data were derived by different teams who used different techniques and different assumptions (e. g. concerning the limb-darkening). Hence, their mixed compilations cannot usually serve as reliable benchmarks. We need a long record of TTV data, obtained mostly by the same team. In our target list, only WASP-12 perfectly suits our needs. It has ∼ 200 observed transit lightcurves in total, and about half of them were processed by the same group  (Maciejewski et al. 2013 (Maciejewski et al. , 2016 (Maciejewski et al. , 2018b)). Simultaneously, these data were obtained at quite different telescopes located in different astroclimate conditions.
Therefore, they have different quality characteristics, offering the necessary degree of physical heterogeneity. Here we used the transit times from  (Maciejewski et al. 2018b ) and  (Maciejewski et al. 2016 ) that included the most reprocessed lightcurves of  (Maciejewski et al. 2013 ). We did not include the timings from the 2013 paper not reprocessed in the 2016 one (since they would be statistically different). After that we sampled the same transits from our homogeneous data release. Thus we obtained two similar TTV time series to be compared, each containing 73 data points at the same epochs. We fitted both them with a quadratic trend model (  5 ), resulting in almost identical trend fits. We then computed the resulting RMS: 22. 1 sec for the Maciejewski et al. data and 24. 9 sec for our data release.
The Maciejewski et al. data winning with a slightly smaller scatter of the residuals, though this 12 per cent difference is comparable to the probable statistical uncertainty (1/ √ N for N = 73 yields the same 12 per cent). Therefore, the intrinsic statistical accuracy of our processing pipeline appears similar high quality TTV data sets available in the literature. However, our data reveal important difference in another aspect. The value of the reduced χ 2 for the quadratic TTV model is 0. 79 for Maciejewski et al. , implying that they tend to overestimate their TTV uncertainties by 12 per cent on average. On the contrary, our data imply the reduced χ 2 of 1. 49, which means that our TTV uncertainties appear underestimated by the factor √ 1. 49, or by ∼ 22 per cent on average. We notice that it is quite frequent that the uncertainties reported for some measurements have a remarkable systematic bias.