# [MBI] Ćwiczenie 1: Asemblacja *de novo* DNA **Data: 31.03.2022** ## Skład zespołu * Rafał Kulus (300249) * Kamil Przybyła (300254) ## Genom referencyjny $$ 300254 \mod 150 = 104 $$ ``` ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/001/506/805/GCF_001506805.1_ASM150680v1/GCF_001506805.1_ASM150680v1_genomic.fna.gz ``` ## Generowanie odczytów Przygotowano katalog na volume'y kontenerów o nazwie `mbi_lab1`, w którym umieszczony został plik z genomem referencyjnym. Wszystkie kontenery korzystały z tego samego katalogu. ```console $ docker run --rm -v ~/mbi_lab1:/tmp -w /tmp wkusmirek/pirs pirs simulate -x 50 -m 400 -v 20 -l 100 --subst-error-rate=0.01 GCF_001506805.1_ASM150680v1_genomic.fna ``` Wyjście: ``` [pIRS] Bases in reference sequences: 1643938 [pIRS] Read pairs simulated: 410984 [pIRS] Bases in reads: 82196800 [pIRS] Coverage: 50.00 [pIRS] Substitution error count: 1104830 [pIRS] Average substitution error rate: 1.344% [pIRS] Insertion count: 365 [pIRS] Deletion count: 864 [pIRS] Average insertion rate: 0.00044% [pIRS] Average deletion rate: 0.00105% [pIRS] Average insertion length: 1.08 [pIRS] Average deletion length: 1.03 [pIRS] Fragments affected by GC bias: 7.97% [pIRS] Bases masked by EAMSS algorithm: 0 ``` ### Q&A: **Ile zostało wygenerowanych odczytów? Jakiej długości?** Wygenerowano 821968 odczytów (410984 par) o łącznej długości 82196800. Każdy z odczytów składał się z 100 nukleotydów. **Oblicz wygenerowaną głębokość pokrycia genomu odczytami. Czy wynik jest przybliżony do zakładanego poziomu 50x?** Tak, bo $$ \frac{82196800}{1643938} \approx 49.99 $$ **W jaki sposób można znaleźć odczyty, które zawierają błędy?** W pliku `Sim_100_400.read.info` znajduje się informacja: ``` # This file is a log of every read that was simulated. It shows # exactly where each read came from and the substitution errors, # insertions, and deletions (if any) that were made to it. ``` Opis przykładowego odczytu: ``` read_400_1/2 GCF_001506805.1_ASM150680v1_genomic.fna NZ_CP010506.1 Campylobacter jejuni strain CJ677CC527 chromosome, complete genome 522281 - 431 0 51,C->T;56,C->G; - - ``` Ostatnie 3 kolumny zawierają informacje o podstawieniach, wstawieniach i usunięciach. **Odczytaj z wygenerowanych plików odległośći pomiędzy sparowanymi odczytami. Czy wartości zgadzają się z ustawianymi parametrami aplikacji?** W pliku `Sim_100_400.insert_len.distr` znajduje się informacja: ``` # This file shows the length distribution of the simulated inserts. # We were trying to simulate inserts with a mean length of 400 and a # standard deviation of 20. The actual mean is 399.524, and the actual # standard deviation is 19.987. ``` A zatem wartości są bardzo zbliżone do ustawionych parametrów aplikacji. ## Asemblacja de novo ```console $ docker run --rm -v ~/mbi_lab1:/tmp -w /tmp wkusmirek/dnaasm dnaasm -assembly -k 55 -genome_length 1643938 -insert_size_mean_inward 400 -insert_size_std_dev_inward 20 -single_edge_counter_threshold 5 -i1_1 Sim_100_400_1.fq -i1_2 Sim_100_400_2.fq ``` Część logów z `dnaasm/dnaasm_calc_0.log`: ``` [2022-Apr-04 20:26:56.954641] [info] - num of sequences: 328 [2022-Apr-04 20:26:56.954700] [info] - sum: 3284345 [2022-Apr-04 20:26:56.954712] [info] - max: 113244 [2022-Apr-04 20:26:56.954725] [info] - average: 10013.247070 [2022-Apr-04 20:26:56.954734] [info] - median: 165.000000 [2022-Apr-04 20:26:56.954757] [info] - N50: 48593 ``` ### Q&A: **Czy suma długości wygenerowanych sekwencji jest w przybliżeniu równa długości badanego genomu? Dlaczego?** $$ 3284345 / 1643938 \approx 1.997 $$ A zatem suma długości wygenerowanych sekwencji stanowi niespełna dwukrotność długości badanego genomu. Jest to oczekiwany wynik, gdyż assemblowane są 2 nici DNA, a wartość jest równa około 2. Jest ona mniejsza, gdyż czasem idealne dopasowanie nie jest możliwe. **Czy plik z sekwencjami wynikowymi jest w formacie FASTA czy FASTQ?** FASTA, ponieważ zawiera tylko sekwencje i ich identyfikatory (które są poprzedzone znakiem `>`). **Czy możliwa jest konwersja pliku w formacie FASTA na FASTQ? Jaką informacje należy wówczas sztucznie wygenerować?** Tak, należy wówczas sztucznie wygenerować wartości jakości (*quality values*) dla każdej sekwencji. **Czy możliwa jest konwersja pliku w formacie FASTQ na FASTA? Jaka informacja przy takiej konwersji zostaje utracona?** Tak, zostanie wtedy utracona informacja o wartościach jakości (*quality values*) sekwencji. ## Sprawdzenie wyników ```console $ docker run --rm -v ~/mbi_lab1:/tmp -w /tmp wkusmirek/quast quast.py -R GCF_001506805.1_ASM150680v1_genomic.fna out ``` Zawartość pliku `quast_results/latest/report.txt`: ``` All statistics are based on contigs of size >= 500 bp, unless otherwise noted (e.g., "# contigs (>= 0 bp)" and "Total length (>= 0 bp)" include all contigs). Assembly out # contigs (>= 0 bp) 328 # contigs (>= 1000 bp) 110 # contigs (>= 5000 bp) 94 # contigs (>= 10000 bp) 76 # contigs (>= 25000 bp) 46 # contigs (>= 50000 bp) 20 Total length (>= 0 bp) 3284345 Total length (>= 1000 bp) 3250691 Total length (>= 5000 bp) 3199118 Total length (>= 10000 bp) 3060380 Total length (>= 25000 bp) 2557982 Total length (>= 50000 bp) 1580324 # contigs 120 Largest contig 113244 Total length 3257087 Reference length 1643938 GC (%) 30.29 Reference GC (%) 30.29 N50 48593 NG50 92440 N75 26170 NG75 59133 L50 21 LG50 8 L75 42 LG75 14 # misassemblies 0 # misassembled contigs 0 Misassembled contigs length 0 # local misassemblies 0 # unaligned mis. contigs 0 # unaligned contigs 0 + 0 part Unaligned length 0 Genome fraction (%) 99.150 Duplication ratio 1.998 # N's per 100 kbp 0.00 # mismatches per 100 kbp 0.12 # indels per 100 kbp 0.00 Largest alignment 113244 Total aligned length 3257087 NA50 48593 NGA50 92440 NA75 26170 NGA75 59133 LA50 21 LGA50 8 LA75 42 LGA75 14 ``` ### Q&A: **Czy asemblacja de novo genomu pozwoliła uzyskać satysfakcjonujące wyniki? Dlaczego?** Asemblacja de novo pozwoliła uzyskać bardzo satysfakcjonujące wyniki (genom został zrekonstruowany w 99.150%), ponieważ pokrycie genomu odczytami było 50-krotne, a same odczyty nie były zbyt długie w porównaniu do całkowitej długości genomu. Ponadto poziom wprowadzonych błędów też jest dosyć niski. **Czym są translokacje w genomie? Czy aplikacja QUAST dostarcza informacji o liczbie translokacji w wynikach asemblacji de novo względem genomu referencyjnego? Jeśli tak, to skąd mogą być takie informacje odczytane?** Translokacje to błędy polegające na tym, że fragment sekwencji znajdzie się na innym chromosomie niż w genomie referencyjnym. Aplikacja QUAST dostarcza informacji na ten temat i można je odczytać z pliku `misassemblies_report.txt` znajdującym się w katalogu `contigs_reports`: ``` All statistics are based on contigs of size >= 500 bp, unless otherwise noted (e.g., "# contigs (>= 0 bp)" and "Total length (>= 0 bp)" include all contigs). Assembly out # misassemblies 167 # relocations 167 # translocations 0 # inversions 0 # misassembled contigs 3 Misassembled contigs length 399526 # local misassemblies 101 # unaligned mis. contigs 1 # mismatches 3 # indels 0 # indels (<= 5 bp) 0 # indels (> 5 bp) 0 Indels length 0 ``` ## Zadanie implementacyjne Wskaźnik GC mówi o stosunku liczby wystąpień nukleotydów G lub C do całej długości sekwencji genomu. Według aplikacji QUAST wartość ta wynosiła 30.29% (`reference GC`). Kod źródłowy skryptu: ```python from Bio import SeqIO def main(): gc_count = 0 total = 0 for seq_record in SeqIO.parse('GCF_001506805.1_ASM150680v1_genomic.fna', 'fasta'): seq = seq_record.seq total += len(seq_record) for letter in seq: if letter == 'G' or letter == 'C': gc_count += 1 print(f'GC count: {gc_count}') print(f'Total length: {total}') print(f'Percent: {100.0 * gc_count / total}%') if __name__ == '__main__': main() ``` Wynik działania skryptu: ``` GC count: 498007 Total length: 1643938 Percent: 30.29353905074279% ``` Wynik zgadza się z wartością zwróconą przez QUAST.