Lezione 2: Struttura delle Read Illumina (PDF)

Lezione 2: struttura delle read Illumina, analisi della qualità dei dati di sequenziamento Illumina Sequenziamento Nella scorsa lezione abbiamo visto come funzionano le principali tecnologie di sequenziamento In questa lezione vedremo cosa significa, per un biologo computazionale, avere in m...

Lezione 2: struttura delle read Illumina, analisi della qualità dei dati di sequenziamento Illumina Sequenziamento Nella scorsa lezione abbiamo visto come funzionano le principali tecnologie di sequenziamento In questa lezione vedremo cosa significa, per un biologo computazionale, avere in mano i dati che vengono prodotti da un sequenziatore. Prenderemo in considerazione i risultati di sequenziamento della tecnologia Illumina e li vedremo nel dettaglio. Il sequenziamento di un campione produce diverse migliaia di reads e le inserisce in un solo file, che potrà essere utilizzato per l’analisi dei dati Come è fatto questo file? Che formato ha? Ci sono delle regole nella sua formattazione che possano farci capire quale è la qualità del sequenziamento? Risponderemo a tutte queste domande, ma prima dobbiamo parlare ancora un po’ della struttura dei nostri dati di sequenziamento: le read! Librerie e adattatori nel sequenziamento Illumina Il primo passo del sequenziamento consiste nel costruire una libreria di DNA o RNA. Una libreria è un insieme di molecole di DNA o RNA di cui si vuole effettuare un sequenziamento, affiancate da adattatori. Gli adattatori contengono una serie di sequenze necessarie per effettuare il legame delle sequenze alla piastra di sequenziamento e per le fasi di amplificazione 3’ P5 Index2 i5 Rd1 SP Inserto di DNA Rd2 SP Index1 i7 P7 5’ 5’ P5 Index2 i5 Rd1 SP Inserto di DNA Rd2 SP Index1 i7 P7 3’ Le sequenze che si formano sono costituite da: 1) adattatori, sequenze necessarie per effettuare il sequenziamento e di cui la composizione nucleotidica è, di conseguenza e 2) inserti di DNA, sequenze la cui composizione nucleotidica non è nota e deve essere ricavata Queste sequenza è una read Gli adattatori saranno rimossi dalle read, in quanto sono serviti per il sequenziamento, ma non rappresentano molecole che dovranno essere sequenziate Librerie e adattatori nel sequenziamento Illumina Il primo passo del sequenziamento consiste nel costruire una libreria di DNA o RNA. Una libreria è un insieme di molecole di DNA o RNA di cui si vuole effettuare un sequenziamento, affiancate da adattatori. Gli adattatori contengono una serie di sequenze necessarie per effettuare il legame delle sequenze alla piastra di sequenziamento e per le fasi di amplificazione 3’ P5 Index2 i5 Rd1 SP Inserto di DNA Rd2 SP Index1 i7 P7 5’ 5’ P5 Index2 i5 Rd1 SP Inserto di DNA Rd2 SP Index1 i7 P7 3’ L’inserto di DNA (in grigio) rappresenta la molecola da sequenziare Librerie e adattatori nel sequenziamento Illumina Il primo passo del sequenziamento consiste nel costruire una libreria di DNA o RNA. Una libreria è un insieme di molecole di DNA o RNA di cui si vuole effettuare un sequenziamento, affiancate da adattatori. Gli adattatori contengono una serie di sequenze necessarie per effettuare il legame delle sequenze alla piastra di sequenziamento e per le fasi di amplificazione 3’ P5 Index2 i5 Rd1 SP Inserto di DNA Rd2 SP Index1 i7 P7 5’ 5’ P5 Index2 i5 Rd1 SP Inserto di DNA Rd2 SP Index1 i7 P7 3’ L’inserto di DNA (in grigio) rappresenta la molecola da sequenziare L’adattatore è costituito da: Librerie e adattatori nel sequenziamento Illumina Il primo passo del sequenziamento consiste nel costruire una libreria di DNA o RNA. Una libreria è un insieme di molecole di DNA o RNA di cui si vuole effettuare un sequenziamento, affiancate da adattatori. Gli adattatori contengono una serie di sequenze necessarie per effettuare il legame delle sequenze alla piastra di sequenziamento e per le fasi di amplificazione 3’ P5 Index2 i5 Rd1 SP Inserto di DNA Rd2 SP Index1 i7 P7 5’ 5’ P5 Index2 i5 Rd1 SP Inserto di DNA Rd2 SP Index1 i7 P7 3’ L’inserto di DNA (in grigio) rappresenta la molecola da sequenziare L’adattatore è costituito da: Due sequenze, chiamate P5 (rossa) e P7 (verde), che consentono alla molecola di potersi legare alla piattaforma Illumina Due sequenze, chiamate ‘’Index’’. Queste sequenze servono per associare correttamente ogni read dal campione da cui è stata generata. Se vengono analizzati 6 campioni contemporaneamente, le read di ogni campione avranno una diversa coppia di sequenze Index. Attenzione: l’Index1 è sempre presente, mentre l’Index2 potrebbe essere assente Una coppia di sequenze primer, chiamate Rd1 SP e Rd2 SP (‘’Rd SP’’ = read sequence primer). Rappresentano le sequenze primer che fungono da innesco per la PCR Librerie e adattatori nel sequenziamento Illumina 3’ P5 Index2 i5 Rd1 SP Inserto di DNA Rd2 SP Index1 i7 P7 5’ 5’ P5 Index2 i5 Rd1 SP Inserto di DNA Rd2 SP Index1 i7 P7 3’ La combinazione di sequenze di adattatori, indici e primer variano a seconda del kit di sequenziamento utilizzato per preparare la libreria. Se consultiamo la guida e vediamo il kit TruSeq possiamo verificare che la struttura dell’adattatore è la seguente: Librerie e adattatori nel sequenziamento Illumina 3’ P5 Index2 i5 Rd1 SP Inserto di DNA Rd2 SP Index1 i7 P7 5’ 5’ P5 Index2 i5 Rd1 SP Inserto di DNA Rd2 SP Index1 i7 P7 3’ La combinazione di sequenze di adattatori, indici e primer variano a seconda del kit di sequenziamento utilizzato per preparare la libreria. Se consultiamo la guida e vediamo il kit TruSeq possiamo verificare che la struttura dell’adattatore è la seguente: Rd2 SP P7 P5 Rd1 SP Librerie e adattatori nel sequenziamento Illumina 3’ P5 Index2 i5 Rd1 SP Inserto di DNA Rd2 SP Index1 i7 P7 5’ 5’ P5 Index2 i5 Rd1 SP Inserto di DNA Rd2 SP Index1 i7 P7 3’ P5 Rd1 SP Librerie e adattatori nel sequenziamento Illumina 3’ P5 Index2 i5 Rd1 SP Inserto di DNA Rd2 SP Index1 i7 P7 5’ 5’ P5 Index2 i5 Rd1 SP Inserto di DNA Rd2 SP Index1 i7 P7 3’ Possiamo valutare, sempre per il kit TruSeq, tutte le possibili coppie di index Ogni index viene associato ad un campione differente, in modo tale che sequenze di campioni diversi possano essere sequenziate contemporaneamente sulla stessa macchina (si crea un pool di sequenze) Librerie e adattatori nel sequenziamento Illumina Nell’immagine ci troviamo nella situazione di sequenziare il DNA proveniente da due campioni Nel kit utilizzato, verranno scelti due diversi adattatori, che saranno diversi solo per la coppia di index utilizzata. Nell’esempio, le sequenze del campione 1 vengono associate all’adattatore azzurro, mentre le sequenze del campione 2 all’adattatore viola Le sequenze dei due campioni vengono ‘’unite’’ in un unico ‘’pool’’. Le sequenze del ‘’pool’’ vengono quindi caricate nel sequenziatore La tecnologia Illumina consente di poter ottenere la sequenza dei filamenti delle sequenze di entrambi i campioni Librerie e adattatori nel sequenziamento Illumina Una volta che il campione è stato sequenziato, come si associano le sequenze al campione di appartenza? Si sfruttano gli index delle sequenze in un processo definito demultiplexing Nella figura si mostra che, grazie a questo processo, le sequenze di uno stesso campione sono state raggruppate e riassociate al campione di appartenenza Questo processo è, generalmente, effettuato direttamente dai sequenziatori Illumina, ma in passato i bioinformatici dovevano occuparsi di questa operazione Librerie e adattatori nel sequenziamento Illumina Il sequenziatore elimina completamente gli adattatori, fornendo le read già ‘’pulite’’ Conviene comunque conoscere la struttura dell’adattatore perché: Ci potrebbero essere delle contaminazioni di adattatori e potrebbe essere necessario rimuoverli. Questo è dovuto a causa di problemi tecnici durante il sequenziamento Se si scaricano dei file online per delle analisi di dati secondari questi potrebbero non aver eseguito il demultiplexing e gli adattatori potrebbero non essere stati rimossi (soprattutto se il file è stato ottenuto con una tecnologia di sequenziamento che non andava a rimuovere gli adattatori). 3’ P5 Index2 i5 Rd1 SP Inserto di DNA Rd2 SP Index1 i7 P7 5’ 5’ P5 Index2 i5 Rd1 SP Inserto di DNA Rd2 SP Index1 i7 P7 3’ Librerie e adattatori nel sequenziamento Illumina L’intero adattatore viene rimosso dal sequenziatore Illumina. Dal manuale vediamo queste informazioni: Librerie e adattatori nel sequenziamento Illumina L’intero adattatore viene rimosso dal sequenziatore Illumina. Dal manuale vediamo queste informazioni: Sequenze utilizzate durante la procedura di preparazione della libreria Librerie e adattatori nel sequenziamento Illumina L’intero adattatore viene rimosso dal sequenziatore Illumina. Dal manuale vediamo queste informazioni: Sequenze identificate dal software per il demultiplexing e per la rimozione dell’adattatore Sequenze utilizzate durante la procedura di preparazione della libreria Librerie e adattatori nel sequenziamento Illumina L’intero adattatore viene rimosso dal sequenziatore Illumina. Dal manuale vediamo queste informazioni: - Come mai c’è un’adenina in più da rimuovere? - La sequenza che verrà rimossa (operazione di trimming) presenta un’adenina in più rispetto a quella dell’adattatore Librerie e adattatori nel sequenziamento Illumina Come mai c’è un’adenina in più da rimuovere? - Quando si legano gli adattatori si susseguono diversi passaggi - Il DNA frammentato viene tagliato dalla T4 DNA Polymerase, un’enzima con attività esonucleotidica - Una fosforilasi, la T4 Polymerase Kinase, fosforila l’estremità 5’ delle estremità tronche, riparandole - L’enzima Taq Polymerase aggiunge un’Adenina all’estremità 3’ (questa operazione è chiamata A-tailing) - Sugli Adattatori è presenta una Timina in più. Tramite l’azione della T4 Ligase le Timine all’estremità 5’ degli Adattatori possono legarsi alle Timine all’estremità 3’ delle molecole da sequenziare - Anche se non figurano nell’adattatore, le adenine saranno presenti nella sequenza finale. Per questo la sequenza da rimuovere è più ‘’lunga’’ di un’Adenina Librerie e adattatori nel sequenziamento Illumina L’intero adattatore viene rimosso dal sequenziatore Illumina. Dal manuale vediamo queste informazioni: - Come mai c’è un’adenina in più da rimuovere? - La sequenza che verrà rimossa (operazione di trimming) presenta un’adenina in più rispetto a quella dell’adattatore - Risposta: a causa dell’A-tailing Librerie e adattatori nel sequenziamento Illumina L’intero adattatore viene rimosso dal sequenziatore Illumina. Dal manuale vediamo queste informazioni: - Come mai c’è un’adenina in più da rimuovere? - La sequenza che verrà rimossa (operazione di trimming) presenta un’adenina in più rispetto a quella dell’adattatore - Risposta: a causa dell’A-tailing - La sequenza legata all’adattatore per la read 2 sembra diversa, come mai? Librerie e adattatori nel sequenziamento Illumina L’intero adattatore viene rimosso dal sequenziatore Illumina. Dal manuale vediamo queste informazioni: - Come mai c’è un’adenina in più da rimuovere? - La sequenza che verrà rimossa (operazione di trimming) presenta un’adenina in più rispetto a quella dell’adattatore - Risposta: a causa dell’A-tailing - La sequenza legata all’adattatore per la read 2 sembra diversa, come mai? - In realtà si tratta della stessa sequenza, complementare e con direzione opposta (reverse complement). Questo capita perché la sequenza viene letta in una direzione quando si effettua la lettura ILLUMINA. Questa direzione va da 5’ a 3’! Librerie e adattatori nel sequenziamento Illumina L’intero adattatore viene rimosso dal sequenziatore Illumina. Dal manuale vediamo queste informazioni: - Come mai c’è un’adenina in più da rimuovere? - La sequenza che verrà rimossa (operazione di trimming) presenta un’adenina in più rispetto a quella 5’ -> 3’ dell’adattatore - Risposta: a causa dell’A-tailing 5’ -> 3’ - La sequenza legata all’adattatore per la read 2 sembra diversa, come mai? 5’ -> 3’ - In realtà si tratta della stessa sequenza, complementare e con direzione opposta (reverse complement). Questo capita perché la sequenza viene letta in una direzione quando si effettua la lettura 3’ -> 5’ ILLUMINA. Questa direzione va da 5’ a 3’! Librerie e adattatori nel sequenziamento Illumina L’intero adattatore viene rimosso dal sequenziatore Illumina. Dal manuale vediamo queste informazioni: - Come mai c’è un’adenina in più da rimuovere? - La sequenza che verrà rimossa (operazione di trimming) presenta un’adenina in più rispetto a quella 5’ -> 3’ dell’adattatore - Risposta: a causa dell’A-tailing 5’ -> 3’ * - La sequenza legata all’adattatore per la read 2 sembra diversa, come mai? 5’ -> 3’ - In realtà si tratta della stessa sequenza, complementare e con direzione opposta (reverse complement). Questo capita perché la sequenza viene letta in una direzione quando si effettua la lettura 3’ -> 5’ ILLUMINA. Questa direzione va da 5’ a 3’! Dato che una read viene polimerizzata sempre in direzione 5’ -> 3’ il software * utilizza il reverse complement della sequenza utilizzata per la preparazione della libreria Sequenziamento Single-End vs Paired-End Prima di vedere come sono strutturati i file, bisogna fare una distinzione su due tipi di tecnologie di sequenziamento: single-end e paired-end Single-end: i frammenti vengono sequenziati a partire da una sola estremità. Il risultato del sequenziamento sarà costituito da un solo file di read. Paired-end: frammenti vengono sequenziati da entrambe le estremità. Il risultato del sequenziamento è costituito da due file. Un file contiene le read del primo filamento, definite read forward, mentre l’altro file contiene le read del secondo filamento, definite read reverse. Il sequenziamento a paired-end reads è più sensibile e accurato di quello a single-end reads, perché facilita notevolmente alcune operazioni, come l’assemblaggio dei genomi (lo vedremo nella prossima lezione) Sequenziamento Paired-End https://firas.io/post/illumina_sbs/ Sequenziamento Single-End vs Paired-End Come viene sequenziato l’index per poter identificare il campione da cui proviene la read? Link: https://support.illumina.com/content/dam/illumina- support/documents/documentation/system_documentation/miseq/indexed-sequencing-overview- guide-15057455-08.pdf Sequenziamento Single-End vs Paired-End Come viene sequenziato l’index per poter identificare il campione da cui proviene la read? Link: https://support.illumina.com/content/dam/illumina- support/documents/documentation/system_documentation/miseq/indexed-sequencing- overview-guide-15057455-08.pdf 1. Read 1—Read 1 follows the standard Read 1 sequencing protocol using SBS reagents. The Read 1 sequencing primer is annealed to the template strand during the cluster generation step. 2. Index Read preparation—The Read 1 product is removed and the Index 1 (i7) sequencing primer is annealed to the same template strand. 3. Index 1 (i7) Read—Following Index Read preparation, the Index 1 (i7) Read performs up to 20 cycles of sequencing. The maximum number of cycles in each Index Read depends on the system and run parameters. Sequenziamento Single-End vs Paired-End Come viene sequenziato l’index per poter identificare il campione da cui proviene la read? Link: https://support.illumina.com/content/dam/illumina- support/documents/documentation/system_documentation/miseq/indexed-sequencing- overview-guide-15057455-08.pdf 4. Index 2 (i5) Read—The Index 1 (i7) Read product is removed and the template anneals to the grafted P5 primer on the surface of the flow cell. The run proceeds through an additional seven chemistry-only cycles (no imaging occurs), followed by up to 20 cycles of sequencing. 5. Read 2 resynthesis—The Index Read product is removed and the original template strand is used to regenerate the complementary strand. The original template strand is then removed to allow hybridization of the Read 2 sequencing primer. 6. Read 2—Read 2 follows the standard paired-end sequencing protocol using SBS reagents. Link utili per approfondire: https://seekdeep.brown.edu/illumina_paired_info.html https://www.cytivalifesciences.com/en/us/news- center/enzymes-in-ngs-library-prep-10001 Organizzazione dei dati: formato fastq – dettagli ILLUMINA Organizzazione dei dati: formato fastq Per ogni campione sono forniti due file in formato fastq: un file contiene tutte le read forward e l’altro file contiene tutte le read reverse In un formato fastq ogni sequenza è caratterizzata da 4 righe Organizzazione dei dati: formato fastq Per ogni campione sono forniti due file in formato fastq: un file contiene tutte le read forward e l’altro file contiene tutte le read reverse. In un formato fastq ogni sequenza è caratterizzata da 4 righe In questo caso sono presenti 12 righe  3 sequenze! La quarta riga contiene informazioni sulla qualità delle sequenze: questa informazione può essere utilizzata dai software per poter effettuare delle operazioni di ‘’controllo qualità’’ Qualità del sequenziamento Per ogni nucleotide viene calcolata la probabilità (indicata con la lettera P) che sia stato sequenziato in maniera errata. Tale probabilità viene calcolata con il Phred Score (indicato con la lettera Q), già precedentemente utilizzato per valutare la qualità di un sequenziamento Sanger 𝑄 = −10 𝑙𝑜𝑔10 𝑃 Philip Green Qualità del sequenziamento Per ogni nucleotide viene calcolata la probabilità (indicata con la lettera P) che sia stato sequenziato in maniera errata. Tale probabilità viene calcolata con il Phred Score (indicato con la lettera Q), già precedentemente utilizzato per valutare la qualità di un sequenziamento Sanger 𝑄 = −10 𝑙𝑜𝑔10 𝑃 𝑃 = 10−𝑄Τ10 Philip Green Qualità del sequenziamento Per ogni nucleotide viene calcolata la probabilità (indicata con la lettera P) che sia stato sequenziato in maniera errata. Tale probabilità viene calcolata con il Phred Score (indicato con la lettera Q), già precedentemente utilizzato per valutare la qualità di un sequenziamento Sanger 𝑄 = −10 𝑙𝑜𝑔10 𝑃 𝑃 = 10−𝑄Τ10 Philip Green Accuratezza in Valore Q Valore P Valore P in percentuale percentuale 10 0,1 10 90 20 0,01 1 99 30 0,001 0,1 99,9 40 0,0001 0,01 99,99 50 0,00001 0,001 99,999 60 0,000001 0,0001 99,9999 https://en.wikipedia.org/wiki/Phred_quality_score Qualità del sequenziamento Per ogni nucleotide viene calcolata la probabilità (indicata con la lettera P) che sia stato sequenziato in maniera errata. Tale probabilità viene calcolata con il Phred Score (indicato con la lettera Q), già precedentemente utilizzato per valutare la qualità di un sequenziamento Sanger Per ogni nucleotide di ogni 𝑄 = −10 𝑙𝑜𝑔10 𝑃 sequenza ho a disposizione un valore che mi indica la qualità del sequenziamento! 𝑃 = 10−𝑄Τ10 Philip Green Accuratezza in Valore Q Valore P Valore P in percentuale percentuale 10 0,1 10 90 20 0,01 1 99 30 0,001 0,1 99,9 40 0,0001 0,01 99,99 50 0,00001 0,001 99,999 60 0,000001 0,0001 99,9999 https://en.wikipedia.org/wiki/Phred_quality_score Il formato fastq Il formato fastq Riga 1) Inizia sempre con @, è un’identificativo della sequenza e riporta le caratteristiche del sequenziamento Riga 2) Sequenza nucleotidica Riga 3) Inizia sempre con il simbolo ‘’+’’. Potrebbe contenere una descrizione, ma non è quasi mai presente Riga 4) Qualità delle basi (Phred Score) codificato in ASCII. Si tratta di un valore probabilistico che indica la probabilità che la base nucleotidica sia sequenziata correttamente. Invece di riportare il valore probabilistico si riporta un solo carattere La qualità delle basi è codificata con un singolo carattere del codice ASCII. In questo modo la qualità dell’i-esima base sulla riga 2 è codificata dall’i-esimo carattere della riga 4. Esempio: quale è la sesta base? Quale è la sua qualità? Il formato fastq Riga 1) Inizia sempre con @, è un’identificativo della sequenza e riporta le caratteristiche del sequenziamento Riga 2) Sequenza nucleotidica Riga 3) Inizia sempre con il simbolo ‘’+’’. Potrebbe contenere una descrizione, ma non è quasi mai presente Riga 4) Qualità delle basi (Phred Score) codificato in ASCII. Si tratta di un valore probabilistico che indica la probabilità che la base nucleotidica sia sequenziata correttamente. Invece di riportare il valore probabilistico si riporta un solo carattere La qualità delle basi è codificata con un singolo carattere del codice ASCII. In questo modo la qualità dell’i-esima base sulla riga 2 è codificata dall’i-esimo carattere della riga 4. Esempio: quale è la sesta base? Quale è la sua qualità? Base posizione 6 = Guanina  Qualità in posizione 6 = F Ma quale è il valore di probabilità associato? Il formato fastq Il formato fastq Quindi la base in posizione 6 è la Guanina Il valore di Phred Score riportato è F Dalla tabella del codice ASCII abbiamo visto che il valore associato ad F è… 37! Il formato fastq Quindi la base in posizione 6 è la Guanina Il valore di Phred Score riportato è F Dalla tabella del codice ASCII abbiamo visto che il valore associato ad F è… 37! 𝑃 = 10−𝑄Τ10 Il formato fastq Quindi la base in posizione 6 è la Guanina Il valore di Phred Score riportato è F Dalla tabella del codice ASCII abbiamo visto che il valore associato ad F è… 37! 𝑃 = 10−𝑄Τ10 𝑃 = 10−37Τ10 = 10−3.7 = 0.0001995262 Questa è la probabilità che la Guanina sia stata sequenziata non correttamente. In percentuale il valore è… 0.0001995262 * 100 = 0.01995262 Il formato fastq Quindi la base in posizione 6 è la Guanina Il valore di Phred Score riportato è F Dalla tabella del codice ASCII abbiamo visto che il valore associato ad F è… 37! 𝑃 = 10−𝑄Τ10 𝑃 = 10−37Τ10 = 10−3.7 = 0.0001995262 Questa è la probabilità che la Guanina sia stata sequenziata non correttamente. In percentuale il valore è… 0.0001995262 * 100 = 0.01995262 La probabilità (in percentuale) che la base sia stata sequenziata correttamente è: 100 − 0.01995262 = 99.98005 Il Phred Score – codifica in codice ASCII ATTENZIONE ALLA CODIFICA! https://en.wikipedia.org/wiki/FASTQ_format Esercitazione 1: comandi bash per analizzare un file fastq apriamo il terminale e creiamo una cartella: mkdir genomica_computazionale scarichiamo la cartella «Fastq_files.zip»: wget http://beagle.bio.uniroma2.it/Fastq_files.zip spostiamola in./genomica_computazionale con il comando mv mv Fastq_files.zip genomica_computazionale/ Esercitazione 1: comandi bash per analizzare un file fastq spostiamoci nella cartella con il comando: cd genomica_computazionale decomprimiamo la cartella con il comando: unzip Fastq_files.zip spostiamoci nella cartella InputData/ Esercitazione 1: comandi bash per analizzare un file fastq Nella directory sono presenti tre campioni di sequenziamento di genoma di Salmonella enterica, sequenziati con tecnologia Paired-End Aprendo un terminale Ubuntu ci muoviamo nella directory contenente i file e usiamo il comando ls -la per visualizzare i file Esercitazione 1: comandi bash per analizzare un file fastq Nella directory sono presenti tre campioni di sequenziamento di genoma di Salmonella enterica, sequenziati con tecnologia Paired-End Aprendo un terminale Ubuntu ci muoviamo nella directory contenente i file e usiamo il comando ls -la per visualizzare i file Domanda: come mai abbiamo 6 file per 3 campioni? Esercitazione 1: comandi bash per analizzare un file fastq Nella directory sono presenti tre campioni di sequenziamento di genoma di Salmonella enterica, sequenziati con tecnologia Paired-End Aprendo un terminale Ubuntu ci muoviamo nella directory contenente i file e usiamo il comando ls -la per visualizzare i file Domanda: come mai abbiamo 6 file per 3 campioni? Esercitazione 1: comandi bash per analizzare un file fastq Nella directory sono presenti tre campioni di sequenziamento di genoma di Salmonella enterica, sequenziati con tecnologia Paired-End Aprendo un terminale Ubuntu ci muoviamo nella directory contenente i file e usiamo il comando ls -la per visualizzare i file Domanda: come mai abbiamo 6 file per 3 campioni? Perché sono stati sequenziati con tecnologia paired-end! Esercitazione 1: comandi bash per analizzare un file fastq Nella directory sono presenti tre campioni di sequenziamento di genoma di Salmonella enterica, sequenziati con tecnologia Paired-End Aprendo un terminale Ubuntu ci muoviamo nella directory contenente i file e usiamo il comando ls -la per visualizzare i file Esercitazione 1: comandi bash per analizzare un file fastq ISZ-2_S5_L001_R1_001.fastq file forward ISZ-2_S5_L001_R2_001.fastq file reverse Per ogni campione è presente una coppia di file: Identificativo Un file contiene le read forward campione Un file contiene le read reverse Esercitazione 1: comandi bash per analizzare un file fastq E’ possibile ricavare il numero di reads presenti in un file, senza aprirlo? Esercitazione 1: comandi bash per analizzare un file fastq E’ possibile ricavare il numero di reads presenti in un file, senza aprirlo? Esercitazione 1: comandi bash per analizzare un file fastq E’ possibile ricavare il numero di reads presenti in un file, senza aprirlo? I file sono in formato.gz (sono dei file compressi!) Possiamo decomprimerli utilizzando il comando gunzip –k nome.file.fastq.gz Esercitazione 1: comandi bash per analizzare un file fastq E’ possibile ricavare il numero di reads presenti in un file, senza aprirlo? I file sono in formato.gz (sono dei file compressi!) Possiamo decomprimerli utilizzando il comando gunzip –k nome.file.fastq.gz File decompresso Esercitazione 1: comandi bash per analizzare un file fastq E’ possibile ricavare il numero di reads presenti in un file, senza aprirlo? I file sono in formato.gz (sono dei file compressi!) Possiamo decomprimerli utilizzando il comando gunzip –k nome.file.fastq.gz File decompresso Per aprirli tutti con un solo comando si poteva usare la seguente sintassi: gunzip –k *.fastq.gz Esercitazione 1: comandi bash per analizzare un file fastq E’ possibile ricavare il numero di reads presenti in un file, senza aprirlo? I file sono in formato.gz (sono dei file compressi!) Il parametro –k mantiene il Possiamo decomprimerli utilizzando il comando gunzip –k nome.file.fastq.gz file originale Per aprirli tutti con un solo comando si poteva usare la seguente sintassi: gunzip –k *.fastq.gz Esercitazione 1: comandi bash per analizzare un file fastq E’ possibile ricavare il numero di reads presenti in un file, senza aprirlo? Utilizzando il comando wc –l è possibile contare il numero di righe di un file Dato che il numero di righe in un file fastq dipende dal numero di sequenze nel file è possibile dividere tale numero per ricavare il numero delle sequenze (basta dividerlo per 4!) Esercitazione 1: comandi bash per analizzare un file fastq E’ possibile ricavare il numero di reads presenti in un file, senza aprirlo? Utilizzando il comando wc –l è possibile contare il numero di righe di un file Dato che il numero di righe in un file fastq dipende dal numero di sequenze nel file è possibile dividere tale numero per ricavare il numero delle sequenze (basta dividerlo per 4!) Quindi: 3011452/4 = 752863 reads! Esercitazione 1: comandi bash per analizzare un file fastq Quindi: 3011452/4 = 752863 reads! Posso inserire tutta l’operazione in un codice bash wc -l ISZ-2_S5_L001_R2_001.fastq | awk '{print $1}' Esercitazione 1: comandi bash per analizzare un file fastq Quindi: 3011452/4 = 752863 reads! Posso inserire tutta l’operazione in un codice bash wc -l ISZ-2_S5_L001_R2_001.fastq | awk '{print $1}' pipe (simbolo) awk (linguaggio) Prende l’output del comando Linguaggio utilizzato per manipolare i testi. La wc –l e lo trasforma sintassi, in questo caso, indica di prendere direttamente nell’input del una stringa, dividerla in automatico per gli programma ‘awk’ spazi e prendere l’elemento in posizione 1 Esercitazione 1: comandi bash per analizzare un file fastq Quindi: 3011452/4 = 752863 reads! Posso inserire tutta l’operazione in un codice bash wc -l ISZ-2_S5_L001_R2_001.fastq | awk '{print $1}' pipe (simbolo) awk (linguaggio) Prende l’output del comando Linguaggio utilizzato per manipolare i testi. La wc –l e lo trasforma sintassi, in questo caso, indica di prendere direttamente nell’input del una stringa, dividerla in automatico per gli programma ‘awk’ spazi e prendere l’elemento in posizione 1 Il risultato consiste nel solo numero e non riporta il nome del file Esercitazione 1: comandi bash per analizzare un file fastq Quindi: 3011452/4 = 752863 reads! Posso inserire tutta l’operazione in un codice bash rows=$(wc -l ISZ-2_S5_L001_R2_001.fastq | awk '{print $1}’) Esercitazione 1: comandi bash per analizzare un file fastq Quindi: 3011452/4 = 752863 reads! Posso inserire tutta l’operazione in un codice bash rows=$(wc -l ISZ-2_S5_L001_R2_001.fastq | awk '{print $1}’) nome della variabile sintassi necessaria per inserire il risultato di un comando bash in una variabile expr $rows / 4 Effettua la divisione Esercitazione 1: comandi bash per analizzare un file fastq Quindi: 3011452/4 = 752863 reads! Posso inserire tutta l’operazione in un codice bash rows=$(wc -l ISZ-2_S5_L001_R2_001.fastq | awk '{print $1}’) nome della variabile sintassi necessaria per inserire il risultato di un comando bash in una variabile expr $rows / 4 Effettua la divisione Esercitazione 1: comandi bash per analizzare un file fastq Abbiamo valutato un solo file… possiamo inserire l’intera procedura in un ciclo for per valutare tutti i file presenti in una directory for file in `ls./*.fastq`; do row=$( wc -l $file | awk '{print $1}' ) seq=$( expr $row / 4 ) printf "Filename: %s\tRead number: %s\n" $file $seq done Esercitazione 1: comandi bash per analizzare un file fastq Abbiamo valutato un solo file… possiamo inserire l’intera procedura in un ciclo for per valutare tutti i file presenti in una directory Per ogni file fastq presente nella directory for file in `ls./*.fastq`; do row=$( wc -l $file | awk '{print $1}' ) seq=$( expr $row / 4 ) printf "Filename: %s\tRead number: %s\n" $file $seq done Esercitazione 1: comandi bash per analizzare un file fastq Abbiamo valutato un solo file… possiamo inserire l’intera procedura in un ciclo for per valutare tutti i file presenti in una directory Per ogni file fastq presente nella directory for file in `ls./*.fastq`; do row=$( wc -l $file | awk '{print $1}' ) seq=$( expr $row / 4 ) printf "Filename: %s\tRead number: %s\n" $file $seq done Esegui le operazioni: - Calcolo numero righe - Calcolo sequenze - Stampa a schermo il risultato Esercitazione 1: comandi bash per analizzare un file fastq Abbiamo valutato un solo file… possiamo inserire l’intera procedura in un ciclo for per valutare tutti i file presenti in una directory for file in `ls./*.fastq`; do row=$( wc -l $file | awk '{print $1}' ) seq=$( expr $row / 4 ) printf "Filename: %s\tRead number: %s\n" $file $seq done Esercitazione 1: comandi bash per analizzare un file fastq Abbiamo valutato un solo file… possiamo inserire l’intera procedura in un ciclo for per valutare tutti i file presenti in una directory for file in `ls./*.fastq`; do row=$( wc -l $file | awk '{print $1}' ) seq=$( expr $row / 4 ) printf "Fielname: %s\tRead number: %s\n" $file $seq done Tecnologia Paired-End: I file ‘’R1’’ e ‘’R2’’ a coppie presentano lo stesso numero di reads Esercitazione 1: comandi bash per analizzare un file fastq Stesso comando, ma eseguito solo sui file ‘’R1’’ for file in `ls./*R1*.fastq`; do row=$( wc -l $file | awk '{print $1}' ) seq=$( expr $row / 4 ) printf "Filename: %s\tRead number: %s\n" $file $seq done