# [MBI] Ćwiczenie 3: Resekwencjonowanie genomu człowieka

**Data: 15.05.2022**

## Skład zespołu

* Rafał Kulus (300249)
* Kamil Przybyła (300254)

## Mapowanie

### Q&A:

**Sprawdź zawartość wygenerowango pliku SAM. Jaka jest typowa
długość odczytów?**

Odczyty są długości wynoszącej około 76 symboli.

**Jak jest różnica w wielkości plików FASTQ, BAM, SAM?**

| Plik | Rozmiar |
| -------- | -------- |
| FASTQ | 57MB |
| BAM | 14.1MB |
| SAM | 75.9MB |

**Uwaga:** Wielkości plików podane zostały w megabajtach ($10^6$ bajtów), a nie mebibajtach ($2^20$ bajtów).

## Wizualizacja zawartości pliku BAM w programie IGV

### Q&A:

**Jaka jest pozycja tego wariantu? Ile odczytów wskazuje na wariant a ile na referencje? Czy jest to wariant homo czy heterozygotyczny?**

Pozycja wariantu: chr1:156.518.379

Na wariant wskazuje 230 odczytów.

Na referencję wskazuje 1 odczyt.

Jest to wariant homozygotyczny, ponieważ na danej pozycji znajdują się identyczne allele.


## Wykrywanie wariantów

### Q&A:

**Ile wariantów zawiera plik?**

Plik zawiera 6050 wariantów.

**Ile wariantów zostało po filtracji? Jakich innych parametrów możemy użyć do dalszej filtracji liczby wariantów?**

Po filtracji zostało 241 wariantów.

Do dalszej filtracji można wykorzystać dowolny z paramerów wypisanych z listy, tworząc wyrażenia filtrujące, które spełniać będą musiały interesujące nas warianty:

##INFO=<ID=INDEL,Number=0,Type=Flag,Description="Indicates that the variant is an INDEL.">
##INFO=<ID=IDV,Number=1,Type=Integer,Description="Maximum number of reads supporting an indel">
##INFO=<ID=IMF,Number=1,Type=Float,Description="Maximum fraction of reads supporting an indel">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Raw read depth">
##INFO=<ID=VDB,Number=1,Type=Float,Description="Variant Distance Bias for filtering splice-site artefacts in RNA-seq data (bigger is better)",Version="3">
##INFO=<ID=RPB,Number=1,Type=Float,Description="Mann-Whitney U test of Read Position Bias (bigger is better)">
##INFO=<ID=MQB,Number=1,Type=Float,Description="Mann-Whitney U test of Mapping Quality Bias (bigger is better)">
##INFO=<ID=BQB,Number=1,Type=Float,Description="Mann-Whitney U test of Base Quality Bias (bigger is better)">
##INFO=<ID=MQSB,Number=1,Type=Float,Description="Mann-Whitney U test of Mapping Quality vs Strand Bias (bigger is better)">
##INFO=<ID=SGB,Number=1,Type=Float,Description="Segregation based metric.">
##INFO=<ID=MQ0F,Number=1,Type=Float,Description="Fraction of MQ0 reads (smaller is better)">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="List of Phred-scaled genotype likelihoods">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##INFO=<ID=ICB,Number=1,Type=Float,Description="Inbreeding Coefficient Binomial test (bigger is better)">
##INFO=<ID=HOB,Number=1,Type=Float,Description="Bias in the number of HOMs number (smaller is better)">
##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes for each ALT allele, in the same order as listed">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
##INFO=<ID=DP4,Number=4,Type=Integer,Description="Number of high-quality ref-forward , ref-reverse, alt-forward and alt-reverse bases">
##INFO=<ID=MQ,Number=1,Type=Integer,Description="Average mapping quality">

Jak widać filtrować możemy na przykład po jakości mapowania, wynikach różnych testów statystycznych czy liczbie alleli.

## Adnotacje wariantów

### Q&A:

**Jaki typ wariantu przeważa?**

Przeważa wariant typu *intron_variant* (35%), czyli wariant występujący w intronie (transcript variant occurring within an intron).


**Załącz do sprawozdania wiersze odpowiadające temu wariantowi. Czy
jest to wariant w części kodującej?**

.	1:156518379-156518379	A	missense_variant	MODERATE	IQGAP3	ENSG00000183856	Transcript	ENST00000361170.2	protein_coding	17/38	-	-	-	1998	1987	663	R/C	Cgt/Tgt	rs744224,COSV63249352	-	-1	-HGNC	20669	-	-	-	-	tolerated(0.13)	benign(0)	0.4071	-	0,1	0,1	-	-	-	-	-	-

.	1:156518379-156518379	A	missense_variant,NMD_transcript_variant	MODERATE	IQGAP3	ENSG00000183856	Transcript	ENST00000491900.1	nonsense_mediated_decay	16/38	-	-	-	1974	1858	620	R/C	Cgt/Tgt	rs744224,COSV63249352	-	-1	-	HGNC	20669	-	-	-	-	tolerated(0.13)	benign(0)	0.4071	-	0,1	0,1	-	-	-	-	-	-

Jest to wariant w części kodującej, o czym świadczy istnienie wartości dla pól `Exon`, `cDNA position`, `CDS position`, `Protein position`, `Amino acids` oraz `Codons`.

## Zadanie implementacyjne

import pandas as pd
import pyranges as pr


def main():
    header_names = ['geneName', 'name', 'chrom', 'strand', 'txStart', 'txEnd',
                    'cdsStart', 'cdsEnd', 'exonCount', 'exonStarts', 'exonEnds']
    df = pd.read_csv('refFlat.txt', header=None, names=header_names, sep="\t")
    df = df.rename(columns={'chrom': 'Chromosome', 'txStart': 'Start', 'txEnd': 'End'})
    prdf = pr.PyRanges(df)

    vcf_header_names = ['CHROM', 'POS', 'ID', 'REF', 'ALT', 'QUAL',
                        'FILTER', 'INFO', 'FORMAT', 'coriell_chr1.bam']
    df_vcf = pd.read_csv('coriell_chr1.vcf', header=None,
                         names=vcf_header_names, sep="\t", comment='#')
    df_vcf = df_vcf.rename(columns={'CHROM': 'Chromosome', 'POS': 'Start'})
    df_vcf['End'] = df_vcf['Start']
    prdf_vcf = pr.PyRanges(df_vcf)

    # We have only 1. chromosome
    prdf_chr1 = prdf[prdf.Chromosome == 'chr1']

    genes = prdf_chr1.geneName.unique()
    counts = len(genes) * [0]

    for idx, gene in enumerate(genes):
        prdf_gene = prdf_chr1[prdf_chr1.geneName == gene]
        intersections = prdf_vcf.intersect(prdf_gene, how='containment')
        if(len(intersections)):
            counts[idx] = len(intersections.Start.unique())

    result = pd.DataFrame({'Gene': genes, 'Count': counts})
    result.to_csv('result.csv', header=None, index=None, sep='\t')
    print(result)


if __name__ == '__main__':
    main()

3669

3290

MIR6859-1	0
MIR6859-2	0
MIR6859-3	0
MIR6859-4	0
MIR1302-2	0
MIR1302-9	0
MIR1302-10	0
MIR1302-11	0
FAM138A	3
FAM138F	3
FAM138C	3
LOC729737	0
LOC100132287	0
LOC100132062	0
LOC101928626	0
FAM87B	0
LINC00115	0
FAM41C	0
LINC02593	1
CENPS	1
CD5L	1
LOC100288175	0
MIGA1	2
LINC01342	0
(...)

W sumie w poszczególnych genach znajduje się 3669 wariantów (z 6050 wariantów dostępnych w pliku `coriell_chr1.vcf`), ale tylko 3290 unikalnych (wartość wyliczona inną wersją skryptu, "na piechotę"). Oznacza to zatem, że współrzędne (Start, End) niektórych genów na siebie nachodzą.

W `counts[idx] = len(intersections.Start.unique())` są dodawane tylko unikalne warianty z tego względu, że `geneName` nie jest wartością unikalną i w efekcie jeden wariant potrafi zostać zaliczony do kilku przedziałów o tym samym `geneName`. Gdyby zliczać także duplikaty, suma wystąpień wariantów we wszystkich genach wyniosłaby ponad 11 tysięcy, co przekracza prawie dwukrotnie liczbę dostępnych wariantów w pliku `coriell_chr1.vcf`.