CB12
Z MBI
Obsah
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
- 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 count | 10,350,729 |
| item bases | 162,179,256 (5.32%) |
| item total | 162,179,256 (5.32%) |
| smallest item | 1 |
| average item | 16 |
| biggest item | 3,732 |
| smallest score | 186 |
| average score | 333 |
| biggest score | 1,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:
- Selecton, clanok
- Data monkey clanok
- Skusili sme na Selecton poslat CLCA4 zo 7 cicavcov, subor tu: [8]
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
