Výpočetní metody v bioinformatice a systémové biologii

Týden 2 - Dalsi zdroje genomovych a anotacnich dat, uvod do programovani v R - 24. 2. 2020

Running IGV (or any other graphical application) on FI biolinux virtual computer:

1) Login (use your FI username and secondary password)
    $ ssh -X username@biolinux.fi.muni.cz
2) Run IGV (after login you will be in your FI home directory)
    $ igv
See also: gggenes, Gviz

Vzorove riesenie zalozene na odovzdanej ulohe:

1) Vizualnym prieskumom stopy RepeatMasker v UCSC Genome Browser, Human assembly hg38 najdeme napr. chr1:18,151,080-18,152,180 s L2 elementom prerusenym dvoma inymi typmi transpozonov (spoliehame sa na kontinuitu nazvu L2a v lavej aj pravej casti intervalu; ak by sme chceli mat istotu, ze nejde o dve nezavisle kopie L2a, museli by sme robit dalsie testy).

2) Pomocou menu View - DNA ziskame DNA sekvenciu oblasti:

>hg38_dna range=chr1:18151080-18152180 5'pad=0 3'pad=0 strand=+ repeatMasking=none
ACCTGAGCGTCCCAGGGCCTCTACCCAGGCTGTATCCTCTGTGTGGAGGG
ACCTTTCATCACTTCATTTAGGTCTCTGTTCAAATCTTCCCTTGTCAGAG
AGTCCTTCATGAGCCACTTGTCACAGTCAGGGTAATTAACACTAGCTGCT
GTAACAAACAATTCCAACCTCTTCTTGACACAGCACAGTAAAAGTGTATT
TTTTGTTCATGCAAAATCAAATGCATTTATTCCTGGTCAGAGGCTCCCCT
GAGCACTCACATCCAGGCAATGGTTCTGGGGGCCCAGCCTCTTTCATTTT
GTGCCTCTGTCATCCCTGGGAGGTTCCCCGCTGGATCCTCTGTATTCAGT
CTGCTGAAGAGTGAAGAGAAAGATCACATGGAGGCTCTGACCTAGAGTTT
CAAGGGCCAGGTCCAGAGTGGCATATACCATCGACATCCACATTCTAGTG
GCCAGAGCAGGTCACATGGCCAAATCCAACACTGATGGGGAGGGACAGTG
CATGCCACATATAACAGGCGGGGGAGAGGAGTGAATAATTCACACCATGG
TAGAGGATTTCACAGGATGTCTGGGGCCAGGCCTGGGCGACACTTATTCC
CTGCTTCCATGGACCACAACTGACTCACACATCCTCCTTAACTGCAGGGG
GCTGTGGGGGGCGGGTCTTGGAAATACGGAGAAGCTCATGGATATTGGTG
AAAACTGACAGTTTCTGCCCCACAACCTAGTAATTAAAGAGATCCTAGTA
GTGCTTCTCATCATGTGGGTGACCGTATAATTCATCTTCCAAATGGGACC
ATTTTAAGGTGAAAGGGGTCACTAATGATTATCCTCTGAACACTGGGTTT
CAAGCAGGTTGCTCCAGGGATTTGGGGCACGTGGTCACCCTACTCATCAC
TCTGCATCCTCTCAGCCTGCTTTATTGTTCTTTTTCACACTTAATCACTG
CCTAGGAGTGGACTCCGTGAGGACAGGGACCTTTCTCCCTCTGTTGCTAT
ATCCTCAGTGCCTAGAACTGGGCCTGGCACATTAGAGGTGCCCATTATAT
GCCTGGTAGATGAACGAATGAATTAACTCAACTGTCTCACAGATTCCTTT
T

Cez menu Tools - Table Browser ziskame GTF subor popisujuci danu oblast (interval sa automaticky prenesie z okna GB, premenime nastavenie stop na Repeats/RepeatMasker):

chr1    hg38_rmsk    exon    18151087    18151196    237.000000    +    .    gene_id "L2a"; transcript_id "L2a";
chr1    hg38_rmsk    exon    18151204    18151572    962.000000    -    .    gene_id "MLT1J2"; transcript_id "MLT1J2";
chr1    hg38_rmsk    exon    18151629    18151801    260.000000    -    .    gene_id "MLT1J2"; transcript_id "MLT1J2_dup1";
chr1    hg38_rmsk    exon    18151846    18151971    423.000000    -    .    gene_id "MER94"; transcript_id "MER94";
chr1    hg38_rmsk    exon    18151985    18152148    472.000000    +    .    gene_id "L2a"; transcript_id "L2a_dup1";

3) Upravime oba subory. Da sa samozrejme robit manualne, ale efektivnejsie hlavne u GTF je nejaka automatizovana verzia. Mne sa na taketo ad-hoc upravy osvedcilo skriptovanie cez perl na prikazovom riadku a kombinacia s nastrojmi prikazoveho riadku Linux. Upravou potrebujeme dostat z GTF -> GFF3 a upravit suradnice v nom, aby odpovedali stiahnutej sekvencii a nie celemu chr1.

Najprv skratime hlavicku fasta suboru, aby koncila na medzere:

perl -p -i -e "s/(>[^ ]+) .+/\$1/" ukol1.fa

Jedna sa o tzv. "perl one-liner" a prepinac -e sposobi, ze prikaz v uvodzovkach (substitute) sa vykona na kazdom riadku suboru (prepinace -p -i)

V GTF subore zmenime odkaz na chr1 na nazov nasej sekvencie "hg38_dna":

perl -p -i -e "s/^chr1/hg38_dna/" ukol1.gtf

Dalej podobnymi nastrojmi odcitame od vsetkych cisel zacinajucich cislicami 1815 v GTF subore hodnotu 18151079, cim ziskame suradnice odpovedajuce fasta suboru:

cat ukol1.gtf | perl -p -e "m/(1815[0-9]+).+(1815[0-9]+)/; \$beg=\$1-18151079; \$end=\$2-18151079; s/\$1\t\$2/\$beg\t\$end/;" > ukol1_recalc.gtf

Prepinac -i nepouzivame, pretoze nechceme menit subor "inline", ale chceme vytvorit pomocou vystupu novy subor. Nakoniec premenime GTF na GFF3, co znamena vytvorit hlavicku a pouzit v deviatom stlpci u premennych miesto medzery symbol "=":

echo "##gff-version 3" > ukol1.gff3
echo "##sequence-region hg38_dna 1 1100" >> ukol1.gff3
cat ukol1_recalc.gtf | perl -p -e "s/id \\\"/id=\\\"/g;" >> ukol1.gff3

4) Skontrolujeme platnost vytvoreneho GFF3 formatu online na:

http://genometools.org/cgi-bin/gff3validator.cgi

5) Vytvorime si verzie GFF3 suboru so ziadanymi prvkami. Napriklad na ziskanie vsetkych kuskov fragmentovanej sekvencie L2a, upravime takto:

grep "L2a\|\#\#" ukol1.gff3 > ukol1_l2a.gff3

6) Sekvenciu kompletneho L2a ziskate pouzitim BEDTOOLS:

bedtools getfasta -fi ukol1.fa -bed ukol1_l2a.gff3 -fo ukol1_l2a.fa

7) Sekvenciu L2a z cvicenia porovname s pripadnou sekvenciou L2a z databazy Dfam na adrese http://www.dfam.org/search