16S rRNA gene sequencing sparse count matrices: a count data simulator and optimal pre-processing pipelines

The study of microbial communities has deeply changed since it was firstly introduced in the 17th century. In the late 1970s, a breakthrough in the way bacterial communities were studied was brought by the discovery that ribosomal RNA (rRNA) genes could be used as molecular markers to perform organisms classification. Some decades later, the advent of DNA sequencing technology revolutionized the study of microbial communities, permitting a culture-independent view on the overall community contained within a sample. Today, one of the most widely used approaches for microbial communities profiling is based on the sequencing of the gene that codes for the 16S subunit of prokaryotic ribosome (16S rRNA gene), that being ubiquitous to all bacteria, but having an exact DNA sequence unique to each species, is used as a sort of molecular fingerprint for assigning to each community member a taxonomic characterization. The advent of Next-Generation Sequencing (NGS) platforms ensured 16S rRNA gene sequencing (16S rDNA-Seq) an increasing growth in election rate as preferred methodology to perform microbiome studies. Despite this, the continuous development of both experimental and computational procedures for 16S rDNA-Seq caused an unavoidable lack in standardization concerning sequencing output data treatment and analysis. This is further complicated by the very peculiar characteristics that distinguish the matrix in which samples information is summarized after sequencing. In fact, the instrumental limit on the maximum number of obtainable sequences makes 16S rDNA-Seq data compositional, i.e. they are data in which the detected abundance of each bacterial species is dependent from the level of presence of other populations in the sample. Additionally, 16S rDNA-Seq-derived matrices are typically highly sparse (70-95% of null values). These peculiarities make the commonly adopted loan of bulk RNA sequencing tools and approaches inappropriate for 16S rDNA-Seq count matrices analyses. In particular, unspecific pre-processing steps, such as normalization, risk to introduce biases in case of highly sparse matrices. The main objective of this thesis was to identify optimal pipelines that filled the above gaps in order to assure solid and reliable conclusions from 16S rRNA-Seq data analyses. Among all the analysis steps included in a typical pipeline, this project was focused on the pre-processing of count data matrices obtained from 16S rDNA-Seq experiments. This task was carried out through several steps. first, state of the art methods for 16S rDNA-Seq count data pre-processing were identified performing a thorough literature search, which revealed a minimal availability of specific tools and the complete lack in the usual 16S rDNA-Seq analysis pipeline of a pre-processing step in which the information loss due to sequencing is recovered (zero-imputation). At the same time, the literature search highlighted that no specific simulators were available to directly obtain synthetic 16S rDNA-Seq count data on which perform the analysis to identift optimal pre-processing pipelines. Thus, a 16S rDNA-Seq sparse count matrices simulator that considers the compositional nature of this data was developed. Then, a comprehensive benchmark analysis of forty-nine pre-processing pipelines was designed and performed to assess currently used and most-recen tpre-processing approaches performance and to test for appropriateness in including zero-imputation step into 16S rDNA-Seq analysis framework. Overall, this thesis considers the 16S rDNA-Seq data pre-processing problem and provide a useful guide for a robust data pre-processing when performing a 16S rDNA-Seq analysis. Additionally, the simulator proposed in this work could be a spur and valuable tool for researchers involved in developing and testing bioinformatics methods, thus helping in filling the lack of specific tools for 16S rDNA-Seq data.

Lo studio delle comunità microbiche è profondamente cambiato da quando fu per la prima volta proposto nel XVII secolo. Quando il ruolo fondamentale dei microbi nel regolare e causare malattie umane divenne evidente, i ricercatori iniziarono a sviluppare una varietà di tecniche per isolare e coltivare i batteri in laboratorio con l'obiettivo di caratterizzarli e classificarli. Alla fine degli anni '70, una svolta in come venivano studiate le comunità batteriche fu apportata dalla scoperta che i geni che codificano per l'RNA ribosomale (rRNA) potevano essere utilizzati come marcatori molecolari per la classificazione degli organismi. Alcuni decenni più tardi, l'avvento della tecnologia di sequenziamento del DNA ha rivoluzionato lo studio delle comunità microbiche, consentendo una visione complessiva coltura-indipendente della comunità contenuta in un campione. Oggi, uno degli approcci più diffusi per profilazione di comunità microbiche si basa sul sequenziamento del gene che codifica per la subunità 16S del ribosoma procariotico (gene dell'rRNA 16S). Poiché il ribosoma svolge un ruolo essenziale nella vita procariotica, esso è onnipresente in tutti i batteri, ma la sua esatta sequenza di DNA è unica per ogni specie. Per questo motivo, esso viene utilizzato come una sorta di impronta molecolare per assegnare a ciascun membro della comunità una caratterizzazione tassonomica. L'avvento delle piattaforme di Next Generation Sequencing (NGS), in grado di produrre un'enorme mole di dati riducendo tempi e costi, ha assicurato alla tecnica di sequenziamento del gene rRNA 16S (16S rDNA-Seq) una crescita nel tasso di elezione come metodologia preferita per eseguire studi sul microbioma. Nonostante ciò, il continuo sviluppo di procedure sia sperimentali che computazionali per 16S rDNA-Seq ha causato una inevitabile mancanza di standardizzazione riguardo al trattamento e all'analisi dei dati di sequenziamento. Ciò è ulteriormente complicato dalle caratteristiche molto peculiari che contraddistinguono la matrice in cui tipicamente le informazioni dei campioni sono riassunte dopo il sequenziamento. Infatti, il limite strumentale sul numero massimo di sequenze ottenibili rende i dati 16S rDNA-Seq composizionali, cioè dati in cui l'abbondanza rilevata di ogni specie batterica dipende dal livello di presenza di altre popolazioni nel campione. Inoltre, le matrici derivate da 16S rDNA-Seq sono in genere molto sparse (70-95% di valori nulli). Ciò è dovuto sia alla diversità biologica tra i campioni sia alla perdita di informazione sulle specie rare durante il sequenziamento, un effetto che è fortemente dipendente sia dalla distribuzione solitamente asimmetrica delle abbondanze delle specie presenti nei microbiomi, sia dal numero di campioni sequenziati nella stessa corsa di sequenziamento (il cosiddetto livello di multiplexing). Le suddette peculiarità rendono la comunemente adottata mutuazione di tool e approcci dall’ambito del sequenziamento di tipo bulk RNA inadeguata per analisi di matrici di conte derivanti da 16S rDNA-Seq. In particolare, fasi di pre-elaborazione non specifiche, come la normalizzazione, rischiano di introdurre forti bias in caso di matrici molto sparse. L'obiettivo principale di questa tesi era quello di identificare delle pipeline di analisi ottimali che riempissero le suddette lacune al fine di ottenere conclusioni solide e affidabili dall'analisi dei dati dell'rRNA-Seq 16S. Tra tutte le fasi di analisi incluse in una tipica pipeline, questo progetto si è concentrato sulla pre-elaborazione di matrici di conte ottenute da esperimenti di 16S rDNA-Seq. Questo scopo è stato raggiunto attraverso diversi passaggi. In primo luogo, sono stati identificati metodi all'avanguardia per la pre-elaborazione dei dati di conte di 16S rDNA-Seq eseguendo un'accurata ricerca bibliografica, che ha rivelato una minima disponibilità di strumenti specifici e la completa mancanza nella consueta pipeline di analisi 16S rDNA-Seq di una fase di pre-elaborazione in cui venga recuperata la perdita di informazioni dovuta al sequenziamento (zero-imputation). Allo stesso tempo, la ricerca bibliografica ha evidenziato che non erano disponibili simulatori specifici per ottenere direttamente dati di conte 16S rDNA-Seq sintetici su cui eseguire l'analisi per identificare pipeline di pre-elaborazione ottimali. Di consequenza, è stato sviluppato un simulatore di matrici di conte sparse derivanti da 16S rDNA-Seq che considera la natura composizionale di questi dati. In seguito, un'analisi comparativa completa di quarantanove pipeline di pre-elaborazione è stata progettata ed eseguita con lo scopo di valutare le prestazioni degli approcci di pre-elaborazione più comunemente utilizzati e più recenti e per verificare l’appropriatezza dell’inclusione di una fase di zero-imputation nel contesto delle analisi di 16S rDNA-Seq. Nel complesso, questa tesi considera il problema della pre-elaborazione dei dati provenienti da 16S rDNA-Seq e fornisce una guida utile per una pre-elaborazione dei dati robusta quando durante un'analisi 16S rDNA-Seq. Inoltre, il simulatore proposto in questo lavoro potrebbe essere uno stimolo e uno strumento prezioso per i ricercatori coinvolti nello sviluppo e nel test dei metodi di bioinformatica, contribuendo così a colmare la mancanza di strumenti specifici per i dati di rDNA-Seq 16S.