Týden 2 - Dalsi zdroje genomovych a anotacnich dat, uvod do programovani v R - 24. 2. 2020
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