---
# System prepended metadata

title: '[MBI] Ćwiczenie 1: Asemblacja *de novo* DNA'

---

# [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.