Manually curated and harmonised transcriptomics datasets of psoriasis and atopic dermatitis patients

We present manually curated transcriptomics data of psoriasis and atopic dermatitis patients retrieved from the NCBI Gene Expression Omnibus and EBI ArrayExpress repositories. We collected 39 transcriptomics datasets, deriving from DNA microarrays and RNA-Sequencing technologies, for a total of 1677 samples. We provide quality-checked, homogenised and preprocessed gene expression matrices and their corresponding metadata tables along with the estimated surrogate variables. These data represent a ready-made valuable source of knowledge for translational researchers in the dermatology field.

We present manually curated transcriptomics data of psoriasis and atopic dermatitis patients retrieved from the NCBI Gene Expression Omnibus and EBI ArrayExpress repositories. We collected 39 transcriptomics datasets, deriving from DNA microarrays and RNA-Sequencing technologies, for a total of 1677 samples. We provide quality-checked, homogenised and preprocessed gene expression matrices and their corresponding metadata tables along with the estimated surrogate variables. These data represent a ready-made valuable source of knowledge for translational researchers in the dermatology field.

Background & Summary
Psoriasis (PSO) and Atopic dermatitis (AD) are among the most common inflammatory skin disorders associated with immunologic impairment. While the first signs of AD tend to appear in the early childhood, the manifestation of PSO is most common during the third decade of life 1 . Both the diseases have a substantial negative impact on the quality of life of affected patients. Although a number of therapeutic approaches have been developed in the last two decades to mitigate PSO and AD symptoms, their pathophysiology is still not completely understood 2,3 . AD is believed to be driven by epidermal barrier disruption, activation of specific T-cell subsets, and dysbiosis of the commensal skin microbiome 2 while psoriatic inflammation is sustained by uncontrolled responses of the innate and adaptive cutaneous immune system, which lead to intense keratinocyte proliferation and dysfunctional differentiation 4 .
Transcriptomics technologies, such as DNA microarray and RNA Sequencing (RNA-Seq), have been used to characterise the molecular alterations of human diseases 5 , including PSO and AD. To date, only marginal efforts have been carried out in order to collect, quality-check and harmonize PSO-and AD-related transcriptomics data in order to make them easily reusable by the research community. Therefore, the motivation behind this study was to create a source of ready-to-use data of gene expression profiles of PSO and AD patients derived from both DNA microarray and RNA-Seq publicly available datasets.
The preprocessed and harmonized microarray data provided in this study were collected from the NCBI Gene Expression Omnibus (GEO) and EBI ArrayExpress public repositories, while the RNA-Seq datasets were retrieved from the European Nucleotide Archive (ENA). Overall, 26 microarrays datasets were collected, for a total of 991 samples, 632 of which from patients affected by psoriasis and 70 by atopic dermatitis. Some of the microarray datasets contain samples collected from patients affected by other skin diseases such as psoriatic arthritis, psoriasis sebaceous hyperplasia, palmoplantar pustulosis, lichen planus and discoid lupus. These datasets were generated with commercially available Affymetrix and Agilent platforms. All of the analytical steps performed in this work were carried out through the use of the eUtopia software 6 . We also retrieved 13 RNA-Seq datasets, for a total of 686 samples, 392 of which from patients affected by psoriasis and 94 by atopic dermatitis. RNA-seq data were mostly produced through Illumina platforms, while a minority of datasets were produced through other platforms. All the datasets underwent meta-data curation and harmonisation, data quality check and preprocessing with standardised procedures. The curation and harmonisation of the meta-data consisted in the definition and usage of a common data model for all of the collected datasets. The data models, to which the raw meta-data were mapped to, are reported in the data dictionary files (enclosed with the preprocessed data). The data dictionary describes all the variables reported in the final metadata tables. For each variable, the description, type and allowed values are reported. At the same time, this work is aimed at homogenising the preprocessing procedures in order to improve the comparability of the gene expression data across different studies and platforms. Therefore, in this work we provide meta-data tables, along with the inferred surrogate batch variables, as well as the preprocessed gene expression estimates.
Our analysis significantly increases the FAIRness 7 of publicly available PSO and AD transcriptomics data and represents a valuable "ready-to-use" resource available to the scientific community.

Methods
Microarray data. Data collection and homogenization. Transcriptomics data generated by DNA microarrays of psoriasis and atopic dermatitis patients were retrieved from NCBI GEO 8 (GEO -https://www.ncbi.nlm. nih.gov/geo/) and EBI ArrayExpress (https://www.ebi.ac.uk/arrayexpress/) repositories by using the R packages GEOquery 9 and ArrayExpress 10,11 , respectively. For each dataset, a table specifying the disease (psoriasis/atopic dermatitis) and the origin of biopsy (lesional/non-lesional sample) in addition to other phenotypic information was also retrieved. Since the phenotypic information was heterogeneous across the datasets, rigorous harmonization procedure was performed. The GEO and Array Express identifiers of the retrieved datasets are reported in Tables 1-3.
Data quality check. The retrieved datasets were thoroughly quality checked. In particular, each sample was evaluated by visual inspection of the array pseudo-images, quality check reports and density plots of probe intensities by using the eUTOPIA software 6 . Further, outlier detection step, based on the sample distributions, was performed within each dataset by using ad hoc R scripts (see Code Availability section).
The distributions of the values of these three metrics were investigated by means of boxplots and the sample outlierness was evaluated for each measure based on the data distribution. Eventually, a concordance outlierness score was computed across the three metrics. In particular, a sample was removed from the analysis if considered an outlier in at least two out of three metrics, one of them being the RNA degradation curve.
Normalization. Data normalization was performed by using the eUTOPIA software. Affymetrixbased studies were normalized by using the justRMA from the R affy v1.66.0 package 14 . Agilent-based studies were quantile normalized with the normalizeQuantiles function from the limma v3.44.3 package 15 .
Surrogate variable analysis. In order to investigate the effect of unknown batches that might mask biological variability, Surrogate Variable Analysis (SVA) was performed with the eUtopia software, which implements the sva R package 16 . The analysis was performed by using origin of biopsy or diagnosis as variable of interest. The other biological variables (if present and if not confounded with the variable of interest) were used as covariates 6 . The estimated surrogate variables for each dataset are included in the meta-data tables.
Probe annotation. Custom annotation files (CDF files) were downloaded from Brainarray (http://brainarray. mbni.med.umich.edu/Brainarray/Database/CustomCDF/CDF\_download.asp) for Affymetrix-based microarrays. The latest version of Agilent probe annotation was retrieved from the Agilent website (https: //earray.chem. agilent.com/earray/). The probesets were mapped to the Ensembl gene IDs and the expression matrix was aggregated by computing the median of the expression of the Agilent probes mapping to the same Ensembl transcript ID. The entire DNA microarray data preprocessing is depicted in Fig. 1. RNA Sequencing data. Data collection and homogenization. Raw files in ".fastq" format were retrieved from the European Nucleotide Archive (ENA). Along with the raw data files, the metadata tables reporting the samplewise clinical features for each dataset were also collected. As for the DNA microarray data, the metadata tables of RNASeq data were carefully harmonized to improve the across-datasets comparability. Phenotypic www.nature.com/scientificdata www.nature.com/scientificdata/ information for each dataset is reported along with the gene expression tables. GEO and ENA identifiers of the retrieved datasets are reported in Tables 1-3.
Quality control. All the RNA-Seq datasets underwent quality check through the use of FastQC v0.11.7 (https:// www.bioinformatics.babraham.ac.uk/projects/fastq c/). Reads were trimmed for low-quality ends in addition to adapters removal by TrimGalore v0.4.4_dev (http://www.bioinformatics.babraham.ac.uk/ projects/trim_galore/). In particular, the reads were trimmed if the Phred score was lower than 20 and discarded if the number of undetected nucleotides was greater than 50. The trimmed and adapter-clipped raw files were further quality checked with FastQC v0.11.7.
Read alignment. RNA Sequencing reads were then aligned against the human reference genome assembly GRCh38. The alignment was performed through the use of the HISAT2 algorithm 17   Low counts filtering. In order to filter out the transcripts with low expression levels in all the samples of each dataset, the proportion test strategy was used as implemented in the function filtered.data of the R package NOISeq v2.31.0 21 .
Normalization. RNASeq expression data were normalized using the upper quantile method 22 implemented in the R package NOISeq v2.31.0.
Surrogate Variable Analysis. As for the DNA microarray data, in order to identify unknown sources of technical variability, a Surrogate Variable Analysis (SVA) was performed through the use of the svaseq function implemented in the sva v3.36.0 R/Bioconductor package 16 . The analysis was performed by using disease state or diagnosis as variable of interest. The other biological variables (if present and if not confounded with the variable of interest) were used as covariates 6 . The estimated surrogate variables for each dataset are included in the meta-data tables, along with the gene expression tables. The entire RNA-Seq data preprocessing is depicted in Fig. 2.

Data Records
The complete list of DNA microarray and RNA Sequencing datasets discussed in this work is reported in Tables 1-3. All the preprocessed transcriptomics data, along with harmonised meta-data, were submitted to Zenodo 23 .

Technical Validation
DNA microarray and RNA-Seq data are linked to clinical meta-data, reporting multiple information such as gender, age or the treatment (including e.g. drug dose). Additionally, sample meta-data is recorded, such as the tissue type a sample was taken from, or whether the tissue derives from a lesional or nonlesional sample.
In order to ensure that the data is recorded in a consistent and well-formed way, we created data dictionaries describing each of these variables. The data dictionaries contain detailed information describing the content of a variable, the data type (numeric, categorical, text, date, etc), the allowed values of categorical data or ranges of numeric variables.
The data was validated by checking compliance with the rules encoded in the data dictionaries. Data that was found not to comply with the rules was manually curated by consulting the original data sources. In fact, a large proportion of the datasets were found not to meet the requirements encoded in the data dictionaries. For www.nature.com/scientificdata www.nature.com/scientificdata/ instance, big heterogeneity was found in the description of the skin status. "Involved skin", "psoriatic skin" were reported in order to describe the "lesional" status of the skin. "Normal", "ctrl", "Non-involved skin of healthy individual" were used to describe the "healthy control" samples. Yet, to define the gender, "m", "f ", "male" and "female" were used across the datasets. All of these variables were mapped to the allowed values reported in the data dictionaries to improve the comparability across the datasets.

Usage Notes
The transcriptomics data presented in this article is an unprecedented source of preprocessed, harmonized, "ready-to-use" and FAIR datasets, made available to the scientific community. Data derived from both DNA microarray and RNASeq technologies can be exploited in order to uncover the molecular mechanisms underlying psoriasis and atopic dermatitis. Differential expression analysis can be carried for instance by the limma package 15 for the microarray data, and the edgeR 24 , DESeq 2 25 or NOISeq 21 packages for RNA-Seq data, respectively. Functional analysis of differentially expressed genes can be performed by using FunMappOne 26 , the R/ Bioconductor package ReactomePA 27 or Ingenuity Pathway Analysis (Qiagen, http://www.ingenuity.com/products/ipa). The inference and analysis of co-expression networks can be performed, for instance, by using the INfORM tool 28 . Altogether, these analyses can aid the stratification of PSO and AD patients, the identification of relevant biomarkers and novel therapeutic targets.

Code availability
R scripts for the analysis of DNA microarray and RNA-Seq transcriptomics data are available for download at: https://github.com/Greco-Lab/psoriasis-dermatitis-analysis.