Analisi dei dati WGS 🧬
In questo tutorial, forniamo gli strumenti e i passaggi importanti per eseguire l'assemblaggio De Novo dei dati WGS.
Dati
Analizziamo i dati WGS ottenuti tramite un sequenziamento paired-end dei genomi dei microbi presenti in campioni di latte e liquami:
sewage_R1.fastq.gz
sewage_R2.fastq.gz
milk_R1.fastq.gz
milk_R2.fastq.gz
1. Verifica della qualità delle letture:
Innanzitutto, dobbiamo esplorare le letture e ispezionarle per capire quale tipo di pre-elaborazione debba essere eseguita per migliorarne la qualità complessiva e l'usabilità. Ad esempio, gli adattatori utilizzati per il sequenziamento devono essere rimossi prima dell'analisi a valle poiché non appartengono al profilo genomico dei microbi. Le letture che utilizzeremo sono già state filtrate dagli adattatori dal sistema di sequenziamento. Tuttavia, mettiamo qui lo strumento di controllo qualità chiamato falco che può essere utilizzato come segue:
falco milk_R1.fastq.gz
Gli output di questo comando sono file .txt e .html che riassumono i diversi aspetti delle letture, come presenza/assenza di adattatori, lunghezze delle letture... Consigliamo di aprire il file .html tramite un browser, poiché contiene grafici facili da interpretare (la descrizione di ciascun grafico/metrica può essere trovata qui).
2. Filtraggio delle letture brevi:
Più lunghe sono le letture, più accurato sarà l'assemblaggio. Come regola generale, le letture più brevi di 50 pb devono essere scartate. BBTools è una suite di strumenti ampiamente utilizzata per l'analisi dei dati WGS. Qui, utilizziamo bbduk.sh per filtrare le letture brevi come segue:
bbduk.sh in=milk_R1.fastq.gz in2=milk_R2.fastq.gz out1=filtered_milk_R1.fastq.gz. out2=filtered_milk_R2.fastq.gz minlen=50
La lunghezza minima di lettura è specificata dall'argomento "minlen".
2. Rimozione del genoma dell'ospite:
I campioni di latte e liquami contengono altri genomi provenienti, ad esempio, dalle cellule delle mucche. Questi genomi eucariotici devono essere rimossi per conservare solo i genomi microbici:
1. Dobbiamo scaricare il genoma bovino di riferimento:
wget "https://hgdownload.soe.ucsc.edu/[…]/bosTau9.fa.gz" > bosTau9.fa.gz
2. Dobbiamo indicizzare il genoma della mucca usando bowtie2:
bowtie2-build --large-index --threads 30 bosTau9.fa.gz indexed_genome
bowtie2 --time --threads 30 -x indexed_genome -1 filtered_milk_R1.fastq.gz -2 filtered_milk_R2.fastq.gz -S milk.sam (da fare anche per filtered_sewage R1 e R2)
3. Convertire i file SAM (Sequence Alignment Map) in BAM (Binary Alignment Map) utilizzando samtools:
samtools view -bS milk.sam > milk.bam (da fare anche per le fognature)
4. Ordinare i file BAM perché i file ordinati vengono analizzati (letti) più velocemente, il che velocizza l'analisi a valle:
samtools sort milk.bam -o milk_sorted.bam (da fare anche per le fognature)
5. Estrarre le letture che non sono state mappate sul genoma bovino. Queste letture corrispondono al genoma microbico:
samtools view -b -f 12 -F 256 milk_sorted.bam > milk_microbial_reads.bam (da fare anche per le fognature)
-f significa includere le letture con flag 12. Entrambe le versioni della lettura (forward e reverse) non sono state mappate sul genoma bovino.
-F significa escludere le letture con flag 256: escludere le letture non mappate da un allineamento primario.
Il significato dei flag è disponibile su questo sito web.
6. Ora possiamo dividere le letture microbiche in avanti e indietro per ripristinare il formato iniziale (R1.fastq.gz e R2.fastq.gz):
samtools fastq -@ 30 milk_microbial_reads.bam -1 milk_clean_R1.fastq.gz -2 milk_clean_R2.fastq.gz -n (da fare anche per le fognature)
Assemblaggio:
Ora che i dati grezzi sono stati elaborati, possiamo iniziare l'assemblaggio delle letture. I due strumenti più utilizzati per l'assemblaggio de novo sono MEGAHIT e METASPADES. MEGAHIT è più veloce di SPADES e richiede meno memoria. D'altra parte, SPADES è più lento, ma produce assemblaggi migliori in termini di qualità. Assembleremo solo due campioni (1 di latte e 1 di liquame), quindi possiamo utilizzare METASDAPDES.
1. Eseguiamo semplicemente l'assemblaggio per le letture pre-elaborate (>= 50 pb e pulite dal genoma bovino):
metaspades.py -t 30 -1 milk_clean_R1.fastq.gz -2 milk_clean_R2.fastq.gz -o milk (da fare anche per le fognature)
2. Per quanto riguarda le letture, più lungo è il contig, migliore è l'assemblaggio. Quindi, dobbiamo filtrare i contig brevi (< 1000pb):
bbduk.sh in=scaffolds.fasta out=scaffolds_1000.fasta minlen=1000 (da fare per impalcature per latte e fognature scaffolds.fasta)
scaffolds.fasta è l'output di metaspades.py.
3. Una volta ottenuti i contig (ovvero gli scaffold), possiamo verificare la qualità dell'assemblaggio utilizzando QUAST:
quast.py scaffolds.fasta -o quality_reports (da fare per latte e fognature)
Binning
Il binning cerca di raggruppare contig simili in base al loro organismo di origine.
1. In primo luogo, dobbiamo mappare le letture sui contig:
bowtie2-build --large-index --threads 30 milk/scaffolds.fasta indexed_milk (da fare per fognature)
2. Poi, mappiamo le letture sui contig:
bowtie2 --time --threads 30 -x indexed_milk -1 milk_clean_R1.fastq.gz -2 milk_clean_R2.fastq.gz -S milk_binning.sam (da fare per fognature)
samtools view -bS milk_binning.sam > milk_binning.bam (da fare per fognature)
samtools sort --write-index milk_binning.bam -o milk_binning_sorted.bam (da fare per fognature)
3. Ora possiamo iniziare a suddividere i dati utilizzando MetaBAT:
jgi_summarize_bam_contig_depths milk_binning_sorted.bam --outputDepth milk.dep (da fare per fognature)
metabat2 -t 30 --verbose -i milk/scaffolds_1000.fasta -a milk.dep -m 1500 -o bins/milk (da fare per fognature)
Nota: il raffinamento dei bin può essere eseguito quando si utilizzano diversi algoritmi di binning. Questo metodo unisce tutti gli output di questi algoritmi per ottenere bin raffinati. In questo caso, abbiamo utilizzato un solo algoritmo di binning, quindi non è possibile eseguire il raffinamento dei bin. Cura dei bin: gli utenti possono intervenire manualmente per migliorare ulteriormente i bin raffinati utilizzando strumenti come KBase, ggKbase, Galaxy e Anvio.
4. Verificare la qualità del binning:
checkm2 predict --threads 30 --extension fa --input *.fa --output-directory quality (da fare per latte e fognature)
Manteniamo solo i bin con Completezza > 90% e Contaminazione < 5%. Ad esempio, una Completezza del 100% significa che il bin (il genoma di bozza) contiene tutti i geni che dovrebbero essere presenti. Una Contaminazione del 2% significa che il 2% dei geni presenti appartiene a un altro organismo.
5. Dereplicazione dei bin: alcuni bin possono essere molto simili in quanto rappresentano un unico organismo. Per ridurre la ridondanza, possiamo eseguire la dereplicazione dei bin per unire bin molto simili utilizzando dRep:
dRep dereplicate . -p 30 -comp 90 -con 5 -g .bins/*.fa --genomeInfo quality_report.csv (quality_report.tsv è stato creato da checkm2)
Annotazione:
Ora che abbiamo le bozze dei genomi (ovvero i bin), possiamo eseguire annotazioni tassonomiche e funzionali per determinare rispettivamente la classificazione tassonomica e le funzioni dei geni.
1. L'annotazione tassonomica può essere eseguita utilizzando GTDB-Tk:
gtdbtk classify_wf --genome_dir bins/dereplicated_genomes --extension "fa" --out_dir . --cpus 30
2. Annotazione funzionale tramite DRAM:
DRAM-setup.py prepare_databases --output_dir dram_db
DRAM.py annotate --threads 30 -i dereplicated_genomes/*.fa -o annotation
