1-BIN-301, 2-AIN-501 Methods in Bioinformatics, 2021/22

Introduction · Rules · Tasks and dates · Materials · Moodle · Discussion
Cvičenia vo štvrtok o 14:00 sú určené pre študentov BIN, INF, mINF, mAIN, DAV. Cvičenia vo štvrtok o 17:20 sú pre študentov z PriFUK a z fyzikálnych odborov. Obidvoje cvičenia sa budú konať už v prvom týždni semestra.


CB12: Rozdiel medzi revíziami

Z MBI
Prejsť na: navigácia, hľadanie
(Vytvorená stránka „=CB12= ==Zhrnutie semestra== * vid prezentacia k cviceniu ==Uvod do teorie grafov== * vid prezentacia k cviceniu ==Ukážka práce v Linuxe== ===Prvá časť - prípr...“)
 
(CB12)
 
Riadok 1: Riadok 1:
=CB12=
+
==Zhrnutie semestra a úvod do teórie grafov==
==Zhrnutie semestra==
+
* vid prezentacie k cviceniu
* vid prezentacia k cviceniu
+
  
==Uvod do teorie grafov==
+
==Uniprot==
* vid prezentacia k cviceniu
+
* Prehladnejsi pohlad na proteiny, vela linkov na ine databazy, cast vytvarana rucne
 +
** Pozrieme sa na známy koronavírusový proteín Spike
 +
** Nájdime ho na stránke http://www.uniprot.org/ pod názvom SPIKE_SARS2
 +
** Pozrime si podrobne jeho stránku, ktoré časti boli predpovedané bioinformatickými metódami z prednášky?
 +
** Všimnime si niektorú Pfam doménu a pozrime si jej stránku
 +
 
 +
==PSI BLAST==
 +
* Toto cvičenie je z časti inšpirované stránkou [http://www.ncbi.nlm.nih.gov/Class/FieldGuide/problem_set.html]
 +
* Budeme uvažovať vzdialene podobné enzýmy
 +
** Bis(5'-adenosyl)-triphosphatase ([http://www.uniprot.org/uniprot/P49789 Uniprot])
 +
** Galactose-1-phosphate uridylyltransferase (GALT/GAL7) ([http://www.uniprot.org/uniprot/P31764 Uniprot])
 +
** Ich domény patria v databáze Pfam do toho istého klanu
 +
* Skúsme nájsť túto podobnosť v BLASTe: http://blast.ncbi.nlm.nih.gov/ v časti proteíny, zvoľme databázu Swissport, ako Query zadajme Accesion nášho proteínu P49789, spustime program PSI-BLAST
 +
* V prvom kole PSI-BLAST spúšťa bežný BLASTP
 +
* GAL gén (konkrétne GAL7_HAEIN, accession P31764) sa nachádza medzi výsledkami, ale má príliš vysokú E-value
 +
* Spustíme teraz druhú iteráciu PSI-BLAST, ktorá zostaví profil z proteínov s nízkou E-value v prvej iterácii
 +
* Aká je E-value nájdeného zarovnania?
 +
* Ak by výpočet dlho trval, výsledky sú tu: [http://compbio.fmph.uniba.sk/vyuka/mbi-data/cb08/psi-blast1.html 1. kolo], [http://compbio.fmph.uniba.sk/vyuka/mbi-data/cb08/psi-blast2.html 2. kolo]
 +
 
 +
==Sekvenovanie v UCSC genome browseri==
 +
* Vráťte sa na UCSC genome browser http://genome-euro.ucsc.edu/
 +
* Pozrieme si niekoľko vecí týkajúcich sa sekvenovania a skladania genómov
 +
* Hore v modrom menu zvoľte Genomes, časť Other
 +
* Na ďalšej stránke zvoľte človeka a pomocou menu Human Assembly '''zistite, kedy boli pridané posledné dve verzie ľudského genómu (hg19 a hg38)'''
 +
* Na tej istej stránke dole nájdete stručný popis zvolenej verzie genómu. '''Pre ktoré oblasti genómu máme v hg38 najviac alternatívnych verzií? (haplotypes)'''
 +
* Prejdite na región chr21:31,250,000-31,300,000 v hg19  touto linkou: [http://genome-euro.ucsc.edu/cgi-bin/hgTracks?db=hg19&position=chr21%3A31250000-31300000]
 +
* Zapnite si tracky Mapability a RepeatMasker na "full"
 +
* Mapability: nakoľko sa daný úsek opakuje v genóme a či teda vieme jednoznačne jeho čítania namapovať pri použití Next generation sequencing
 +
* Ako a prečo sa  pri rôznych dĺžkach čítaní líšia? (Keď kliknete na linku "Mapability", môžete si prečítať bližšie detaily.)
 +
* Približne v strede zobrazeného regiónu je pokles mapovateľnosti. '''Akému typu opakovania zodpovedá?''' (pozrite track RepeatMasker)
 +
* Zapnite si tracky "Assembly" a "Gaps" a pozrite si región chr2:110,000,000-110,300,000 v hg19: [http://genome-euro.ucsc.edu/cgi-bin/hgTracks?db=hg19&position=chr2%3A110000000-110300000] '''Aká dlhá je neosekvenovaná medzera (gap) v strede tohto regiónu?''' Približnú veľkosť môžete odčítať z obrázku, presnejší údaj zistíte kliknutím na čierny obdĺžnik zodpovedajúci tejto medzere (úplne presná dĺžka aj tak nebola známa, nakoľko nebola osekvenovaná).
 +
* Cez menu položku View, In other genomes si pozrite, ako zobrazený úsek vyzerá vo verzii hg38. Ako sa zmenila dĺžka z pôvodných 300kb?
 +
* Prejdite na genóm Rhesus, verzia rheMac2, región chr7:59,022,000-59,024,000 [http://genome-euro.ucsc.edu/cgi-bin/hgTracks?db=rheMac2&position=chr7%3A59022000-59024000], zapnite si tracky Contigs, Gaps, Quality scores
 +
* '''Aké typy problémov v kvalite sekvencie v tomto regióne vidíte?'''
 +
** Opäť si môžete pozrieť, či sa problémy odstránili a ako sa zmenila dĺžka sekvencie v najnovšej verzii rheMac8
 +
 
 +
==Gény, evolúcia a komparatívna genomika v UCSC genome browseri (cvičenie pri počítači)==
 +
 
 +
* Zobrazme si gén CLCA4 [http://genome-euro.ucsc.edu/cgi-bin/hgTracks?db=hg38&position=chr1%3A86538658-86589173]
 +
* Zapnite si štandardnú sadu track-ov
 +
* Po kliknutí na gén si môžete prečítať o jeho funkcii, po kliknutí na ľavú lištu alebo na názov tracku v zozname na spodku stránky si môžete prečítať viac o tracku a meniť nastavenia
 +
* V tracku RefSeq genes si všimnite, že v tejto databáze má tento gén dve formy zostrihu, jedna z nich sa považuje za nekódujúcu, pretína sa aj s necharakterizovanou nekódujúcou RNA na opačnom vlákne
 +
** Track RefSeq a jeho subtrack RefSeq Curated treba zapnut na pack
 +
* Nižšie vidíte track H3K27Ac Mark (Often Found Near Regulatory Elements) on 7 cell lines from ENCODE, kde bola táto histónová modifikácia v okolí génu detegovaná?
 +
* Všimnite si aj track DNase I Hypersensitivity, ktorý zobrazuje otvorený chromatin, prístupný pre viazanie transkripčných faktorov. Všimnite si jeho súvis s H3K27Ac trackom
 +
* Obidva tracky sú súčasťou tracku ENCODE regulation, v ktorom si môžete zapnúť aj ďalšie pod-tracky
 +
 
 +
* Vsimnime si track Vertebrate Multiz Alignment & Conservation (100 Species)
 +
** v spodnej casti tracku vidime zarovnania s roznymi inymi genomami
 +
** v nastaveniach tracku zapnite Element Conservation (phastCons) na full a Conserved Elements na dense
 +
** v tomto tracku vidíme PhyloP, co zobrazuje uroven konzerovanosti danej bazy len na zaklade jedneho stlpca zarovnania a dva vysledky z phyloHMM phastCons, ktory berie do uvahy aj okolite stlpce
 +
* Konkretne cast Conserved elements zobrazuje konkretne useky, ktore su najvac konzervovane
 +
** Ak chceme zistit, kolko percent genomu tieto useky pokryvaju, ideme na modrej liste do casti Tools->Table browser, zvolime group Comparative genomics, track Conservation, table 100 Vert. El, region zvolime genome (v celom genome)  a stlacime tlacidlo Summary/statistics, dostaneme nieco taketo:
 +
<TABLE border=1>
 +
<TR><TD>item count</TD><TD ALIGN=RIGHT>10,350,729</TD></TR>
 +
<TR><TD>item bases</TD><TD ALIGN=RIGHT>162,179,256 (5.32%)</TD></TR>
 +
<TR><TD>item total</TD><TD ALIGN=RIGHT>162,179,256 (5.32%)</TD></TR>
 +
<TR><TD>smallest item</TD><TD ALIGN=RIGHT>1</TD></TR>
 +
<TR><TD>average item</TD><TD ALIGN=RIGHT>16</TD></TR>
 +
<TR><TD>biggest item</TD><TD ALIGN=RIGHT>3,732</TD></TR>
 +
<TR><TD>smallest score</TD><TD ALIGN=RIGHT>186</TD></TR>
 +
<TR><TD>average score</TD><TD ALIGN=RIGHT>333</TD></TR>
 +
<TR><TD>biggest score</TD><TD ALIGN=RIGHT>1,000</TD></TR>
 +
</TABLE>
 +
** Ak by nas zaujimali iba velmi dlhe "conserved elements", v Table browser stlacime tlacidlo Filter a na dalsej obrazovke do policka Free-form query dame '''chromEnd-chromStart>=1500'''
 +
** Potom mozeme skusit Summary/Statistics alebo vystup typu Hyperlinks to genome browser a Get output - dostaneme zoznam tychto elementov a kazdy si mozeme jednym klikom pozriet v browseri, napr. taketo
 +
*** [http://genome-euro.ucsc.edu/cgi-bin/hgTracks?db=hg38&position=chr1:50201403-50203312 lod=24051 at chr1:50201403-50203312]
 +
*** [http://genome-euro.ucsc.edu/cgi-bin/hgTracks?db=hg38&position=chr1:55663689-55667047 lod=1899 at chr1:55663689-55667047] atd
 +
 
 +
* Pozrime si teraz ten isty gen CLCA4 v starsej verzii genomu hg18 [http://genome-euro.ucsc.edu/cgi-bin/hgTracks?db=hg18&position=chr1%3A86776929-86827444]
 +
** V casti Genes and Gene Prediction Tracks zapnite track Pos Sel Genes, ktory obsahuje geny s '''pozitivnym vyberom''' (cervenou, pripadne slabsie fialovou a modrou)
 +
** Ked kliknete na cerveny obdlznik pre tento gen, uvidite, v ktorych castiach fylogenetickeho stromu bol detegovany pozitivny vyber
 +
** Po priblizeni do jedneho z exonov [http://genome-euro.ucsc.edu/cgi-bin/hgTracks?db=hg18&position=chr1%3A86805823-86805917] vidite dosledky nesynonymnych mutacii
 +
 
 +
Poznamka: Existuju aj webservery na predikciu pozitivneho vyberu, napriklad tieto dva:
 +
* [http://selecton.tau.ac.il/ Selecton], [http://www.tau.ac.il/~talp/publications/selecton2007.pdf clanok]
 +
* [http://www.datamonkey.org/ Data monkey] [http://mbe.oxfordjournals.org/cgi/content/abstract/22/5/1208 clanok]
 +
* Skusili sme na Selecton poslat CLCA4 zo 7 cicavcov, subor tu: [http://compbio.fmph.uniba.sk/vyuka/mbi-data/cb07/clca4.mfa]
 +
** vysledky [http://compbio.fmph.uniba.sk/vyuka/mbi-data/cb07/clca4-selecton.html] a [http://compbio.fmph.uniba.sk/vyuka/mbi-data/cb07/clca4-omega.txt] (metoda ale odporuca aspon 10 homologov)
  
 
==Ukážka práce v Linuxe==
 
==Ukážka práce v Linuxe==

Aktuálna revízia z 14:42, 10. december 2020

Zhrnutie semestra a úvod do teórie grafov

  • vid prezentacie k cviceniu

Uniprot

  • Prehladnejsi pohlad na proteiny, vela linkov na ine databazy, cast vytvarana rucne
    • Pozrieme sa na známy koronavírusový proteín Spike
    • Nájdime ho na stránke http://www.uniprot.org/ pod názvom SPIKE_SARS2
    • Pozrime si podrobne jeho stránku, ktoré časti boli predpovedané bioinformatickými metódami z prednášky?
    • Všimnime si niektorú Pfam doménu a pozrime si jej stránku

PSI BLAST

  • Toto cvičenie je z časti inšpirované stránkou [1]
  • Budeme uvažovať vzdialene podobné enzýmy
    • Bis(5'-adenosyl)-triphosphatase (Uniprot)
    • Galactose-1-phosphate uridylyltransferase (GALT/GAL7) (Uniprot)
    • Ich domény patria v databáze Pfam do toho istého klanu
  • Skúsme nájsť túto podobnosť v BLASTe: http://blast.ncbi.nlm.nih.gov/ v časti proteíny, zvoľme databázu Swissport, ako Query zadajme Accesion nášho proteínu P49789, spustime program PSI-BLAST
  • V prvom kole PSI-BLAST spúšťa bežný BLASTP
  • GAL gén (konkrétne GAL7_HAEIN, accession P31764) sa nachádza medzi výsledkami, ale má príliš vysokú E-value
  • Spustíme teraz druhú iteráciu PSI-BLAST, ktorá zostaví profil z proteínov s nízkou E-value v prvej iterácii
  • Aká je E-value nájdeného zarovnania?
  • Ak by výpočet dlho trval, výsledky sú tu: 1. kolo, 2. kolo

Sekvenovanie v UCSC genome browseri

  • Vráťte sa na UCSC genome browser http://genome-euro.ucsc.edu/
  • Pozrieme si niekoľko vecí týkajúcich sa sekvenovania a skladania genómov
  • Hore v modrom menu zvoľte Genomes, časť Other
  • Na ďalšej stránke zvoľte človeka a pomocou menu Human Assembly zistite, kedy boli pridané posledné dve verzie ľudského genómu (hg19 a hg38)
  • Na tej istej stránke dole nájdete stručný popis zvolenej verzie genómu. Pre ktoré oblasti genómu máme v hg38 najviac alternatívnych verzií? (haplotypes)
  • Prejdite na región chr21:31,250,000-31,300,000 v hg19 touto linkou: [2]
  • Zapnite si tracky Mapability a RepeatMasker na "full"
  • Mapability: nakoľko sa daný úsek opakuje v genóme a či teda vieme jednoznačne jeho čítania namapovať pri použití Next generation sequencing
  • Ako a prečo sa pri rôznych dĺžkach čítaní líšia? (Keď kliknete na linku "Mapability", môžete si prečítať bližšie detaily.)
  • Približne v strede zobrazeného regiónu je pokles mapovateľnosti. Akému typu opakovania zodpovedá? (pozrite track RepeatMasker)
  • Zapnite si tracky "Assembly" a "Gaps" a pozrite si región chr2:110,000,000-110,300,000 v hg19: [3] Aká dlhá je neosekvenovaná medzera (gap) v strede tohto regiónu? Približnú veľkosť môžete odčítať z obrázku, presnejší údaj zistíte kliknutím na čierny obdĺžnik zodpovedajúci tejto medzere (úplne presná dĺžka aj tak nebola známa, nakoľko nebola osekvenovaná).
  • Cez menu položku View, In other genomes si pozrite, ako zobrazený úsek vyzerá vo verzii hg38. Ako sa zmenila dĺžka z pôvodných 300kb?
  • Prejdite na genóm Rhesus, verzia rheMac2, región chr7:59,022,000-59,024,000 [4], zapnite si tracky Contigs, Gaps, Quality scores
  • Aké typy problémov v kvalite sekvencie v tomto regióne vidíte?
    • Opäť si môžete pozrieť, či sa problémy odstránili a ako sa zmenila dĺžka sekvencie v najnovšej verzii rheMac8

Gény, evolúcia a komparatívna genomika v UCSC genome browseri (cvičenie pri počítači)

  • Zobrazme si gén CLCA4 [5]
  • Zapnite si štandardnú sadu track-ov
  • Po kliknutí na gén si môžete prečítať o jeho funkcii, po kliknutí na ľavú lištu alebo na názov tracku v zozname na spodku stránky si môžete prečítať viac o tracku a meniť nastavenia
  • V tracku RefSeq genes si všimnite, že v tejto databáze má tento gén dve formy zostrihu, jedna z nich sa považuje za nekódujúcu, pretína sa aj s necharakterizovanou nekódujúcou RNA na opačnom vlákne
    • Track RefSeq a jeho subtrack RefSeq Curated treba zapnut na pack
  • Nižšie vidíte track H3K27Ac Mark (Often Found Near Regulatory Elements) on 7 cell lines from ENCODE, kde bola táto histónová modifikácia v okolí génu detegovaná?
  • Všimnite si aj track DNase I Hypersensitivity, ktorý zobrazuje otvorený chromatin, prístupný pre viazanie transkripčných faktorov. Všimnite si jeho súvis s H3K27Ac trackom
  • Obidva tracky sú súčasťou tracku ENCODE regulation, v ktorom si môžete zapnúť aj ďalšie pod-tracky
  • Vsimnime si track Vertebrate Multiz Alignment & Conservation (100 Species)
    • v spodnej casti tracku vidime zarovnania s roznymi inymi genomami
    • v nastaveniach tracku zapnite Element Conservation (phastCons) na full a Conserved Elements na dense
    • v tomto tracku vidíme PhyloP, co zobrazuje uroven konzerovanosti danej bazy len na zaklade jedneho stlpca zarovnania a dva vysledky z phyloHMM phastCons, ktory berie do uvahy aj okolite stlpce
  • Konkretne cast Conserved elements zobrazuje konkretne useky, ktore su najvac konzervovane
    • Ak chceme zistit, kolko percent genomu tieto useky pokryvaju, ideme na modrej liste do casti Tools->Table browser, zvolime group Comparative genomics, track Conservation, table 100 Vert. El, region zvolime genome (v celom genome) a stlacime tlacidlo Summary/statistics, dostaneme nieco taketo:
item count10,350,729
item bases162,179,256 (5.32%)
item total162,179,256 (5.32%)
smallest item1
average item16
biggest item3,732
smallest score186
average score333
biggest score1,000
    • Ak by nas zaujimali iba velmi dlhe "conserved elements", v Table browser stlacime tlacidlo Filter a na dalsej obrazovke do policka Free-form query dame chromEnd-chromStart>=1500
    • Potom mozeme skusit Summary/Statistics alebo vystup typu Hyperlinks to genome browser a Get output - dostaneme zoznam tychto elementov a kazdy si mozeme jednym klikom pozriet v browseri, napr. taketo
  • Pozrime si teraz ten isty gen CLCA4 v starsej verzii genomu hg18 [6]
    • V casti Genes and Gene Prediction Tracks zapnite track Pos Sel Genes, ktory obsahuje geny s pozitivnym vyberom (cervenou, pripadne slabsie fialovou a modrou)
    • Ked kliknete na cerveny obdlznik pre tento gen, uvidite, v ktorych castiach fylogenetickeho stromu bol detegovany pozitivny vyber
    • Po priblizeni do jedneho z exonov [7] vidite dosledky nesynonymnych mutacii

Poznamka: Existuju aj webservery na predikciu pozitivneho vyberu, napriklad tieto dva:

Ukážka práce v Linuxe

Prvá časť - príprava

  • Prihláste sa na server podľa pokynov.
  • Potom spúšťajte jednotlivé príkazy podľa pokynov nižšie.
  • Odporúčame príkazy kopírovať myšou (v internetovom prehliadači vysvietiť, stlačiť Ctrl-C, v konzole Ctrl-Shift-V)


# riadky začínajúce mrežou # sú komentáre, netreba ich spúšťať

# Dôležité: v príkazoch nižšie xx nahraďte vašim číslom skupiny, napr. 01
mkdir xx
cd xx
# príkaz mkdir (make directory) vytvoril priečinok
# príkaz cd (change directory) zmenil váš aktuálny priečinok na tento nový

# v konzole by ste mali mať user@server:~/xx$
# kde xx je číslo vašej skupiny, napr. 01

# stiahneme si súbor s dátami zo stránky
wget http://compbio.fmph.uniba.sk/vyuka/mbi-data/cb12.zip
# rozzipujeme ho
unzip cb12.zip

Druhá časť - skladanie genómov, mapovanie čítaní, zarovnanie

# prejdeme na priečinok s prvou časťou ohľadom sekvenovania
cd 1-seq

# ls vypíše zoznam súborov v priečinku
ls
# ls -l vypíše dlhšiu informáciu (long)
ls -l
# ls -lSh usporiada súbory podľa veľkosti (Size) a veľkosti vypíše priateľskejšie pre ľudí (human)
ls -lSh

# mali by sme vidieť kúsok sekvencie z E.coli (prípona .fasta)
# a 2 súbory zo sekvenovania prístrojom Illumina Miseq  (prípona .fastaq.gz)
# tieto súbory obsahujú čítania z vyššie uvedeného kúsku genómu


# ideme skladať genóm, bude to trvať dlho, preto to chceme spustiť na pozadí
# aby sme mohli medzitým robiť niečo iné
screen # stlačte Enter
# spustite skladanie programom spades
spades.py -t 1 -m 1 --pe1-1 miseq_R1.fastq.gz --pe1-2 miseq_R2.fastq.gz -o spades > spades.log
# stlačte naraz Ctrl-a potom d
# spades teraz beží na pozadí

# príkaz top zobrazí bežiace procesy
# ukončíte ho stlačením q (quit)
top

# príkaz less umožňuje prezerať si obsah textového súboru
# aj príkaz less ukončíte stlačením q, šípkami sa pohybujete po súbore
less ref.fasta
# čítania sú komprimované, preto namiesto less použijeme zless
zless miseq_R1.fastq.gz
# tieto príkazy spočítajú počet riadkov - ako z toho zistíme počet čítaní?
zcat miseq_R1.fastq.gz | wc -l 
zcat miseq_R2.fastq.gz | wc -l 

# keď spades skončí, vrátime sa do screen a ukončíme ho
screen -r
# exit ukončí screen
exit

# spades dal výstup do podpriečinku spades, pozrime si ho
ls spades
# skopírujeme si hlavný výsledok do nášho priečinka (cp = copy)
cp -ip spades/contigs.fasta spades.fasta
less spades.fasta
# pozrime si hlavičky jednotlivých sekvencií vo fasta súbore
grep '>' spades.fasta

# programom last si spravíme dotplot referencia vs. naše skladanie
# 1) vytvorenie indexu pre referenciu
lastdb ref.fasta ref.fasta 
# 2) samotné zarovnanie
lastal -f TAB ref.fasta spades.fasta > aln.tab
# 3) vytvorenie obrázku s dotplotom
last-dotplot aln.tab aln.png

# a ešte dotplot referencia vs. referencia
# 2) samotné zarovnanie (index už máme)
lastal -f TAB ref.fasta ref.fasta > aln2.tab
# 3) vytvorenie obrázku s dotplotom
last-dotplot aln2.tab aln2.png

# pozrieme si dotploty programom eog
eog aln.png &
eog aln2.png &


# zarovnajme čítania k referenčnému genómu v 4 krokoch
# 1) indexovanie fasta súboru
bwa index ref.fasta
# 2) samotné zarovnávanie čítaní programom bwa
bwa mem ref.fasta miseq_R1.fastq.gz miseq_R2.fastq.gz > ref-miseq.sam
# 3) zmeníme textový sam formát na binárny bam formát
samtools view -S -b ref-miseq.sam | samtools sort - ref-miseq
# 4) vytvoríme index bam súboru
samtools index ref-miseq.bam

# pozrime sa na zoznam súborov od najnovšieho po najstarší
ls -lth
# sam súbor so zarovnaniami sa dá pozrieť, ale nie je veľmi prehľadný
less ref-miseq.sam


# vytvoríme aj zarovnanie nášho poskladaného genómu k referencii vo formáte bam
samtools faidx ref.fasta
lastal ref.fasta spades.fasta -E1e-20 | maf-convert sam > ref-spades.sam
samtools view -S -b -t ref.fasta.fai ref-spades.sam | samtools sort - ref-spades
samtools index ref-spades.bam

# výsledky si zobrazíme v grafickom prehliadači igv 
# obdoba genome browsera, ktorú si môžete nainštalovať na vašom počítači
# POZOR: POTREBUJE VEĽA PAMÄTE, SPUSTÍME IBA JEDEN NARAZ
igv -g ref.fasta
# pomocou Menu->File->Load from File otvorte ref-spades.bam a ref-miseq.bam
# pozrime si región ecoli-frag:224,000-244,000
#   Vidíte jednotlivé kontigy? Sedí tento pohľad s dotplotom? 
# a potom bližšie ecoli-frag:227,300-227,600
#   Všimnite si sekvenačné chyby rozdiely medzi referenciou a kontigmi

Tretia časť - hľadanie génov, RNA-seq

# v druhom cvičení si vyskúšame hľadanie génov
# najskôr sa presuňme do druhého priečinku
cd ../2-genes

# pozrime si, aké máme súbory
ls -lSh
# mali by sme mať kúsok referenčného genómu huby Aspergillus nidulans 
# fastq súbor s čítaniami z RNA-seq pre tento kúsok referencie
# gff súbor s anotáciou génov z databázy

# spustíme hľadač génov Augustus 2x:
# raz s parametrami priamo pre A.nidulans a raz s parametrami pre ľudský genóm
augustus --species=anidulans ref2.fasta > augustus-anidulans.gtf
augustus --species=human ref2.fasta > augustus-human.gtf

# RNA-seq zarovnáme k sekvencii nástrojom tophat2 (podporuje intróny)
bowtie2-build ref2.fasta ref2.fasta
tophat2 -i 10 -I 10000 --max-multihits 1 --output-dir rnaseq ref2.fasta rnaseq.fastq
samtools sort rnaseq/accepted_hits.bam rnaseq
samtools index rnaseq.bam

# predikcie génov a RNA-seq si pozrieme v igv
igv -g ref2.fasta
# v igv si otvorte annot.gff, augustus-anidulans.gtf, augustus-human.gtf, rnaseq.bam
# - ktoré parametre Augustusu dali presnejšie predpovede (za predpokladu, že anotácia je správna)
# - pozrite si zblízka niektorý gén s vysokou expresiou (napr. druhy gen sprava), 
#   mali by ste vidieť čítania podporujúce intróny