ML-DSP: Machine Learning with Digital Signal Processing for ultrafast, accurate, and scalable genome classification at all taxonomic levels

Dopo aver progettato e implementato il prototipo dello strumento di classificazione delle sequenze genomiche ML-DSP, abbiamo studiato quale tipo di normalizzazione della lunghezza e quale tipo di distanza fossero più adatti alla classificazione dei genomi con questo metodo. Abbiamo quindi condotto un’analisi completa delle varie rappresentazioni numeriche delle sequenze di DNA utilizzate in letteratura e abbiamo determinato i tre migliori risultati. Dopo aver impostato i parametri principali (metodo di normalizzazione della lunghezza, distanza e rappresentazione numerica), abbiamo testato la capacità di ML-DSP di classificare i genomi mtDNA a livelli tassonomici che vanno dal livello di dominio fino al livello di genere, e abbiamo ottenuto livelli medi di accuratezza di classificazione di >97%. Infine, abbiamo confrontato ML-DSP con altri metodi di classificazione dei genomi basati sull’allineamento e senza allineamento, e abbiamo dimostrato che ML-DSP ha raggiunto una maggiore accuratezza e una velocità significativamente maggiore.

Analisi delle distanze e degli approcci di normalizzazione della lunghezza

Per decidere quale misura di distanza e quale metodo di normalizzazione della lunghezza fossero più adatti al confronto dei genomi con ML-DSP, abbiamo utilizzato nove diversi sottoinsiemi di sequenze mtDNA complete dal nostro dataset. Questi sottoinsiemi sono stati selezionati per includere la maggior parte dei genomi mtDNA completi disponibili (dataset Vertebrati di 4322 sequenze mtDNA), così come sottoinsiemi contenenti sequenze simili, di lunghezza simile (dataset Primati di 148 sequenze mtDNA), e sottoinsiemi contenenti genomi mtDNA che mostrano grandi differenze di lunghezza (dataset Piante di 174 sequenze mtDNA).

I punteggi di accuratezza della classificazione ottenuti utilizzando le due misure di distanza considerate (Euclidea e coefficiente di correlazione di Pearson) e due diversi approcci di normalizzazione della lunghezza (normalizzazione alla lunghezza massima e normalizzazione alla lunghezza mediana) su diversi dataset sono elencati nella tabella 2. I punteggi di accuratezza di classificazione sono leggermente più alti per PCC, ma sufficientemente vicini a quelli ottenuti quando si usa la distanza euclidea da essere inconcludenti.

Nel resto di questo documento abbiamo scelto il coefficiente di correlazione di Pearson perché è indipendente dalla scala (a differenza della distanza euclidea, che è, ad es, sensibile all’offset del segnale, per cui segnali con la stessa forma ma punti di partenza diversi sono considerati dissimili), e la normalizzazione della lunghezza alla lunghezza mediana perché è economica in termini di utilizzo della memoria.

Analisi di varie rappresentazioni numeriche di sequenze di DNA

Abbiamo analizzato l’effetto sulla precisione di classificazione ML-DSP di tredici diverse rappresentazioni numeriche unidimensionali per sequenze di DNA, raggruppate come: Mappature fisse rappresentazioni numeriche del DNA (Tabella 1 rappresentazioni #1, #2, #3, #6, #7, vedi , e rappresentazioni #10, #11, #12, #13 – che sono varianti unidimensionali della rappresentazione binaria proposta in ), mappature basate su alcune proprietà fisico-chimiche dei nucleotidi (Tabella 1 rappresentazione #4, vedi , e rappresentazione #5, vedi ), e mappature basate sui valori più vicini (Tabella 2 rappresentazioni #8, #9, vedi ).

I set di dati utilizzati per questa analisi erano gli stessi della tabella 2. I classificatori di apprendimento automatico supervisionato utilizzati per questa analisi erano i sei classificatori elencati nella sezione Metodi e implementazione, con l’eccezione dei dataset con più di 2000 sequenze dove due classificatori (Subspace Discriminant e Subspace KNN) sono stati omessi perché troppo lenti. I risultati e i punteggi medi di accuratezza per tutte queste rappresentazioni numeriche, classificatori e set di dati sono riassunti nella tabella 3.

Come si può osservare dalla tabella 3, per tutte le rappresentazioni numeriche, i punteggi medi di precisione della tabella (ultima riga: media delle medie, prima sui sei classificatori per ogni set di dati, e poi su tutti i set di dati), sono alti. Sorprendentemente, anche usando una rappresentazione numerica a singolo nucleotide, che tratta tre dei nucleotidi come se fossero uguali e ne individua solo uno (“Just-A”), si ottiene un’accuratezza media del 91,9%. La migliore accuratezza, per questi set di dati, si ottiene quando si usa la rappresentazione “PP”, che produce un’accuratezza media del 92,3%.

Per gli esperimenti successivi abbiamo selezionato le tre migliori rappresentazioni in termini di punteggi di accuratezza: “PP”, “Just-A”, e rappresentazioni numeriche “Real”.

ML-DSP per tre classi di vertebrati

Come applicazione di ML-DSP utilizzando la rappresentazione numerica “PP” per le sequenze di DNA, abbiamo analizzato l’insieme dei genomi mtDNA dei vertebrati (lunghezza mediana 16.606 bp). La MoDMap, cioè la visualizzazione 3D in scala multidimensionale delle interrelazioni dei genomi come descritte dalle distanze nella matrice delle distanze, è illustrata in Fig. 3. Il dataset contiene 3740 genomi mtDNA completi: 553 genomi di uccelli, 2313 genomi di pesci e 874 genomi di mammiferi. Quantitativamente, il punteggio di precisione di classificazione ottenuto dal classificatore Quadratico SVM è stato del 100%.

Classificazione dei genomi con ML-DSP, a tutti i livelli tassonomici

Abbiamo testato la capacità di ML-DSP di classificare sequenze mtDNA complete a vari livelli tassonomici. Per ogni set di dati, abbiamo testato usando le rappresentazioni numeriche “PP”, “Just-A” e “Real”.

Il punto di partenza era il dominio Eukaryota (7396 sequenze), che è stato classificato in regni, poi il regno Animalia è stato classificato in phyla, ecc. Ad ogni livello, abbiamo scelto il cluster con il più alto numero di sequenze e poi lo abbiamo classificato nei sotto-cluster del livello tassonomico successivo. Il livello più basso classificato era la famiglia Cyprinidae (81 sequenze) nei suoi sei generi. Per ogni set di dati, abbiamo testato tutti e sei i classificatori, e il massimo di questi sei punteggi di precisione di classificazione per ogni set di dati sono riportati nella tabella 4.

Nota che, ad ogni livello tassonomico, i punteggi massimi di accuratezza di classificazione (tra i sei classificatori) per ciascuna delle tre rappresentazioni numeriche considerate sono alti, che vanno dal 91.4% al 100%, con solo tre punteggi inferiori al 95%. Poiché anche questa analisi non ha rivelato un chiaro vincitore tra le prime tre rappresentazioni numeriche, è sorta la domanda se la rappresentazione numerica che usiamo sia importante. Per rispondere a questa domanda, abbiamo eseguito due esperimenti aggiuntivi, che sfruttano il fatto che il coefficiente di correlazione di Pearson è indipendente dalla scala, e cerca solo un modello mentre confronta i segnali. Per il primo esperimento abbiamo selezionato le tre migliori rappresentazioni numeriche (“PP”, “Just-A”, e “Real”) e, per ogni sequenza in un dato dataset, una rappresentazione numerica tra queste tre è stata scelta casualmente, con uguale probabilità, per essere il segnale digitale che la rappresenta. I risultati sono mostrati sotto la colonna “Random3” nella tabella 4: il punteggio massimo di accuratezza su tutti i set di dati è del 96%. Questo è quasi lo stesso dell’accuratezza ottenuta quando è stata usata una particolare rappresentazione numerica (1% in meno, che è ben all’interno dell’errore sperimentale). Abbiamo poi ripetuto questo esperimento, questa volta scegliendo a caso una qualsiasi delle tredici rappresentazioni numeriche considerate. I risultati sono mostrati sotto la colonna “Random13” nella tabella 4, con un punteggio medio di precisione dell’88,1%.

In generale, i nostri risultati suggeriscono che tutte e tre le rappresentazioni numeriche “PP”, “Just-A”, e “Real” hanno punteggi di accuratezza di classificazione molto alti (media >97%), e anche una scelta casuale di una di queste rappresentazioni per ogni sequenza nel dataset non influenza significativamente il punteggio di accuratezza di classificazione di ML-DSP (media 96%).

Si nota anche che, oltre ad essere altamente accurato nelle sue classificazioni, ML-DSP è ultraveloce. Infatti, anche per il più grande set di dati nella Tabella 2, subphylum Vertebrata (4322 genomi mtDNA completi, lunghezza media 16.806 bp), il calcolo della matrice di distanza (che è il grosso del calcolo di classificazione) è durato meno di 5 s. Classificare un nuovo genoma mtDNA di primate ha richiesto 0,06 s quando addestrato su 148 genomi mtDNA di primate, e classificare un nuovo genoma mtDNA di vertebrato ha richiesto 7 s quando addestrato sui 4322 genomi mtDNA di vertebrato. Il risultato è stato aggiornato con un esperimento in cui QSVM è stato addestrato sui 4322 genomi completi di vertebrati nella tabella 2, e interrogato sui 694 nuovi genomi di mtDNA di vertebrati caricati su NCBI tra il 17 giugno 2017 e il 7 gennaio 2019. L’accuratezza della classificazione è stata del 99,6%, con solo tre genomi mtDNA di rettili erroneamente classificati come genomi di anfibi: Bavayia robusta, bavayia robusta della foresta – una specie di geco, NC_034780, Mesoclemmys hogei, tartaruga testa di rospo di Hoge, NC_036346, e Gonatodes albogularis, geco dalla testa gialla, NC_035153.

visualizzazioneMoDMap vs. ML-DSP risultati di classificazione quantitativa

L’ipotesi testata dai prossimi esperimenti era che l’accuratezza quantitativa della classificazione delle sequenze di DNA tramite ML-DSP sarebbe stata significativamente più alta di quella suggerita dal raggruppamento visivo dei taxa nella MoDMap prodotta con la stessa matrice di distanza a coppie.

Come esempio, la MoDMap in Fig. 4a, visualizza la matrice di distanza dei genomi mtDNA della famiglia Cyprinidae (81 genomi) con i suoi generi Acheilognathus (10 genomi), Rhodeus (11 genomi), Schizothorax (19 genomi), Labeo (19 genomi), Acrossocheilus (12 genomi), Onychostoma (10 genomi); solo i generi con almeno 10 genomi sono considerati. La MoDMap sembra indicare una sovrapposizione tra i cluster Acheilognathus e Rhodeus, che è biologicamente plausibile in quanto questi generi appartengono alla stessa sottofamiglia Acheilognathinae. Tuttavia, quando si zooma tracciando una MoDMap solo di questi due generi, come mostrato in Fig. 4b, si può vedere che i cluster sono chiaramente separati visivamente. Questa separazione è confermata dal fatto che il punteggio di precisione del classificatore SVM quadratico per il set di dati in Fig. 4b è del 100%. Lo stesso punteggio quantitativo di accuratezza per la classificazione del dataset in Fig. 4a con Quadratic SVM è del 91,8%, che intuitivamente è molto meglio di quanto suggerirebbe la MoDMap corrispondente. Questo è probabilmente dovuto al fatto che la MoDMap è un’approssimazione tridimensionale delle posizioni dei punti che rappresentano il genoma in uno spazio multidimensionale (il numero di dimensioni è (n-1), dove n è il numero di sequenze).

Detto questo, MoDMaps può ancora servire per scopi esplorativi. Per esempio, la MoDMap in Fig. 4a suggerisce che le specie del genere Onychostoma (sottofamiglia elencata “unknown” in NCBI) (giallo), possono essere geneticamente correlate alle specie del genere Acrossocheilus (sottofamiglia Barbinae) (magenta). Esplorando ulteriormente la matrice di distanza, si scopre che in effetti la distanza tra i centroidi di questi due cluster è inferiore alla distanza tra ciascuno di questi due cluster-centroidi agli altri cluster-centroidi. Questo supporta le ipotesi, basate sull’evidenza morfologica, che il genere Onychostoma appartiene alla sottofamiglia Barbinae, rispettivamente che il genere Onychostoma e il genere Acrossocheilus sono strettamente correlati. Si noti che questa esplorazione, suggerita da MoDMap e confermata dai calcoli basati sulla matrice di distanza, non avrebbe potuto essere avviata sulla base del solo ML-DSP (o altri algoritmi di apprendimento automatico supervisionato), in quanto ML-DSP predice solo la classificazione di nuovi genomi in uno dei taxa su cui è stato addestrato, e non fornisce altre informazioni aggiuntive.

Come altro punto di confronto tra le MoDMap e i risultati dell’apprendimento automatico supervisionato, la Fig. 5a mostra la MoDMap del superordine Ostariophysi con i suoi ordini Cypriniformes (643 genomi), Characiformes (31 genomi) e Siluriformes (107 genomi). La MoDMap mostra i cluster come sovrapposti, ma il classificatore SVM quadratico che classifica quantitativamente questi genomi ha una precisione del 99%. Infatti, la matrice di confusione in Fig. 5b mostra che Quadratic SVM classifica erroneamente solo 8 sequenze su 781 (ricordiamo che, per m cluster, la matrice di confusione m×m ha le righe etichettate dalle classi vere e le colonne etichettate dalle classi predette; la cella (i,j) mostra il numero di sequenze che appartengono alla vera classe i, e sono state predette essere di classe j). Questo indica che quando la rappresentazione visiva in una MoDMap mostra sovrapposizioni di cluster, questo può essere dovuto solo alla riduzione della dimensionalità a tre dimensioni, mentre ML-DSP in realtà fornisce una classificazione quantitativa molto migliore basata sulla stessa matrice di distanza.

Applicazioni ad altri set di dati genomici

I due esperimenti in questa sezione indicano che l’applicabilità del nostro metodo non è limitata alle sequenze di DNA mitocondriale. Il primo esperimento, Fig. 6a, mostra la MoDMap di tutte le 4721 sequenze complete di virus dengue disponibili in NCBI il 10 agosto 2017, classificate nei sottotipi DENV-1 (2008 genomi), DENV-2 (1349 genomi), DENV-3 (1010 genomi), DENV-4 (354 genomi). La lunghezza media di questi genomi virali completi è di 10.595 bp. Nonostante i genomi virali della dengue siano molto simili, la precisione di classificazione di questo set di dati in sottotipi, utilizzando il classificatore Quadratico SVM, è stata del 100%. Il secondo esperimento, Fig. 6b, mostra la MoDMap di 4710 genomi batterici, classificati in tre phyla: Spirochete (437 genomi), Firmicutes (1129 genomi) e Proteobatteri (3144 genomi). La lunghezza media di questi genomi batterici completi è di 104.150 bp, con la lunghezza massima di 499.136 bp e la lunghezza minima di 20.019 bp. L’accuratezza di classificazione del classificatore Quadratico SVM per questo set di dati è stata del 95,5%.

Confronto di ML-DSP con strumenti all’avanguardia basati sull’allineamento e senza allineamento

Gli esperimenti computazionali in questa sezione confrontano ML-DSP con tre metodi all’avanguardia basati sull’allineamento e senza allineamento: lo strumento basato sull’allineamento MEGA7 con allineamento usando MUSCLE e CLUSTALW , e il metodo senza allineamento FFP (Feature Frequency Profiles).

Per questa analisi delle prestazioni abbiamo selezionato tre set di dati. I primi due set di dati sono set di dati di riferimento utilizzati in altri studi di confronto delle sequenze genetiche: Il primo set di dati comprende 38 genomi virali dell’influenza, e il secondo set di dati comprende 41 sequenze complete di mtDNA dei mammiferi. Il terzo dataset, di nostra scelta, è molto più grande, composto da 4.322 sequenze complete di mtDNA di vertebrati, ed è stato selezionato per confrontare la scalabilità.

Per i metodi basati sull’allineamento, abbiamo usato la matrice di distanza calcolata in MEGA7 dalle sequenze allineate con MUSCLE o CLUSTALW. Per il FFP senza allineamento, abbiamo usato il valore predefinito di k=5 per i k-mer (un k-mer è una qualsiasi sequenza di DNA di lunghezza k; qualsiasi aumento del valore del parametro k, per il primo set di dati, ha portato a un punteggio di precisione di classificazione inferiore per il FFP). Per ML-DSP abbiamo scelto la rappresentazione numerica Integer e abbiamo calcolato l’accuratezza media di classificazione su tutti i sei classificatori per i primi due set di dati, e su tutti i classificatori tranne Subspace Discriminant e Subspace KNN per il terzo set di dati.

La tabella 5 mostra il confronto delle prestazioni (accuratezza di classificazione e tempo di elaborazione) di questi quattro metodi. Il tempo di elaborazione include tutti i calcoli, a partire dalla lettura dei set di dati fino al completamento della matrice di distanza – l’elemento comune a tutti e quattro i metodi. I tempi di elaborazione elencati non includono il tempo necessario per il calcolo degli alberi filogenetici, le visualizzazioni MoDMap o la classificazione.

Come si vede nella tabella 5 (colonne 3, 4 e 6) ML-DSP supera di gran lunga il software basato sull’allineamento MEGA7(MUSCLE/CLUSTALW) in termini di tempo di elaborazione. In termini di accuratezza, per i set di dati di riferimento dei virus e dei mammiferi più piccoli, le accuratezze medie di ML-DSP e MEGA7(MUSCLE/CLUSTALW) erano comparabili, probabilmente a causa delle piccole dimensioni del set di allenamento per ML-DSP. Il vantaggio di ML-DSP rispetto agli strumenti basati sull’allineamento è diventato più evidente per il set di dati dei vertebrati più grande, dove le accuratezze di ML-DSP e degli strumenti basati sull’allineamento non potevano nemmeno essere confrontate, poiché gli strumenti basati sull’allineamento erano così lenti che dovevano essere interrotti. Al contrario, ML-DSP ha classificato l’intero set di 4322 genomi mtDNA di vertebrati in 28 s, con una precisione media di classificazione del 98,3%. Questo indica che ML-DSP è significativamente più scalabile di MEGA7(MUSCLE/CLUSTALW), basato sull’allineamento, poiché può classificare velocemente e accuratamente insiemi di dati che gli strumenti basati sull’allineamento non possono nemmeno elaborare.

Come si vede nella tabella 5 (colonne 5 e 6), ML-DSP supera significativamente il software senza allineamento FFP in termini di precisione (precisione media di classificazione 98.3% per ML-DSP vs. 48.3% per FFP, per il grande dataset dei vertebrati), e allo stesso tempo è complessivamente più veloce.

Questo confronto indica anche che, per questi dataset, entrambi i metodi senza allineamento (ML-DSP e FFP) hanno un vantaggio schiacciante sui metodi basati sull’allineamento (MEGA7 (MUSCLE/CLUSTALW)) in termini di tempo di elaborazione. Inoltre, quando si confrontano i due metodi senza allineamento l’uno con l’altro, ML-DSP supera significativamente FFP in termini di precisione di classificazione.

Come altro angolo di confronto, la Fig. 7 mostra le MoDMaps del primo set di dati di riferimento (38 genomi di virus dell’influenza) prodotte dalle matrici di distanza generate rispettivamente da FFP, MEGA7 (MUSCLE), MEGA7 (CLUSTALW) e ML-DSP. La figura 7a mostra che con FFP è difficile osservare una separazione visiva del dataset in cluster di sottotipi. La figura 7b, MEGA7 (MUSCLE), e la figura 7c MEGA7 (CLUSTALW) mostrano sovrapposizioni dei cluster di punti che rappresentano i sottotipi H1N1 e H2N2. Al contrario, la Fig. 7d, che visualizza la matrice di distanza prodotta da ML-DSP, mostra una chiara separazione tra tutti i sottotipi.

Finalmente le Figg. 8 e 9 mostrano gli alberi filogenetici generati da ciascuno dei quattro metodi considerati. La figura 8a, l’albero generato da FFP, ha molti genomi mal classificati, cosa che ci si aspettava data la visualizzazione MoDMap della sua matrice di distanza in Fig. 7a. La figura 9a mostra l’albero filogenetico generato da MEGA7, che era lo stesso sia per MUSCLE che per CLUSTALW: ha solo un genoma H5N1 classificato erroneamente, collocato in mezzo ai genomi H1N1. Le figure 8b e 9b mostrano l’albero filogenetico generato utilizzando la distanza prodotta da ML-DSP (mostrato due volte, in parallelo con gli altri alberi, per facilitare il confronto). ML-DSP ha classificato correttamente tutti i genomi.

Discussione

L’efficienza computazionale di ML-DSP è dovuta al fatto che è privo di allineamenti (quindi non ha bisogno di allineamenti di sequenze multiple), mentre la combinazione di rappresentazioni numeriche 1D, Discrete Fourier Transform e Pearson Correlation Coefficient lo rende estremamente efficiente in termini di tempo di calcolo, e quindi scalabile.

ML-DSP non è senza limiti. Prevediamo che la necessità di sequenze di uguale lunghezza e l’uso della normalizzazione della lunghezza potrebbero introdurre problemi con l’esame di piccoli frammenti di sequenze di genomi più grandi. Di solito i genomi variano in lunghezza e quindi la normalizzazione della lunghezza comporta sempre l’aggiunta (up-sampling) o la perdita (down-sampling) di alcune informazioni. Anche se il coefficiente di correlazione di Pearson può distinguere i modelli di segnale anche in piccoli frammenti di sequenza, e non abbiamo trovato nessuno svantaggio considerevole considerando i genomi completi di DNA mitocondriale con le loro inevitabili variazioni di lunghezza, la normalizzazione della lunghezza può causare problemi quando abbiamo a che fare con i frammenti di genomi, e le sequenze di genoma nucleare molto più grandi.

Infine, ML-DSP ha due svantaggi, inerenti a qualsiasi algoritmo di apprendimento automatico supervisionato. Il primo è che ML-DSP è un metodo black-box che, pur producendo una previsione di classificazione molto accurata, non offre una spiegazione (biologica) per il suo risultato. Il secondo è che si basa sull’esistenza di un insieme di allenamento da cui attinge la sua “conoscenza”, cioè un insieme costituito da sequenze genomiche note e dalle loro etichette tassonomiche. ML-DSP utilizza tale insieme di allenamento per “imparare” a classificare nuove sequenze in una delle classi tassonomiche su cui è stato addestrato, ma non è in grado di assegnarle a un taxon a cui non è stato esposto.