Ottimizzazione dei primer PCR mirati al gene dell’RNA ribosomiale 16S batterico

Vincoli problematici

Come indicato nel paragrafo precedente, un primer-set-pair ottimale dovrebbe massimizzare simultaneamente efficienza e copertura e ridurre al minimo il matching-bias. Di seguito, descriviamo come abbiamo codificato quantitativamente questi vincoli.

Efficienza

Le coppie di primer-set-perfect devono soddisfare diversi vincoli, volti a migliorare l’efficienza e la specificità della PCR ., Tuttavia, soddisfare contemporaneamente tutti i vincoli è spesso poco pratico e la maggior parte dei primer all’avanguardia viola uno o più vincoli . Abbiamo quindi deciso di introdurre l’efficienza come punteggio di ottimizzazione, codificando molti dei vincoli come funzioni di punteggio fuzzy. Più precisamente, abbiamo definito il nostro punteggio di efficienza come la somma di dieci termini di punteggio: sette termini di punteggio fuzzy relativi ai vincoli di efficienza del primer singolo, calcolati in media su tutti i primer nelle coppie primer-set -, più tre termini di punteggio relativi all’efficienza delle coppie primer-set-nel suo complesso., Poiché tutti i termini sono destinati a variare tra 0 e 1, il punteggio di ottimizzazione varia da 0 (efficienza minima) a 10 (efficienza massima).

In linea di massima, il nostro punteggio fuzzy conta 1 per ogni vincolo che è perfettamente soddisfatto, o, in alternativa, un valore compreso tra 0 e 1 a seconda di quanto il primer sia vicino al limite del vincolo. Ad esempio, si consideri la temperatura di fusione del primer, Tm. Tm dovrebbe essere maggiore o uguale a 52 gradi in un primer perfetto , ma 51 è ancora tollerabile, anche se non ideale., In questo caso, la nostra funzione di punteggio fuzzy assegna 1 a temperature di 52 gradi o superiori, 0 a temperature di 50 gradi o meno e considera una funzione crescente lineare tra 50 e 52 gradi. Ogni termine è precisamente descritto in quanto segue.

I 7 termini del punteggio single-primer sono:

1. Temperatura di fusione: la temperatura di fusione Tm di un primer viene calcolata con la formula nearest-neighbor . Il termine del punteggio è 1 se Tm ≥ 52, 0 se Tm ≤ 50 e (Tm-50)/2 se 50 < Tm < 52.,

2. GC-content: GC-content è la frazione fGC di coppie di basi nella sequenza di primer uguale a G (guanina) o C (citosina). Il termine del punteggio è 1 se 0.5 ≤ fGC ≤ 0.7, 0 se fGC > 0.7 o fGC < 0.4 e (0.5 – fGC)/0.1 se 0.4 ≤ fGC < 0,5.

3. 3′-end stability – score term 1: sono definiti due termini di punteggio relativi alla stabilità 3′-end. Il primo termine è 0 se le ultime tre basi del primer sono costituite interamente da As (adenine) e Ts, (timine) e 1 altrimenti.,

4. 3′-end stability – score term 2: il secondo termine di punteggio è 0 se le ultime 5 basi contengono più di 3 Cs o Gs e 1 altrimenti.

5. Omopolimeri: un omopolimero è una sequenza di nucleotidi identici. Il termine del punteggio è 1 se non ci sono omopolimeri più lunghi di 4 nt, 0,5 se non ci sono omopolimeri più lunghi di 5 nt e 0 se c’è almeno un omopolimero più lungo di 5 nt nella sequenza.

6. Auto-dimeri: la presenza di regioni auto-complementari tra coppie di primer identici può portare alla generazione di auto-dimeri., Considerando il numero massimo di corrispondenze in un allineamento privo di lacune tra un primer con il suo complemento inverso, maxM, il termine del punteggio è 1 se maxM ≤ 8, 0 se maxM ≥ 11 e (11-maxM)/3 se 8 < maxM < 11.

7. Forcine: una forcina può essere formata in presenza di auto-complementarità all’interno della sequenza di primer, specialmente alla sua estremità 3′., Il termine punteggio è 0 se, per almeno un allineamento privo di gap tra il primer e il complemento inverso del suo 3 ‘ – end, sia l’ultimo nucleotide e 3 o più dei 4 nucleotidi precedenti corrispondono, e 1 altrimenti.

I 3 termini del punteggio primer-set-pairs sono definiti come segue:

1. Intervallo di temperatura di fusione: l’intervallo di temperatura di fusione ΔTm di una coppia primer-set-pair è calcolato come il massimo meno il minimo delle temperature di fusione di tutti i primer nella coppia impostata., Il termine del punteggio è 1 se ΔTm ≤ 3, 0 se ΔTm ≥ 5 e (5-ΔTm)/2 se 3 < ΔTm < 5.

2. Dimeri: consideriamo il numero massimo di corrispondenze maxM tra tutti i possibili allineamenti tra tutte le possibili combinazioni di primer avanti e indietro da una coppia di primer-set -. Il termine del punteggio è 1 se MAXM ≤ 8, 0 se MAXM ≥ 11 e (11-maxM) / 3 se 8 <maxM< 11.,

3. Gamma di lunghezza amplicone: a causa della riduzione nota dell’efficienza PCR con l’aumento della lunghezza amplicone , vogliamo che le lunghezze degli ampliconi generati si trovino in un intervallo ristretto. In particolare vogliamo evitare ampliconi molto più brevi della lunghezza target, poiché sarebbero troppo amplificati rispetto agli altri. Tuttavia, vogliamo essere in grado di tollerare una piccola frazione di valori anomali, al fine di evitare di penalizzare coppie di primer-set potenzialmente preziose a causa solo di poche sequenze rare., Dato un insieme rappresentativo di sequenze batteriche 16S, chiamato ” set di riferimento” d’ora in poi, consideriamo la differenza Δamplen tra la mediana e il primo percentile di lunghezze di amplicone attraverso tutti i possibili ampliconi, formata abbinando tutte le combinazioni di primer avanti e indietro dalla coppia di set con l’insieme di riferimento. Il termine del punteggio è 1 se Δamplen ≤ 50 nucleotidi, 0 se Δamplen ≥ 100 e (100-Δamplen)/50 se 50 < Δamplen < 100.,

La scelta dei criteri di punteggio e la soglia predefinita si basano sulla letteratura precedente . Tuttavia, sia le soglie che gli intervalli di tolleranza fuzzy possono essere impostati dall’utente in modo diverso da quello predefinito e in base alle sue esigenze sperimentali specificando i valori desiderati come parametri di input quando si chiama lo strumento da riga di comando.

Copertura

Il punteggio di copertura è definito come il numero di sequenze 16S abbinato da almeno un primer., Date le sequenze di un primer e di un 16S batterico, definiamo seed gli ultimi 5 nucleotidi all’estremità 3’di un primer e consideriamo una sequenza 16S come corrispondente al primer se esiste una regione della sequenza 16S che corrisponde esattamente i) al seme del primer; e ii) il resto del primer con al massimo 2 disallineamenti . Una sequenza 16S da un set di riferimento è considerata coperta da una coppia di primer-set-se almeno un primer avanti e uno inverso nella coppia di primer-set-set corrispondono alla sequenza., Poiché l’efficienza della PCR diminuisce con la lunghezza dell’amplicone, imponiamo un ulteriore vincolo: data una coppia di primer-set-e un set di riferimento di sequenze 16S, stimiamo la lunghezza dell’amplicone target come la mediana delle lunghezze di tutti gli ampliconi ottenuti abbinando tutte le combinazioni di primer avanti e indietro dalla coppia di primer-set-con il set di riferimento. Consideriamo quindi non coperte tutte le sequenze di riferimento 16S la cui lunghezza dell’amplicone differisce di oltre 100 nucleotidi (più o meno lunghi) dalla lunghezza del bersaglio.,

Matching-bias

Dato un set di riferimento di sequenze 16S e un primer-set-pair, il terzo punteggio di ottimizzazione misura la variabilità del numero di combinazioni di primer avanti e indietro corrispondenti a ciascuna sequenza di riferimento 16S. La variabilità di copertura dovuta al bias di corrispondenza dovrebbe essere ridotta al minimo, o almeno contabilizzata, quando lo studio ha lo scopo di quantificare le abbondanze relative delle diverse specie batteriche, a causa del bias di amplificazione verso le specie coperte da più combinazioni di primer avanti e indietro., Come misura di matching-bias, sfruttiamo il coefficiente di variazione della copertura attraverso le sequenze target, calcolato come la deviazione standard sulla media del numero di combinazioni corrispondenti a ciascuna sequenza.

Set di riferimento di sequenze 16S, preparazione e annotazione

Per ottimizzare i tre punteggi sopra, ci affidiamo a un set rappresentativo di sequenze 16S batteriche estratte da un database pubblico di sequenze 16S, GreenGenes ., Il database di sequenze GreenGenes 16S è organizzato in Unità tassonomiche operative (OTU), che sono cluster nidificati di sequenze nel database, organizzati a diversi livelli di somiglianza tra cluster. Per ogni livello di somiglianza, una sequenza di riferimento è associata a ciascun cluster, al massimo simile a tutte le altre sequenze nello stesso cluster . L’insieme di sequenze di riferimento può quindi essere considerato un sottoinsieme rappresentativo dell’intero database di sequenze, diventando sempre più accurato per aumentare i livelli di somiglianza tra cluster (e, quindi, il numero di sequenze di riferimento)., Abbiamo scelto un livello di somiglianza inter-cluster dell ‘ 85% come un buon compromesso tra rappresentatività e complessità, corrispondente a un set di 5088 sequenze rappresentative da utilizzare per valutare i criteri di ottimizzazione.

Sebbene molto sensibile nell’annotare i domini dei batteri e degli Archaea, la tassonomia dei GreenGenes non è progettata per distinguere sequenze appartenenti a eucarioti o virus., Per questo motivo, abbiamo deciso di ri-annotare le sequenze batteriche 16S sfruttando la tassonomia originale NCBI per identificare con precisione, tra le sequenze rappresentative, solo quelle appartenenti al dominio dei batteri. Poiché le informazioni sul dominio mancano dall’annotazione NCBI per circa il 20% delle sequenze, abbiamo progettato una procedura ad hoc per identificare le sequenze batteriche tra queste. La procedura è descritta in dettaglio nei Materiali supplementari (vedere il file aggiuntivo 1)., Abbiamo scelto conservativamente di considerare solo le sequenze annotate come batteri sia nella nostra annotazione curata, basata su NCBI che nell’annotazione originale di GreenGenes. Ciò ha provocato una serie di 4573 sequenze rappresentative di 16S appartenenti al dominio dei batteri.,

Algoritmo di ottimizzazione

Poiché il problema della scelta ottimale dei primer richiede l’ottimizzazione simultanea di diversi punteggi concorrenti, può essere lanciato come un problema di ottimizzazione multi-obiettivo, dove lo spazio di ricerca è l’insieme di tutte le possibili coppie di primer-set-e una funzione di punteggio, o criterio di ottimizzazione, può essere definito in modo, Quando più di un criterio deve essere ottimizzato contemporaneamente, ma gli obiettivi da ottimizzare sono contrastanti, di solito non si è interessati a una singola soluzione, ma piuttosto all’insieme di soluzioni ottimali di Pareto, cioè all’insieme di soluzioni per le quali nessuno degli obiettivi può essere migliorato senza sacrificare almeno un altro obiettivo ., Il risultato dell’ottimizzazione multi-obiettivo non è più una coppia di primer-set ottimale unica, come nell’ottimizzazione single-objective, ma piuttosto una raccolta di coppie di primer-set che non sono peggiori di qualsiasi altra coppia di primer-set e rigorosamente migliori secondo almeno uno dei criteri. Più precisamente, per il problema di ottimizzazione tri-obiettivo di massimizzare i punteggi di ottimizzazione di efficienza (E) e copertura (C) e minimizzare il punteggio di matching-bias (M), come definito nella sezione precedente, le coppie di primer-set-candidati vengono valutate secondo un vettore di funzione obiettivo f = (f E ; f C ; fM)., Dato due primer-set-pairs p e p’, diciamo che p domina p ‘(p p p’) se e solo se f (p) f f (p’), fE (p) ≥ fE (p’), fC (p) ≥ fC (p’) e fM (p) ≤ fM (p’). Se non esiste p ‘tale che p’ p p, il primer-set-pair p è chiamato Pareto-ottimale. In questo contesto, l’obiettivo della scelta ottimale dei primer è determinare (o approssimare) l’insieme di tutte le coppie Pareto-optimal primer-set, la cui immagine nello spazio tri-obiettivo è chiamata fronte di Pareto .

Per cercare il fronte di Pareto ottimale ci affidiamo all’approccio di ricerca locale di miglioramento iterato a due fasi proposto da Dubois-Lacoste et al., e sfruttato efficacemente in Sambo et al. e Borrotti et al. per la progettazione ottimale multi-obiettivo di esperimenti.

La ricerca locale parte da una soluzione iniziale e la perfeziona iterativamente applicando piccole modifiche locali e valutando ogni volta il loro effetto sulla qualità della soluzione; si ferma quando non ci sono ulteriori modifiche locali in grado di migliorare la soluzione. Il processo viene iterato da diversi punti di partenza e viene restituita la soluzione migliore mai trovata, come approssimazione dell’optimum sconosciuto ., Un’estensione comune della ricerca locale al caso multi-obiettivo è quella di iniziare da un insieme di soluzioni di Pareto iniziali, campionare una soluzione dal fronte, ottimizzare con la ricerca locale una scalarizzazione casuale del problema, cioè una combinazione lineare dei punteggi di ottimizzazione con pesi campionati uniformemente a caso dal simplesso unitario, aggiornare il fronte di Pareto e iterare fino a quando non viene soddisfatta una condizione di terminazione .,

La procedura MULTI-OBJECTIVE-SEARCH, il cui pseudo-codice è riportato in quanto segue, riceve come input l’intervallo desiderato di lunghezze di amplicone (rangeamplen), un insieme rappresentativo di sequenze 16S (repset), un insieme iniziale di coppie di primer (possibilmente degenerate) (init) e il numero di riavvii (nres). La procedura inizia selezionando da init tutte le possibili coppie di primer con la lunghezza dell’amplicone desiderata, la lunghezza del primer (tra 17 e 21 nucleotidi) e il dominio target (Batteri o universali).,

Le coppie di primer degenerati vengono convertite in coppie di primer-set-set non degenerati e aggiunte a un archivio. La procedura, quindi, scorre nrest volte, ogni volta campionamento casuale primer-set-coppia pstart da Pareto anteriore e un vettore casuale α dei pesi relativi per l’ottimizzazione dei punteggi, con pesi campionati uniformemente dall’unità simplex; la procedura, quindi, risolve un scalarization di multi-problema oggettivo, cioè un singolo problema oggettivo in cui una combinazione lineare dei tre obiettivi, con i relativi pesi α è ingrandita, e aggiunge il risultato di archivio., A questo scopo, efficienza, copertura e punteggi matching-bias sono normalizzati al loro massimo, in modo che ogni punteggio normalizzato varia tra 0 e 1, e matching-bias è ridefinito come 1 – matching-bias, in modo che possa essere massimizzato come gli altri punteggi., amplicone lunghezza in rangeamplen

2 Aggiungi ad un archivio corrispondente non-primer degenerato-set-coppie

3 r = 1 per nrest

4 pf = PARETO-FRONT(archivio)

5 Campione pstart da pf

6 Campione α 3, con Σi ai = 1

7 p = LOCAL-SEARCH(pstart , α , repset)

8 Aggiungere p per archivio

9 archivio

Unico obiettivo di ottimizzazione è ottenuto utilizzando il Miglior Miglioramento algoritmo di Ricerca Locale : a partire da un iniziale primer-set-coppia, il LOCALE-algoritmo di RICERCA dei cicli attraverso il primer della coppia e, per ciascun primer, analizza il quartiere, io.,e. l’insieme di tutte le possibili perturbazioni locali del primer. Le perturbazioni locali consistono in tutti i possibili capovolgimenti di un nucleotide (valutando le altre tre possibili basi) e tutte le possibili aggiunte e rimozioni di un nucleotide alle estremità., La ricerca nello spazio della soluzione viene eseguita con il miglior approccio di ricerca locale di miglioramento: dopo aver generato l’intero quartiere come spiegato sopra, l’algoritmo seleziona la migliore perturbazione del vicino, parte da essa per generare il prossimo quartiere e itera fino a raggiungere una soluzione per la quale non è possibile trovare una perturbazione del vicino migliore. La procedura termina quando non è possibile applicare ulteriori miglioramenti locali a qualsiasi primer nella coppia primer-set-pair., La funzione WEIGHTED-SCORE calcola i tre punteggi di ottimizzazione da una coppia di primer-set e dal set rappresentativo, moltiplica i punteggi per i pesi relativi α e restituisce la somma dei risultati.

Abbiamo sviluppato uno strumento software che implementa il nostro approccio e lo abbiamo rilasciato sotto la licenza GNU General Public Licence come strumento software mopo16S (Multi-Objective Primer Optimization for 16S experiments) pressohttp://sysbiobig.dei.unipd.it/?q=Software#mopo16S., mopo16S è implementato come uno strumento a riga di comando C++ multithreading; lo strumento software si basa sugli algoritmi efficienti e strutture dati dalla libreria SeqAn e utilizza la libreria OpenMP per il multithreading.,>

4 per i = 1 a |pcurr|

5 pri = i-esimo primer di pcurr

6 per pnew = pcurr con tutti i possibili aggiunte e rimozioni di base alle estremità e la sostituzione di una base di pri

7 scorenew = PONDERATO PUNTEGGIO(pnew , α , repset)

8 se scorenew > scorebest

9 pbest = pnew

10 scorebest = scorenew

11 pcurr = pbest

12 ritorno pcurr

Stato-of-the-art coppie di primer iniziale soluzioni

Abbiamo selezionato il database online probeBase come fonte di primer candidati-set-coppie per essere utilizzati come soluzioni iniziali da mopo16S., Il database contiene più di 500 coppie di primer (possibilmente degenerati) e riporta per ogni primer la sua sequenza, il filamento e la posizione in cui corrisponde al gene di riferimento 16S Escherichia coli e il dominio target per il quale è progettato (essendo Batteri, Archaea o Universale).,

Dato un intervallo desiderato per il target di ampliconi di lunghezza come input di mopo16S, abbiamo selezionato tutte le coppie di primer da probeBase database di soddisfare tutte le seguenti proprietà:

• Amplicone lunghezza dell’intervallo desiderato;

• Lunghezza di entrambi i primer maggiore o uguale a 17 nt e minore o uguale a 21 nt;

• i Batteri o Universale dominio di destinazione sia di primer.,

Poiché il nostro approccio è quello di lavorare con set di primer non degenerati, in caso di degenerazioni nel primer avanti o inverso, sostituiamo il primer degenerato con un corrispondente set di primer non degenerati, ottenuto assegnando tutte le possibili combinazioni di valori ai nucleotidi degenerati nel primer. Un esempio di questa procedura è riportato nella Tabella 1.

Abbiamo calcolato i tre punteggi per ciascuna delle coppie primer-set-e identificato, tra queste, le coppie primer-set-che formano il fronte iniziale di Pareto.