Improving the RNA-Seq analysis pipeline: read alignment and expression level quantification

DNA and RNA play an essential role in the life of each living organism. The two molecules have different characteristics and properties but their functions are strictly related. DNA encodes all the genetic instructions needed by the main cell activities in the so-called genome. DNA is related to RNA through the gene expression process, which transcribes the information encoded by DNA into RNAs. Opposite to the static information provided by DNA, the set of transcribed RNAs at a specific instant represents the current state of each cell and, at the end, it provides a dynamic characterization of its activity. For this reason, transcriptome analysis represents a powerful tool to identify the dynamic behavior of an organism, such as the response to environmental stimuli or the pathological mechanisms involved in diseases. In recent years, transcriptomic analyses were revolutionized by the advent of RNA sequencing (RNA-Seq), a new methodology that applies current Next Generation Sequencing (NGS) techniques to RNA molecules. RNA-Seq enables to investigate at high resolution all the RNA species present in a sample, characterizing their sequences and quantifying their abundances at the same time. In practice, millions of short transcript sub-sequences, called reads, are sequenced from random positions of the input RNAs using the same NGS platforms employed in DNA sequencing. Unfortunately, no information is provided about which transcripts have generated the reads or from which part of the transcripts they come from. For this reason, reads represent at the same time the output of the sequencing process and the input of complex RNA-Seq data analysis pipelines. The first task in many RNA-Seq data analysis pipelines consists in identifying the relation between the sequencing output (i.e. reads) and the sequenced transcripts. The most common approach to this problem consists in aligning the reads against a reference genome. Once the reads are positioned in the genome, it is possible to infer which transcripts have generated them analyzing the read locations. The information coming from the positions and the number of reads could be employed in a wide range of downstream analyses. For example, counting the number of reads aligned to a gene could give a measure of its expression level, whereas studying which reads are located across exon junction could identify different isoforms. At first glance, these tasks may seem very simple, but the implementation of both the single steps and the whole analysis workflow are in fact complex and still not well defined. Among all the analysis steps in the pipeline, this thesis is focused on the read alignment problem. Read alignment is identified as one of the most critical steps, both for its almost ubiquitous presence in the different RNA-Seq analysis workflows and for its complexity. The study of this pivotal task was carried out through several steps. First, a complete characterization of the problem was performed, analyzing the alignment challenges both from a methodological and a computational point of view. In addition, the algorithms and data structures employed in the alignment process were analyzed together with different ways of modeling the read alignment problem. Then, state of the art methods for RNA-Seq read alignment were identified performing a thorough literature search about RNA-Seq, which revealed the presence of many available methods. At the same time, the literature search highlighted that the identification of a suitable alignment method for a specific application is challenging, mainly due to the lack of accurate comparative analyses. Thus, a comprehensive benchmark analysis of fourteen splice aware alignment methods and four splice unaware tools was designed and performed. The simulation of several datasets describing real scenarios and the definition of a comprehensive set of accuracy and efficiency metrics were performed in order to assess the different alignment methods. The assessment revealed considerable differences between methods’ performance, highlighting often a poor correlation between accuracy and popularity. Finally, the effect of the alignment accuracy on the reliability of an expression level quantification study was assessed for a subset of alignment methods. Overall, this thesis considers the RNA-Seq read alignment problem and presents a thorough characterization of its characteristics and challenges. In a fast evolving research field such as RNA-Seq, the information resulting from the assessment of state of the art methods provides some valuable guidelines for the definition of robust and accurate analysis pipelines.DNA e RNA giocano un ruolo essenziale nelle vita di ogni organismo. Le due molecole hanno differenti caratteristiche e proprietà ma le loro funzioni sono strettamente legate. Il DNA codifica nel genoma tutte le informazioni genetiche necessarie alle principali attività delle cellula. Il DNA è legato all’RNA tramite il processo della espressione genica, processo che trascrive le informazioni codificate dal DNA nel RNA. Diversamente dalle informazioni statiche fornite dal DNA, l’insieme degli RNA trascritti in un certo istante temporale rappresenta lo stato attuale di ogni cellula e fornisce una caratterizzazione dinamica della sua attività. Per questa ragione, l’analisi del trascrittoma rappresenta un potente strumento per identificare il comportamento dinamico di un organismo, come la risposta a stimoli ambientali o i meccanismi patologici alla base di diverse malattie. Negli ultimi anni, le analisi del trascrittoma sono state rivoluzionate dall’avvento dell’RNA sequencing (RNA-Seq), una nuova metodologia che applica le attuali tecnologie di sequenziamento di nuova generazione (NGS) a molecole di RNA. L’RNA-Seq consente di studiare tutte le specie di RNA presenti nel campione in esame, caratterizzando allo stesso tempo a loro sequenza nucleotidica e la loro quantità. In pratica, milioni di sotto sequenze dei trascritti, chiamate read, vengono sequenziate a partire da posizioni casuali dei trascritti presenti nel campione, utilizzando le medesime piattaforme NGS impiegate nel sequenziamento di DNA. Sfortunatamente le tecnologie NGS producono in output le sono read e nessuna informazione viene quindi fornita riguardo a quali trascritti abbiano generato le read o da quale porzione dei trascritti esse provengano. Per questo motivo le read rappresentano allo stesso tempo l’output del processo di sequenziamento e l’input di complesse pipeline di analisi dati RNA-Seq. Il primo passo in molte pipeline consiste proprio nella identificazione della relazione tra l’output del sequenziamento (le read) e i trascritti che sono stati sequenziati. L’approccio più comune alla risoluzione di questo problema è l’allineamento delle read su un genoma di riferimento. Infatti, identificando la posizione di ogni read nel genoma è possibile inferire quale trascritto la abbia originata analizzando la sua posizione all’interno dei geni. L’informazione derivante dalla posizione e dal numero di read può essere poi utilizzata in un ampio spettro di analisi. Ad esempio, il conteggio del numero di read allineate presso un gene può essere utilizzato come misura del suo livello di espressione, mentre lo studio di quali read si trovino a cavallo di una giunzione può permettere l’identificazione di diverse isoforme. A prima vista queste analisi possono sembrare semplici, ma l’implementazione sia della intera pipeline di analisi sia delle singole fasi che la compongono è invece complessa ed ancora non ben definita. Tra tutte le fasi che compongono la pipeline di analisi dati RNA-Seq, questa tesi si focalizza sulla fase di allineamento delle read. L’allineamento delle read costituisce uno dei passi più critici nella intera analisi di dati RNA-Seq, sia per la sua complessità che per la sua diffusione e presenza nella maggior parte delle pipeline di analisi utilizzate. Lo studio di questa fondamentale operazione è stato effettuato attraverso varie fasi. In primo luogo è stata effettuata una completa caratterizzazione del problema dell’allineamento, analizzando gli aspetti critici e i problemi aperti sia dal punto di vista metodologico che computazionale. In secondo luogo, gli algoritmi e le strutture dati utilizzate nel processo di allineamento sono state analizzate insieme alle diverse strategie di modellazione del problema. Successivamente, i metodi stato dell’arte per l’allineamento di read RNA-Seq sono stati individuati attraverso una approfondita analisi della letteratura, la quale ha evidenziato la presenza di molteplici metodi per la risoluzione di questo problema. Contemporaneamente, l’analisi della letteratura ha evidenziato la difficoltà nella scelta del metodo più accurato per il particolare scenario da analizzare. La difficoltà nella individuazione del corretto metodo è dovuta principalmente per la carenza in letteratura di accurate analisi comparative. Per questa ragione, il passo successivo è stato la progettazione ed esecuzione di una approfondita analisi comparativa di 14 metodi per l’allineamento splice aware e di 4 metodi per l’allineamento splice unaware. A questo scopo, è stata effettua la simulazione di diversi dati a descrizione di molteplici scenari reali. In aggiunta, sono state sviluppate diverse metriche per la valutazione della accuratezza ed efficienza dei singoli metodi analizzati. I risultati di questa analisi hanno rivelato considerevoli differenze tra le prestazioni dei singoli metodi, sottolineando spesso uno scarso legame tra popolarità e accuratezza. L’ultimo passo dello studio è stato l’analisi degli effetti delle diverse accuratezze raggiunge in fase di allineamento sulla precisione e affidabilità delle fasi successive nella pipeline di analisi. Nello specifico, sono state studiate le conseguenze dell’uso di un sottoinsieme dei metodi di allineamento sulla accuratezza della quantificazione del livello di espressione. In conclusione, questa tesi analizza il problema dell’allineamento di read RNA-Seq e presenta una approfondita descrizione delle caratteristiche e delle criticità di questa complessa fase della pipeline. In un campo di ricerca dalla veloce evoluzione come l’RNA-Seq, le informazioni risultanti dalla valutazione comparativa dei metodi stato dell’arte fornisce preziose linee guida per l’aggiornamento e la definizione di accurate e affidabili pipeline di analisi