Un promotore U6 (mU6) del topo espande i siti di targeting genomico
Una progettazione ottimale della libreria gRNA richiede un gran numero di siti target genomici accessibili. La trascrizione del gRNA è comunemente guidata dal promotore umano U6 (hU6) che si ritiene richieda guanina (G) come primo nucleotide del suo trascritto1,2,3. Nel caso in cui il primo nucleotide non è una G, è possibile sostituire il primo nucleotide con una G o aggiungere un extra G all’estremità 5′ di gRNA, con conseguente gRNA-DNA mismatch alla fine 5′. WT-SpCas9 può tollerare gRNA-DNA mismatch all’estremità 5′, così può bersaglio qualsiasi sequenza N20NGG con hU6 promotore 1,2,3. Tuttavia, altamente specifico Cas9 nucleasi come eSpCas9(1.1) e SpCas9-HF1 sono sensibili ai mismatch gRNA-DNA all’estremità 5′29. Essi possono solo bersaglio GN19NGG sequenza quando il promotore hU6 viene utilizzato, limitando la selezione del sito di destinazione.
Uno studio precedente ha dimostrato che il mouse U6 (mU6) promotore può avviare sia adenina (A) o G trascrizione30, che potrebbe potenzialmente espandere la selezione del bersaglio. Abbiamo confrontato l’attività del promotore mU6 e hU6 promotore da espressione transitoria di gRNA per l’editing del genoma in WT-SpCas9-esprimendo cellule HEK293T (Fig. 1a). Il promotore mU6 ha mostrato attività simile al promotore hU6 per 12 gRNA testati iniziati con G (Fig. 1b, c, Tabella 1; Fig. 1a supplementare). Abbiamo testato nove gRNA iniziati con A ma privi di G a 1-4 nucleotidi, evitando gRNA troncati funzionali trascritti da G (Tabella 1). Con nostra sorpresa, entrambi i promotori potrebbero promuovere l’editing del genoma con gRNA iniziati con A (Fig. 1d, e; Fig. 1b supplementare). Tra tutti i nove gRNA testati, solo un gRNA (A8) guidato da hU6 visualizzato bassa efficienza. Abbiamo testato ulteriori nove gRNA iniziato con A, ma conteneva G a 1-5 nucleotidi. In questo caso, il promotore mU6 ha mostrato un’attività generalmente superiore al promotore hU6 (Fig. 2a-c supplementare). Successivamente, abbiamo confrontato l’attività del promotore mU6 e del promotore hU6 in cellule HeLa, e hanno mostrato un’attività simile (Fig. 3a-f supplementare). Inoltre, abbiamo confrontato l’attività del promotore mU6 e del promotore hU6 in un vettore lentivirus, e hanno mostrato un’attività simile in due punti temporali, giorno 3 e giorno 5, dopo la trasduzione (Fig. 4a, b). I nostri risultati sono coerenti con uno studio molto recente che il promotore hU6 può trascrivere piccoli RNA iniziati con A31.
In seguito, abbiamo testato l’attività del promotore mU6 per l’editing del genoma con gRNA iniziati con C o T, ma modificati in A o G (Tabella 1). Dopo 3 giorni di editing del genoma, i gRNA iniziati con G hanno mostrato una maggiore attività per tre dei cinque gRNA testati (Fig. 1f; Fig. 5a supplementare), ma la differenza è stata eliminata dopo 5 giorni (Fig. 1g; Fig. 5b supplementare). Il promotore mU6 è stato scelto nello studio seguente.
Una strategia per il test high-throughput dell’attività dei gRNA
Un recente studio ha dimostrato che una strategia di coppia RNA guida-target permette il test high-throughput dell’attività dei gRNA per Cpf132. In questa strategia, le sequenze sintetizzate di RNA guida-target sono consegnate in cellule che esprimono Cas9 da lentivirus (Fig. 2a). Dopo l’editing del genoma, le sequenze bersaglio sono PCR-amplificato per il sequenziamento profondo, permettendo la misura diretta di inserzione / cancellazione (indel) tassi indotti da nucleasi Cas9. Un ulteriore vantaggio è che lentivirus preferibilmente integrare in regioni trascrizionalmente attive che sono molto più accessibili per il CRISPR / Cas9 macchinari 32,33,34, riducendo al minimo l’influenza di editing del genoma da accessibilità cromatina. Pertanto, il set di dati ottenuti da questa strategia fornisce l’opportunità di chiarire l’attività intrinseca di gRNA basato esclusivamente sulle loro caratteristiche di sequenza.
Doench et al. hanno sviluppato uno strumento online che permette la progettazione di gRNA per il knockout genico con WT-SpCas925. Questo strumento analizza un’intera sequenza di codifica genica e classifica tutti i gRNA in base all’attività e agli effetti off-target. Abbiamo usato questo strumento per progettare gRNA per lo screening della libreria. Quattro top-ranked gRNA iniziati con A o G sono stati selezionati per ogni gene (Fig. 2b). Abbiamo anche progettato gRNAs targeting microRNAs. Poiché le lunghezze delle sequenze codificanti microRNA sono molto più brevi delle regioni codificanti dei geni, abbiamo in genere progettato tre gRNA per ogni microRNA. Un totale di 80.263 oligonucleotidi che contengono gRNA e sequenze bersaglio corrispondenti (75.312 gRNA per 19.037 geni codificanti; 4951 gRNA per 1549 microRNA) sono stati sintetizzati da microarray (dati supplementari 1). Gli oligonucleotidi sono stati amplificati in PCR e clonati nei vettori lentivirus tramite assemblaggio Gibson. Analisi della libreria plasmidica da sequenziamento profondo ha rivelato che il tasso di errore (Una lettura contiene qualsiasi mutazione è stato considerato come un errore) indotta dalla sintesi oligonucleotide o amplificazione PCR a guida RNA-regione di sequenza bersaglio era 36.5%. Questa libreria di plasmidi è stata utilizzata nei seguenti esperimenti di screening in pool per profilare l’attività gRNA per WT-SpCas9, eSpCas9(1.1), e SpCas9-HF1.
La libreria è stata confezionata in lentivirus e trasdotta in cellule HEK293T che esprimono WT-SpCas9, eSpCas9(1.1), o SpCas9-HF1 ad un MOI di 0,3. Dopo 5 giorni di editing del genoma, il DNA genomico è stato estratto, e le regioni bersaglio integrato sono stati amplificati PCR per il sequenziamento profondo (Fig. 2b). Le mutazioni alle regioni di sequenza guida RNA-target possono essere indotte sia da nucleasi Cas9 o dalla costruzione della biblioteca. Se una mutazione può essere trovato nella biblioteca originale, è stato considerato come una mutazione indotta dalla costruzione della biblioteca ed escluso dall’analisi indel. Gli indel potrebbero essere rilevati dal sequenziamento profondo nei siti bersaglio integrati (Fig. 6a supplementare). Abbiamo ottenuto validi tassi indel gRNA (legge numero > 100) di 55.604 (che copre 20.211 geni), 58.167 (che copre 20.315 geni), e 56.888 (che copre 20.270 geni) per WT-SpCas9, eSpCas9(1.1), e SpCas9-HF1, rispettivamente (dati supplementari 2). Per quanto a nostra conoscenza, questo è il più grande gRNA on-target set di attività riportato finora in cellule di mammifero.
Il test di screening è stato ripetuto sperimentalmente due volte, e due repliche indipendenti hanno mostrato un alto livello di correlazione per il tasso di indel (R = 0.92 per WT-SpCas9; R = 0.89 per eSpCas9 (1.1); R = 0.91 per SpCas9-HF1, Fig. 2c). Il tasso di indel di gRNA individuali ha anche una forte correlazione tra tre nucleasi Cas9 (Fig. 2c), indicando che alcune caratteristiche di sequenza sono favoriti per questi tre nucleasi Cas9. La distribuzione delle attività di gRNA varia notevolmente, da nessuna attività a tassi di indelebilità del 100% per queste tre nucleasi Cas9 (Fig. 2d). WT-SpCas9 ha mostrato una maggiore efficienza di editing di eSpCas9(1.1) e SpCas9-HF1 nel nostro screening (Fig. 6b supplementare). Poiché si trattava di cloni derivati da una singola cellula che potrebbero influenzare le efficienze, abbiamo selezionato cinque gRNA e trasfettati insieme con le singole nucleasi Cas9 nelle cellule. WT-SpCas9 e SpCas9-HF1 ha mostrato un’attività simile, ma eSpCas9 (1.1) ha mostrato una minore attività al giorno 3 (Fig. 6c supplementare), coerente con gli studi precedenti 15,16.
È stato riportato che il DNA plasmidico residuo da procedure di confezionamento virale può contaminare le cellule trasdotte 35, con conseguente potenziale imprecisione nella misurazione dell’attività gRNA. Abbiamo progettato una coppia di primer specifici per la spina dorsale dei plasmidi per rilevare plasmidi residui, e una coppia di primer specifici per il DNA genomico lentivirus per rilevare sia plasmidi residui e lentivirus integrato nel genoma (Fig. 7a supplementare). Queste due coppie di primer hanno mostrato un’efficienza di amplificazione simile quando il DNA plasmidico è stato usato come template (Fig. 7b supplementare). Il DNA plasmidico residuo potrebbe essere rilevato in entrambi i virus non concentrati e virus concentrati durante il confezionamento del virus (Fig. 7b supplementare). Abbiamo trasdotto le cellule HEK293T ed estratto il DNA genomico al giorno 1 e al giorno 5 dopo la trasduzione, rispettivamente. Il DNA plasmidico residuo potrebbe essere rilevato in entrambi i campioni, ma le bande PCR erano molto deboli al giorno 5 con primer specifici per backbone, indicando che il plasmide residuo degradato nel tempo (Fig. 7b supplementare). Al contrario, bande molto forti potrebbero essere rilevate con primer specifici per i virus al giorno 5, indicando che i lentivirus sempre più integrati nel genoma. Questi risultati hanno suggerito che il DNA plasmidico residuo aveva solo un’influenza minima sullo schermo della libreria.
Caratteristiche della sequenza associate all’attività del gRNA
La caratterizzazione delle caratteristiche della sequenza associate all’attività del gRNA è fondamentale per lo sviluppo di strumenti di progettazione del gRNA. Questo set di dati su larga scala qui generato ci permette di valutare meglio quali caratteristiche hanno contribuito maggiormente all’attività del gRNA. Algoritmi tra cui alberi di regressione gradient-boosted e regressione Lasso sono stati utilizzati per valutare l’importanza delle caratteristiche 36. Tuttavia, alberi di regressione gradiente-boosted fornire punteggi di importanza Gini che riflettono solo il valore assoluto del contributo caratteristica, causando la perdita di informazioni riguardanti la direzione dell’effetto; regressione di lazo mostra scarsa capacità descrittiva. Fortunatamente, un algoritmo recentemente sviluppato SHAP (SHapley Additive exPlanation), un approccio unificato per spiegare l’output di qualsiasi modello di apprendimento automatico, può potenzialmente affrontare queste limitazioni37.
Abbiamo collegato XGBoost con SHAP (chiamato Tree SHAP) per valutare l’importanza di 1031 caratteristiche, comprese le caratteristiche identificate da Doench e Wong et al.25,38, e tutte le accessibilità posizione di gRNA caratteristiche struttura secondaria (dati supplementari 3, 4). Nel complesso, i punteggi previsti sono stati fortemente influenzati dalla posizione-dipendente composizione nucleotidica per tre nucleasi Cas9 (Fig. 3a-c; dati supplementari 4). Il nucleotide più favorito era G in posizione 20 (G_20), il nucleotide immediatamente adiacente alla sequenza PAM. Altre importanti caratteristiche sovrapposte nella top 20 per tre nucleasi sono la temperatura di fusione (Tm, favorita), il numero di dimeri TT (sfavorito), C_18 (favorito), l’energia libera di auto-piegamento (favorito), e G_14 (sfavorito).
Abbiamo poi valutato la composizione nucleotidica dipendente dalla posizione del 25% più alto di gRNA attivi rispetto al 25% più basso di gRNA attivi. I risultati hanno rivelato che G era generalmente favorito, e T era generalmente sfavorito (Fig. 4a-c; Dati supplementari 5). G_20 è stato fortemente favorito per tutte e tre le nucleasi Cas9, coerente con l’analisi SHAP dell’albero. Le differenze di preferenza nucleotide tra WT-SpCas9 e varianti SpCas9 sono stati osservati anche. Le differenze di nucleotide preferito in posizione 3 (C/G vs G), 9 (G vs C/G), 10 (A/G vs A/C), 14 (C vs A/T), 16 (A/C vs C), 17 (G vs A), e 18 (C/G vs C) sono state osservate tra WT-SpCas9 e eSpCas9(1.1). Le differenze di nucleotide preferito in posizione 3 (C/G vs G), 5 (G vs C/T), 7 (C/G vs G), 9 (G vs C/G), 10 (A/G vs C), 12 (A/G vs A), 14 (C vs A/T), 17 (G vs A/C), e 18 (C/G vs C) sono state osservate tra WT-SpCas9 e SpCas9-HF1. Inoltre, abbiamo analizzato la composizione nucleotidica indipendente dalla posizione con il 20% di gRNA attivi (Fig. 4d-f). Il contenuto di G (favorito) e T (sfavorito) ha fortemente influenzato l’attività del gRNA, mentre il contenuto di A e C ha leggermente influenzato l’attività del gRNA per tre nucleasi Cas9, coerentemente con i valori di SHAP dell’albero.
Performance degli algoritmi convenzionali
Oltre a generare set di dati sull’attività dei gRNA, un altro obiettivo di questo lavoro era quello di sviluppare strumenti di predizione per la progettazione di gRNA. Abbiamo valutato le prestazioni di quattro algoritmi convenzionali di previsione dell’attività di gRNA, tra cui la regressione lineare, la regressione lineare regolarizzata L2 (regressione Ridge), la regressione XGBoost, e i modelli multilayer perceptron (MLP) con set di dati generati in questo studio. Per evitare l’over-fitting, abbiamo separato casualmente il set di dati in due sottogruppi con 85% dei dati utilizzati come set di dati di formazione per addestrare i modelli, e il restante 15% tenuto fuori utilizzato per testare la capacità di generalizzazione dei modelli addestrati (Fig. 5a). Per ottenere prestazioni ottimali, le caratteristiche con alti valori di SHAP dell’albero sono state modellate negli algoritmi (Dati supplementari 4).
Di quattro algoritmi testati qui, MLP è il più predittivo, con coefficienti di correlazione Spearman di 0.8416, 0.8457, e 0.8440 per WT-SpCas9, eSpCas9(1.1), e SpCas9-HF1, rispettivamente (Fig. 5c-e; Dati supplementari 6-8). XGBoost è il secondo più predittivo, con coefficienti di correlazione Spearman di 0.8454, 0.8310, e 0.8184 per WT-SpCas9, eSpCas9 (1.1), e SpCas9-HF1, rispettivamente. Anche la regressione lineare e la regressione ridge hanno funzionato bene, ma con un punteggio di correlazione relativamente più basso. Abbiamo anche provato la regressione lasso e la regressione SVM. Ma il coefficiente di penalità per la regressione lasso è quasi zero, il che lo ha reso equivalente al modello lineare. La regressione SVM non è riuscita a finire il benchmarking sulla scala attuale del set di dati entro un tempo ragionevole (3 settimane). Sono stati abbandonati nel confronto finale.
Performance degli algoritmi di deep-learning
Studi recenti hanno dimostrato che due algoritmi basati sul deep-learning, la rete neurale convoluzionale (CNN) e la rete neurale ricorrente (RNN), sono strumenti potenti per l’analisi di sequenze di DNA/proteine39,40,41,42,43. Essi potrebbero ottenere caratteristiche utili da DNA grezzo / sequenza proteica automaticamente senza requisito di ingegneria caratteristica. CNN è stato utilizzato per prevedere l’attività gRNA per Cpf1 e WT-SpCas926,27, mentre RNN non è stato utilizzato per la previsione attività gRNA finora. Abbiamo addestrato sia CNN e RNN per la previsione di attività gRNA. Per evitare l’over-fitting, abbiamo separato in modo casuale il set di dati in tre sottogruppi con 76,5% dei dati utilizzati come set di dati di formazione per addestrare i modelli, 8,5% utilizzato come set di dati di convalida, e il restante 15% tenuto fuori utilizzato per testare la capacità di generalizzazione dei modelli formati (Fig. 5b).
RNN ha superato CNN e altri algoritmi per la previsione di attività gRNA con coefficienti di correlazione Spearman di 0.8555, 0.8491, e 0.8512 per WT-SpCas9, eSpCas9(1.1), e SpCas9-HF1, rispettivamente (Fig. 5c-e; Dati supplementari 6-8). CNN ha ottenuto prestazioni simili a XGBoost con coefficienti di correlazione Spearman di 0.8455, 0.8313, e 0.8343 per WT-SpCas9, eSpCas9(1.1), e SpCas9-HF1, rispettivamente.
Un modello integrato migliora il potere predittivo
La spina dorsale degli algoritmi di apprendimento profondo (CNN o RNN) può solo sfruttare la composizione di k-mer o le sue dipendenze39,44. Studi recenti sulla predizione delle proteine hanno dimostrato che la capacità di predizione dei modelli di deep-learning potrebbe essere potenziata dall’aggiunta di altre caratteristiche, come il peso molecolare, l’idrofobicità e la carica assoluta, che non potrebbero essere ottenute automaticamente dai modelli di deep-learning45,46. Nel nostro lavoro, caratteristiche di sequenza indiretta tra cui accessibilità posizione della struttura secondaria, stem-loop della struttura secondaria, temperatura di fusione, e il contenuto di GC sono fortemente associati con l’attività gRNA (dati supplementari 4), ma non poteva essere ottenuto da deep learning. Considerando che il modello RNN ha ottenuto le migliori prestazioni di tutti gli algoritmi, abbiamo quindi combinato queste caratteristiche biologiche con RNN per la previsione di attività gRNA. È interessante notare che l’aggiunta di queste caratteristiche a RNN aumentato il potere di previsione, con coefficienti di correlazione Spearman di 0.8670, 0.8624, e 0.8603 per WT-SpCas9, eSpCas9 (1.1), e SpCas9-HF1, rispettivamente (Fig. 5c-e). Pertanto, RNN integrato con caratteristiche biologiche (di seguito denominato RNN + Biofeature) è stato utilizzato come il modello finale per la previsione di attività gRNA (Fig. 6).