\chapter{Methods}

In this section, I will review the methodology that was used in the two sections of this doctoral thesis. First, I will explain the components and structure of the research framework that was developed to organize and analyze efficiently the data of the Johns Hopkins Myositis Center longitudinal cohort study. Second, I will detail the most relevant statistical techniques that were used to perform the series of epidemiologic analysis of this clinical cohort. Finally, I will summarize the design of the Myositis muscle RNA sequencing studies, focusing on the bioinformatic analysis of the resulting sequencing data that was obtained. The detailed information about the methods of each of the individual studies that compose this doctoral thesis is included in \autoref{chap:appendix}.

\section{Longitudinal cohort study framework}

\subsection{The Johns Hopkins Myositis Center longitudinal cohort study}
In 2007 the Johns Hopkins Hospital founded a monographic center to treat and do research in patients with different types of myositis. It is a multidisciplinary unit combining the expertise of Rheumatologists, Neurologists, Pulmonologists, and Physical Therapy specialists. Over the years, more than two thousand patients suspected of having an inflammatory myositis have been evaluated and followed over time, becoming the largest longitudinal myositis cohort in the world.

\subsection{Database design}

In 2013, a research database was created to organize the information that was available for the patients of the myositis cohort. The first challenge was to collect all the information, that was originally spread over different files with data of variable quality, included inconsistencies, duplicated observations, was not normalized, and relied heavily on text descriptions.

To organize this heterogeneous collection of data, first, we had to choose what would be the main blocks of information that were relevant from a myositis research standpoint and how to organize them in a way that will be useful not just for a single project but for many studies to come.

It was decided that the minimum set of data that we needed included a unique identifier for each patient and epidemiological data, like the gender, date of birth, or the race. Moreover, it was necessary to include, among others: the general phenotype of the patients; the family history; the medications; the clinical features as they were longitudinally recorded; the results of the laboratory tests, \gls{mri}, muscle biopsy, pulmonary function tests, and chest \gls{ct}. Finally, we would need to include the information of each one of the autoantibodies, both the consensus interpretation and all the individual serologic results.

Also, we had to select the best technical solution to store the data according to the expected level of usage, limited time to maintain it, and availability of the software in the terminals of the different users. It was decided that using Microsoft Access fulfilled the required criteria for its initial intended usage (Figure \ref{fig:db_frontend}). It could be used concurrently by up to 20 users, both the front-end and the back-end could be easily maintained by a single researcher, and it had institutional support in all the user terminals. The database was designed so that only a limited number of people would have access to the back-end, but all the researchers in the protocol could view, query, and modify the data in the front-end. Over the years this framework stopped being compliant with the institutional Hopkins policies and the back-end was migrated to a Microsoft SQL Server, with ongoing plans for the front-end to be moved to a web framework.

A central database with un-duplicated information was shared by all the concurrent research projects, but each researcher was able to work in different subsets of patient and specific sets of variables by defining filters on each project's front-end. By having a central data repository, different projects would review and enter new data for themselves but, by doing so, they would be helping to curate the data for other concurrent and future projects, increasing the overall productivity of the group.

\begin{figure}
	\centering
	\includegraphics[width=1\linewidth]{db_frontend}
	\caption{Appearance of the front-end of the myositis database.}
	\label{fig:db_frontend}
\end{figure}

\subsection{Data mining}
Once the database was built, it was obvious that importing certain types of data could be automated. Originally, it was not feasible to access the back-end of the Hopkins electronic patient record system directly, but users had access to the front-end of the application to collect the data of those patients consented for the study. Also, there was no application programming interface to automate the data mining process. Thus, the only option, if we wanted to automate importing the thousands of registries, would be to do it through mouse and keyboard automation. A first crude data import was successful at automating the extraction of all the relevant clinical and laboratory data (e.g. creatine kinase, aldolase, AST, ALT, and autoantibody reports). This approach was further perfected through a script with modules to include each one of the sequences of key-strokes and mouse-clicks necessary to import each section of the data for every patient. Finally, contingency systems were put in place to reset and tag abnormal data imports when the system did not get the correct set of information.

A portion of the data required expert interpretation from researchers trained to evaluate medical jargon and familiarized with the type of disease. The front-end of the database was progressively modified to simplify these data-entry tasks and include quality control variables that could be used to detect errors during the process.

\subsection{Data wrangling}
Once the information was available in the database, a key part of any data analysis effort was to clean the raw data to adapt it to the specific analyses. It was key to automate these steps in order to be able to update the results if new data was added, be able to design the analysis at the same time that the data was being collected, and reuse the data cleaning steps to optimize productivity. Most analyses of clinical data were performed using Stata and thus, Stata scripts used the raw Access tables as the starting point of the process, which facilitated updating the data by replacing the raw tables. It was first necessary to integrate all the different tables in a single comprehensive dataset while allowing for enough flexibility so the process of integrating the tables was modular. A set of functions was built to merge each auxiliary table to the main one including the epidemiological features of the patients. Also, the variables were labeled and calculated fields were generated. Finally, a limited set of “working datasets” were built for the different families of analysis in order to avoid delays each time it was necessary to modify a piece of the code if the data to be used had not changed.

\subsection{Data analysis}
\subsubsection{Univariate analysis}
Univariate analysis is often the first step in an exploratory analysis. Hypothesis contrast using Fisher's exact test, Chi-squared, Student's T-tests, and Wilcoxon's rank-sum test is the bread and butter of any epidemiological study. However, performing these tests is often extremely labor-intensive due to the massive amount of comparisons that may be needed when multiple variables are being explored simultaneously. Also, generating publication-quality tables is often a source of productivity loss in research. For this reason, a Stata program called \textit{table\textunderscore compare} was built to simplify performing the hypothesis testing, deciding automatically when to perform Fisher's or Chi-square tests and manually to indicate if continuous variables should be considered normally-distributed or nonparametric tests should be applied. Moreover, this tool allows performing multivariate analysis, handles paired data, and can export the results ready for publication. This tool allowed us to speed up performing and updating the different analyses.

\subsubsection{Multilevel regression models}
A complication of observational longitudinal cohort studies, when performed in a clinical practice setting, is that each patient has a different number of observations (e.g. visits or laboratory determinations). In order to avoid bias, it is necessary to use techniques to take into account this uneven length of follow-up. A simplistic solution would be averaging the periods of observation and adjusting for the time of follow-up, but this would result in an unnecessary loss in statistical power. A valid alternative in this situation is to use a family of regression methods called multilevel regression models or random effects mixed models. These models include the fixed effects of conventional regression analysis but then allow for random components grouped in categories (e.g. each individual patient will have a different evolution of the \gls{ck} levels over time). Multilevel regression models allow for random intercepts (e.g. each patient having higher or lower overall \gls{ck}) and random slopes (e.g. each patient can have a faster or slower decrease in the \gls{ck} levels). Strictly speaking, the full model with random intercepts and random slopes is just necessary if it is significantly different from the model with random intercepts, and the model with random intercepts only if it is significantly different from the standard regression model. However, doing these model evaluations adds a layer of complexity that may be eliminated by using the most complex model with random slopes and random intercepts that will control for all possible variation caused by the heterogeneity of the patient population and the differential follow-up.

\subsubsection{Factor analysis of mixed data}
To model the phenotype of the patients we used factor analysis of mixed data. The phenotype is a latent variable, i.e., a variable that cannot be directly observed but has to be inferred through a mathematical model.

Factor analysis of mixed data is a principal component method dedicated to exploring data containing both continuous and categorical variables. The continuous variables are scaled to unit variance and the categorical variables are transformed into a disjunctive data table and then scaled using the specific scaling of multiple correspondence analysis. This ensures to balance the influence of both continuous and categorical variables in the analysis. It means that both variables are on an equal foot to determine the dimensions of variability. We used the package FactoMinerR v.2.1 to perform the factor analysis of mixed data and factoextra v. 1.0.6 to obtain the scree plots and the variable weight plots.

As input for the factor analysis of mixed data we selected a set of clinical, epidemiological, and laboratory parameters that were: 1) well documented in the literature to be associated with the phenotype of patients with myositis, 2) systematically collected in our cohort, and 3) well defined. Thus, we included:

\begin{enumerate}
	\item Epidemiologic variables: Gender, race, age at onset.
	
	\item Clinical variables: presence or absence during the course of the disease of muscle weakness, interstitial lung disease, arthritis, heliotrope or Gottron’s rash, calcinosis, Raynaud’s phenomenon, mechanic’s hands, dysphagia, and fevers.
	
	\item Laboratory values: Maximum \gls{ck}, presence of anti-Ro52 autoantibodies.
\end{enumerate}

Anti-Ro52 autoantibodies were included because they are associated with the severity of the disease and with specific clinical features in patients with myositis.\cite{PinalFernandez2017a,Vancsa2009,Marie2012,Bauhammer2016}

The detailed distribution of the muscle weakness, biopsy features, \gls{mri} patterns, or \gls{emg} findings were excluded from this analysis because they were not available for all patients, and restricting the sample size could bias the study.

\subsubsection{Model comparison}
Both the \gls{aic} and the \gls{bic} are estimators of out-of-sample prediction error and thereby the relative quality of statistical models for a given set of data.\cite{Akaike1974,Schwarz1978} Given a collection of models for the data, \gls{aic} and \gls{bic} estimate the quality of each model, relative to each of the other models. Thus, they provide a means for model selection.

The formula for \gls{aic} is:
\[AIC =  2-2\ln(\hat{L})\]

Whereas \gls{bic} is formally defined as:
\[BIC =  k\ln(n)-2\ln(\hat{L})\]

Where:
\begin{description}
	\item $\hat{L}$=the maximized value of the likelihood function of the model
	
	\item $n$=the number of observations
	
	\item $k$=the number of parameters estimated by the model
\end{description}

\gls{aic} and \gls{bic} are founded on information theory. When a statistical model is used to represent the process that generated the data, the representation will rarely be exact; so some information will be lost by using the model to represent the process. \gls{aic} and \gls{bic} estimate the relative amount of information lost by a given model: the less information a model loses, the higher the quality of that model.

In estimating the amount of information lost by a model, \gls{aic} and \gls{bic} deal with the trade-off between the goodness of fit of the model and the simplicity of the model. In other words, \gls{aic} and \gls{bic} deal with both the risk of overfitting and the risk of underfitting.

Given a set of candidate models for the data, the preferred model is the one with the minimum \gls{aic} or \gls{bic} value. Thus, \gls{aic} and \gls{bic} reward goodness of fit (as assessed by the likelihood function), but it also includes a penalty that is an increasing function of the number of estimated parameters. The penalty discourages overfitting, which is desired because increasing the number of parameters in the model almost always improves the goodness of the fit.

The main difference between \gls{aic} and \gls{bic} is that the penalty for the number of parameters is larger for \gls{bic}. Thus, if ‘k ‘is the number of parameters and ‘n’ is the number of observations, the penalty for \gls{aic} is $2k$, whereas for \gls{bic} is $\ln(n)k$.

Various authors have suggested guidelines to interpret the magnitude of the differences in \gls{aic} and \gls{bic} between two models ($\Delta$AIC and $\Delta$BIC respectively). Thus, Burnham and Anderson\cite{Burnham2004} suggested that models having $\Delta$AIC $\leq$ 2 have substantial support (evidence) to believe that are equivalent, those in which 4 $\leq$ $\Delta$AIC $\leq$ 7 have considerably less support, and models having $\Delta$AIC > 10 have essentially no support. Alternatively, Raftery\cite{Raftery1995} suggested that the evidence that two models are not equivalent would be weak with a $\Delta$BIC between 0-2, positive with $\Delta$BIC between 2-6, strong with $\Delta$BIC between 6-10, and very strong with $\Delta$BIC over 10. Based on this we selected a threshold of $\Delta$AIC and $\Delta$BIC of 10 to consider one model superior to others.

\subsubsection{Graphical analysis of longitudinal data}
A key part of any scientific study is to represent the results visually in a way that is easy and quick to understand even for non-experts. To do so in an efficient manner, functions were build to make the most common graphical analysis, including individual and aggregated patient evolution plots and Kaplan-Meier curves. A useful graphical method to represent the evolution of nonlinear parameters longitudinally is the \gls{lowess}. \gls{lowess} fits simple models to localized segments of the data to build up a function that describes the deterministic part of the variation in the data, point by point.
The advantage of \gls{lowess} over other methods (e.g. quadratic regression) is that it does not need to specify a function to fit the model to the data, making it simple and flexible for complex graphical representations. Alternatively, it requires a dense cloud of observations to be stable and it can be easily biased by outliers if the local density of data is low. Also, \gls{lowess} does not return a simple mathematical function and, thus, is complicated to use for predictive purposes. Finally, \gls{lowess} is relatively computationally intensive, but with the range of observations that we used (in the thousands), this was a negligible issue.

\subsubsection{Survival analysis}
Many clinical questions in our dataset were related to the rate of development of the clinical manifestations from the onset of the disease. Also, death is, in all severe systemic diseases, an important prognostic predictor to analyze. Both the timing to develop a clinical manifestation and the time to death were appropriate questions to analyze using COX regression and Kaplan-Meier curves.

\subsubsection{Standardized cancer and mortality rates}
A recurrent question in cohort studies of patients with myositis is if certain groups of patients have higher cancer or mortality rates. Internal comparisons can be helpful by defining if some groups have higher or lower mortality rates than others or than the rest of the individuals in the cohort. However, the most relevant comparison is often with the general population through indirect standardization. A limitation of these types of comparisons is that they require well-annotated epidemiological data. In the US, both survival and cancer data can be acquired through the Centers for Disease and Control Prevention. Specifically, mortality data can be obtained from the United States Cancer Statistics (UCSC) registry (\url{www.cdc.gov/cancer/uscs/}) while the survival data can be acquired from the Compressed Mortality File (\url{wonder.cdc.gov/wonder/help/cmf.html}).

\section{Myositis muscle RNA sequencing studies}
Obtaining insight into the pathogenesis of the disease is the first basic step to understand the variability in the clinical manifestations, tailor specific treatments for the disease, and estimate how many disease categories there are. In order to do so, Next Generation Sequencing techniques are extremely powerful since they allow us to get great volumes of data on the transcriptional changes of the affected tissues. For this section of the study, we used RNA sequencing from myositis muscle biopsies to understand the mechanisms of muscle damage in these autoimmune diseases. Moreover, for certain studies, we also used RNA sequencing data from differentiated human skeletal muscle cells and injured mouse muscle.

\subsection{Cultured human skeletal muscle cells}
Human skeletal muscle myoblasts (Lonza) were cultured according to the manufacturer’s protocol. When 80\% confluent, the cultures were induced to differentiate into myotubes by replacing the growth medium with differentiation medium (DMEM, 2\% horse serum, and L-glutamine).  Two plates of cells were collected for RNA extraction at 7 separate time points: immediately before differentiation and then daily for 6 days.

\subsection{Mouse Muscle Injury}
Muscle injury and regeneration were induced in mice using cardiotoxin. Briefly, 6 week-old C57BL/6 mice were unilaterally injured by intramuscular injection of 0.1 mL of 10 uM CTX into the tibialis anterior muscle. Injured tibialis anterior muscles were harvested at days 3 (n=2), 5 (n=2), 7 (n=2), 10 (n=4), 14 (n=4), and 28 (n=3) post-injury. Contralateral (uninjured) tibialis anterior muscles were also collected (n=9). Muscle tissue was snap-frozen and stored at -80 degrees Celsius.

\subsection{Human muscle biopsy processing}
Open muscle biopsies were placed in an aluminum foil envelope. 2-methylbutane (isopentane) was pre-chilled using liquid nitrogen and the aluminum foil envelopes were submerged in the isopentane for 15 seconds. After this, the samples were placed in cryovials at -80º for long-term storage. Samples collected at other institutions were shipped in dry ice to the NIH Muscle Disease Unit.

\subsection{RNA extraction}
The first step to perform RNA sequencing is to get the RNA from the samples of interest. To do so we had first to homogenize the muscle tissue using ceramic beads and then separate the RNA from the rest of the components of the cell using the TRIzol protocol. Briefly, muscle biopsies were homogenized in TRIzol using 1.4 mm ceramic bead low-binding tubes, and the RNA was extracted using the regular TRIzol protocol. Concentration and quality of the resulting RNA were assessed using standard NanoDrop and TapeStation protocols, respectively.

\subsection{RNA library preparation}
Libraries were prepared using the NeoPrepTM system according to the TruSeqM Stranded mRNA Library Prep protocol (Illumina) and sequenced using the Illumina HiSeq 2500 or 3000. This methodology had the advantage of being mostly automated, decreasing possible human errors in this step.

\subsection{RNA sequencing analysis}
\subsubsection{Demultiplexing}
Since modern DNA sequencers have a huge sequencing capacity, several samples can be run in each lane of a flow cell. Thus, the first step after we have sequenced the RNAseq library will be to separate all the different experiments in individual folders based on their characteristic index. This was performed using bcl2fastq v.2.17.1 parallelizing it to speed up the process. The output of this step is a file in .fastq format consisting of repeating blocks of 4 lines containing the annotation, sequence, comments, and quality of each one of the reads of the experiment.

\subsubsection{Sequence cleanup}
After the sequences of each sample are demultiplexed it is convenient to clean up the data to eliminate sequences with low quality and remove adapter contamination that may bias the results. This was performed using trimmomatic v.0.36.

\subsubsection{Alignment}
Once the sequences were cleaned each read had to be aligned to the reference genome in order to be able to identify what genes are expressed and how much. A technical complication for this step in RNAseq libraries is that mature RNA has been already spliced and thus, sequences containing sections of two or more exons will not be aligned to any region of the reference DNA since it includes the intronic regions. To overcome this limitation, specific RNA aligners have been developed that work by splitting reads that cannot be aligned confidently and trying to align each fragment of the split separately. This will be repeated recursively until the whole read could be aligned or it was concluded that there was no good match in the reference genome.  Limited studies have compared the multiple available tools to perform RNA alignment, concluding that STAR offers slightly better accuracy compared to other options.\cite{Baruzzo2017} Thus, to perform the alignment we used STAR v.2.5 using the hg19 (GRCh37) or the mm10 (GRCm38) reference genomes.

\subsubsection{Gene expression quantitation}
In RNA sequencing libraries the number of reads of each gene is proportional to the level of expression of that gene in the sample of origin. However, the number of reads will be influenced by the length of the gene and the total number of reads of the library. Thus, to quantify the levels of expression of each gene the resulting number of reads has to be normalized to the total number of reads of the library and the length of each gene. Two basic statistics are used to measure the level of expression of RNAseq experiments depending on how the normalization is performed, \gls{fpkm}, and \gls{tpm}. \glspl{fpkm} are calculated by dividing the read counts of each gene by the total number of read counts of the library and the result is divided by the length of the gene in kilobases. Alternatively, the \gls{tpm} is calculated by dividing first by the length of the gene and later on by the total number of reads of the library. The advantage of \gls{tpm} over \gls{fpkm} would be that the total number of \glspl{tpm} per library will be constant and this will allow us to compare the levels of expression of different genes across different samples. Alternatively, \gls{fpkm} is more popular and is valid for comparing the same gene across different samples, which is one of the most common comparisons to make.

RNAseq libraries can be single-end or pair-end depending on the number of RNA reads that are obtained from each strand of RNA (one in single-end and two at each side in pair-end). In single-end experiments, \gls{fpkm} can be also named RPKM or reads per kilobase million. However, for paired-end experiments the term read is not correct because each pair of reads identifies a fragment, and for this reason, is preferable to use the term \gls{fpkm} that will be correct both for single and pair-end experiments.

There are multiple libraries that allow quantitation of the level of expression of each gene and most of them result in similar results provided that the same gene coordinates are used. In our case, we used Stringtie v.1.3.3, which is a tool contained in the “new Tuxedo package”.\cite{Pertea2016}

\subsubsection{Quality control}
A key part of any Next Generation Sequencing pipeline is to confirm that the quality of the reads is adequate for the purposes of the analysis and there are no sequencing artifacts. This can be performed right after the sequencing, the demultiplexing, the alignment, or after the gene quantification steps. Given that examining the output of the quality control for hundreds of samples can be tedious and time-consuming, in our experience it was more efficient to perform the quality control at the end of the pipeline to capture any failures that the process could have produced. In case of problems with the sample, further steps of quality control after each one of the steps can help to identify the error. For doing the quality control we used fastqc v.0.11.2.

\subsubsection{Differential expression}
The main results of an RNAseq analysis are the expression levels of each one of the ~20,000 genes in a biological sample. However, the relevant questions usually require comparing the differences between conditions. The process of comparing the levels of expression among biologically relevant groups is called differential expression. This is, arguably, one of the most critical steps in the RNAseq pipeline given the multiple models that can be used to do the comparisons and the variety of tools available for the task. Some of these tools use the Poisson distribution, others the negative binomial or the beta-binomial while there are some that use Bayesian non-parametric models. Thus, both edgeR and DESeq2 use a variation of the Fisher exact test adopted for the negative binomial distribution returning exact p-values computed from the derived probabilities. CuffDiff uses the test statistics $T=\dfrac{E[\log(y)]}{Var[\log(y)]}$, where y is the ratio of the normalized counts between two conditions, and this ratio approximately follows a normal distribution. Hence a t-test is used to calculate the p-value for the differential expression. Finally, limma uses a moderated t-statistic to compute p-values in which both the standard error and the degrees of freedom are modified. All of these tools used the Benjamini and Hochberg approach for multiple hypothesis correction.\cite{Benjamini1995}

The results of these tools vary considerably and there is yet no gold standard on which is the best to do these types of analyses. However, there are studies comparing them showing that DESeq may have a slight advantage in terms of detection accuracy.\cite{Seyednasrollah2015} Moreover, our own benchmark tests showed that CuffDiff had a very poor correspondence with previous gene expression data in our field, while both EdgeR, DESeq, and limma had similar and good results. In conclusion, we decided to use the DESeq2 v.1.20 algorithm for differential expression.\cite{Love2014}

DESeq2 performs an internal normalization where the geometric mean is calculated for each gene across all samples. The counts for a gene in each sample are then divided by this mean. The median of these ratios in a sample is the size factor for that sample. This procedure corrects for library size and RNA composition bias, which can arise, for example, when only a small number of genes are very highly expressed in one experiment condition but not in the other.

Additionally, DESeq2 automatically detects count outliers using the Cooks' distance and removes these genes from the analysis. DESeq2 v.1.20 also performs independent filtering which maximizes the number of genes which will have a Benjamini and Hochberg-adjusted p–value\cite{Benjamini1995} less than a critical value set by default to 0.1. Removing these genes with low counts improves the detection power by making the multiple testing adjustment of the p-values less severe. To speed up the computations we prefiltered genes with a total count across conditions below 10. Since these genes would have been excluded from the analysis afterward anyways, this did not influence the calculations at all.

DESeq2 uses shrinkage estimation for dispersions and fold changes. A dispersion value is estimated for each gene through a model fit procedure. Using these estimations, the package fits a negative binomial generalized linear model for each gene and uses the Wald test for significance testing. The Wald test p-values from the subset of genes that pass the independent filtering step are adjusted for multiple testing using the procedure of Benjamini and Hochberg.\cite{Benjamini1995}

To ensure the stability of the central tendency and dispersion values of each biological group between different sections of the study, the normalization process included the totality of the samples even if that specific comparison did not include some of those samples.

We assigned equal weights to each autoantibody subgroups within \gls{dm} and  \gls{imnm} to avoid giving more importance to differentially expressed transcriptomic features of autoantibody subgroups with a higher number of biopsies at this stage of the analysis.

\subsubsection{Pathway analysis}
Once we have identified the list of genes explaining the differences between two or more conditions it is often necessary to find what are the pathways that are activated or suppressed. Doing so has the main difficulty of being based on manual annotations of the pathways. Thus, if the annotation is poor, incomplete, or contains wrong data our pathway analysis will be biased. There are several tools to perform pathway analysis, most of them basing their results on the Fisher’s exact test, the hypergeometric distribution, or a chi-square analysis of the 2x2 contingency tables comparing the number of genes observed in our dataset with the genes present in each pathway. DAVID and MetaCore have similar approaches to the problem by doing a simple frequentist assessment of genes present in the case dataset vs. genes observed in each pathway. Alternatively, IPA Pathway analysis also performs this basic analysis but it also provides information about the activation or suppression of the pathway. Another option is to use \gls{gsea}. This approach has the advantages of not needing an arbitrary threshold to define a list of “significant” genes and that it uses the actual significance of the associations of the genes with the trait being analyzed. However, it is limited to a single pathway comparison at a time and, thus, it is less suitable for pathway discovery purposes.

\subsubsection{RNAseq-based classification}
RNAseq is an extremely powerful technique that quantifies simultaneously the expression of more than twenty thousand genes. This has the potential to be used predictively to classify complex diseases based on its pathogenesis. Multiple algorithms can be used for this task, from simplistic approaches using logistic regression and simple classification trees to bagging or boosting strategies (like the random forest or AdaBoost respectively). Moreover, neural networks and other algorithms like the support vector machines can also be useful for this task. However, given that there are no good studies assessing when to use each algorithm based on the characteristics of the data, we took an agnostic approach to this problem by testing all these mathematical tools in the same training datasets and testing their accuracy in those samples that were not used for training purposes.

In short, the sample was split into a training set containing 2/3 of the observations and a test set containing the remaining 1/3. The training set was used to build the classificatory models and the testing set to evaluate the accuracy of the model. The classifiers that we tested were the linear support vector machine (SVM), the radial basis function SVM, random forest, nearest neighbor, Gaussian process, decision tree, multi-layer perceptron neural network, adaptative boosting (AdaBoost), gaussian naïve Bayes and quadratic discriminant analysis. Models were built in 2/3 random resamples of the data and tested in the remaining 1/3. The accuracy of classifying correctly each one of the myositis subsets was determined based on the mean and 95\% CI of one thousand resampling cycles.

We decided to use all the genes that were significantly different (with a cutoff q-value <0.05) in each group compared to the rest using all the samples. An alternative would have been to include the differentially expressed genes contained in the training set of each cycle. However, this approach was excessively computationally expensive. Nonetheless, to demonstrate the equivalency of these approaches, we modified our pipeline to train 100 cross-validation cycles using only the differentially expressed genes resulting from each training set and the performance of the models was equivalent using both methods.

\subsubsection{Ranking genes in classificatory models}
Oftentimes it is interesting to rank the most important features for the classification model. This can be done, independently of the algorithm, using recursive feature elimination.\cite{Guyon2002} This technique iteratively constructs new models removing the features with low weights and thus, estimating the importance of each feature for the prediction task.