ExBlast-Answers
From teachingmaterials
Contents |
Svar til BLAST øvelsen
Svar af Rasmus Wernersson, marts 2010. Opdateret af Henrik Nielsen, marts 2012, samt af Andreas Porse, marts 2013.
QUESTION 1
Q1: Report the three sequences in FASTA format (give them short UNIQUE names, e.g. "seq1", "seq2", "seq3"). :
>seq_01 CGCCCGACCGTGTAGGAGGTCCGGT >seq_02 TTCCAAGTAGGTTATGTTAAAGCCG >seq_03 CGGTTAATTCGCTGTCTATATAAGG
(of course your particular sequences will not be identical to these)
QUESTION 2
Q2a: How big (in basepairs) is the database we used for the BLAST search?:
44,780,205,193 letters, dvs 44,7 milliarder basepar.
Q2b:
- Do you find any sequences that look like your input sequences (paste in a few example alignments in your report).
Svar: Der er flere 100% homologi hits, fx:
GENE ID: 5438907 BC1G_02961 | hypothetical protein
[Botryotinia fuckeliana B05.10]
Score = 33.7 bits (36), Expect = 6.3
Identities = 18/18 (100%), Gaps = 0/18 (0%)
Strand=Plus/Minus
Query 1 CGCCCGACCGTGTAGGAG 18
||||||||||||||||||
Sbjct 574 CGCCCGACCGTGTAGGAG 557
- What is the typical length of the hits (the alignment length)?:
Svar: Typisk en 17-19 baser.
- What is the typical % identity?:
Svar: 90% - 100%
- In what range is the bit-scores ("max score)?:
Svar: typisk på 31 - 35 bits.
- What is the range of the E-values?:
Svar: 1.8 - 22 (for denne søgning, ofte kan enkelt-hits ligge helt "oppe" på 0.1 - 1). BEMÆRK: vi valgte at bruge en E-value tærskel på 50.0. Som default ligger den på 10.0 - i mit tilfælde ville der alligevel være en del hits som kunne ses. Se fx summary for seq_01:
Sequences producing significant alignments: (Bits) Value ref|XM_001558247.1| Botryotinia fuckeliana B05.10 hypothetica... 33.7 6.3 gb|AF440781.1| Streptomyces cinnamonensis polyether antibioti... 33.7 6.3 emb|AL365194.15| Human DNA sequence from clone RP4-549F15 on ... 33.7 6.3 gb|CP001349.1| Methylobacterium nodulans ORS 2060, complete g... 31.9 22 gb|CP000854.1| Mycobacterium marinum M, complete genome 31.9 22 gb|DQ665955.1| Australoplana sp. MR-2006 28S ribosomal RNA ge... 31.9 22 gb|CP000325.1| Mycobacterium ulcerans Agy99, complete genome 31.9 22 gb|CP000473.1| Solibacter usitatus Ellin6076, complete genome 31.9 22 gb|CP000248.1| Novosphingobium aromaticivorans DSM 12444, com... 31.9 22
QUESTION 3
- What is the biological significance of the hits you found / is there any biological meaning?:
Svar: Dette giver præcis INGEN biologisk mening(!) - hits'ene er rigtige nok, det er faktisk sekvens der findes i databasen. Men vi ved jo at vores sekvenser er helt tilfældige, og derfor ikke har noget slægtskab med dem i databasen. Vi har kun fundet vores hits, fordi databasen er så enormt stor, at vi af rent stokastiske årsager rammer nogle sekvenser der ligner.
E-værdierne fortæller netop dette: Som jeg gennemgik i forelæsningen, vil alignment score ("bit score" i BLAST) følge en ekstremværdi-fordeling for de sekvenser der ikke har noget med vores input-sekvenser at gøre (fordelt omkring en lav bit-score), og de hits der faktisk er relateret fordeler sig omkring en væsentlig højere bit-score.
NB: I skal ikke blive forvirrede over forskellen mellem alignment-score og bit-score; bit-score er blot alignment-score normaliseret med en konstant faktor som giver et resultat der kan udtrykkes i bits.
Ud fra ekstremværdi-fordelingen (som beregnes når BLAST databasen oprettes), er det muligt at estimere hvor mange tilfældige, ikke-relaterede hits, man må forvente ved en given bit-score. Dette udtrykkes som en e-værdi (EXPECT-value) som fortæller, hvor mange ikke-relaterede hits man kan forvente at finde ved en given bit-score (eller højere). Hvis vi fx kigger på vores ovenstående hit med en bit-score på 33.7, giver dette en e-value på 6.3 (i denne database, af 40,3*10^9 baser). Man må altså kunne forvente statistisk set at finde 6.3 hits med en score på 33.7 bits eller højere i denne database - af tilfældige årsager!
QUESTION 4
- 8,989,978,420 basepar (8,9 milliarder).
- Bit-scores: 30 - 35
- 64%-100% homologi
- Længder: 17-25; flere med gaps.
- E-values: 0.35 -15
Eksempel:
Score E Sequences producing significant alignments: (Bits) Value ref|NT_005403.17| Homo sapiens chromosome 2 genomic contig, G... 35.6 0.35 ref|NW_001838863.1| Homo sapiens chromosome 2 genomic contig,... 35.6 0.35 ref|NT_022135.16| Homo sapiens chromosome 2 genomic contig, G... 33.7 1.2 ref|NW_001838859.2| Homo sapiens chromosome 2 genomic contig,... 33.7 1.2 ref|NT_026437.12| Homo sapiens chromosome 14 genomic contig, ... 31.9 4.3 ref|NT_008413.18| Homo sapiens chromosome 9 genomic contig, G... 31.9 4.3 ref|NT_008046.16| Homo sapiens chromosome 8 genomic contig, G... 31.9 4.3 ref|NW_001838113.2| Homo sapiens chromosome 14 genomic contig... 31.9 4.3 ref|NW_001839149.2| Homo sapiens chromosome 9 genomic contig,... 31.9 4.3 ref|NW_001839136.1| Homo sapiens chromosome 8 genomic contig,... 31.9 4.3 ref|NM_019024.1| Homo sapiens HEAT repeat containing 5B (HEAT... 30.1 15 ref|NT_025028.14| Homo sapiens chromosome 18 genomic contig, ... 30.1 15 ref|NT_010859.14| Homo sapiens chromosome 18 genomic contig, ... 30.1 15 ref|NT_009952.14| Homo sapiens chromosome 13 genomic contig, ... 30.1 15 ref|NT_008818.16| Homo sapiens chromosome 10 genomic contig, ... 30.1 15 ref|NT_007933.15| Homo sapiens chromosome 7 genomic contig, G... 30.1 15 ref|NT_006576.16| Homo sapiens chromosome 5 genomic contig, G... 30.1 15 ref|NT_006316.16| Homo sapiens chromosome 4 genomic contig, G... 30.1 15 ref|NT_022517.18| Homo sapiens chromosome 3 genomic contig, G... 30.1 15 ref|NT_022184.15| Homo sapiens chromosome 2 genomic contig, G... 30.1 15 ref|NT_167186.1| Homo sapiens chromosome 1 genomic contig, GR... 30.1 15 ref|NW_001838469.1| Homo sapiens chromosome 18 genomic contig... 30.1 15 ref|NW_001838464.1| Homo sapiens chromosome 18 genomic contig... 30.1 15 ref|NW_001838084.2| Homo sapiens chromosome 13 genomic contig... 30.1 15 ref|NW_001838011.1| Homo sapiens chromosome 10 genomic contig... 30.1 15 ref|NW_001839071.2| Homo sapiens chromosome 7 genomic contig,... 30.1 15 ref|NW_001838929.1| Homo sapiens chromosome 5 genomic contig,... 30.1 15 ref|NW_001838877.2| Homo sapiens chromosome 3 genomic contig,... 30.1 15 ref|NW_001838769.1| Homo sapiens chromosome 2 genomic contig,... 30.1 15 ref|NW_001838537.2| Homo sapiens chromosome 1 genomic contig,... 30.1 15
Eksempel:
Features in this part of subject sequence:
calmodulin-binding transcription activator 1
Score = 33.7 bits (36), Expect = 1.2
Identities = 23/25 (92%), Gaps = 2/25 (8%)
Strand=Plus/Minus
Query 1 CGCCCGACCGTGTAGGAGGTCCGGT 25
|||||||||||||||||| |||||
Sbjct 3502635 CGCCCGACCGTGTAGGAG--CCGGT 3502613
QUESTION 5
Using the pre-calculated examples for this question.
Has the alignment score ("max-score") changed? Would you expect it to? :
- No: The alignment score is identical for both alignments: 33.7 bits (36).
- Since the alignment score is calculated solely from the pairwise alignment, given the scoring matrix and gap penalties, it will never chnage between BLAST databases (for identical alignments).
Has the E-value changed? Why/Why not?:
- Yes: The E-value is lower for the human-only database (Expect = 1.2) compared to the large NR database (Expect = 6.3).
- Conceptually this is easy to understand - getting an alignment with the given score (33.7 bits) is more SIGNIFICANT in the smaller database. In larger database there is a larger chance of randomly picking up matches.
What is the relationship between database size and E-value for hits with identical alignment score?
- The change in E-value i directly related to the change in database size:
- Relative difference in size of databases: 30,408,770,429 / 5,860,289,005 = 5.19
- Relative difference in E-value 6.3 / 1.2 = 5.25 (notice: E-values are rounded to 1 decimal in the BLAST output).
- Both database size and E-values differ by a factor of 5.2 in this example.
- In conclusion: Each time the database size doubles, the E-value doubles as well.
QUESTION 6
Q6: Report the sequences in FASTA format (once again use short UNIQUE names):
>seq_01 aymgsplsflshdhncirdakyvvs >seq_02 danshrsqwnrlgmanwfqhwsewv >seq_03 vsltdgpdkhesqhirqdkyfddia
QUESTION 7
Q7a:
How big is the database this time? :
Number of letters 8,171,867,726 Number of sequences 23,781,854
Svar: 8.17 milliarder bogstaver (aminosyrer).
What is the typical length of the alignment and do they contain gaps?:
Svar: Typisk 17-21. Sjældent gaps, men med en del mismatches.
What is the range of E-values?:
Svar: 200 - >900 - se fx oversigtstabellen for Query Sequence 1:
Score E Sequences producing significant alignments: (Bits) Value ref|YP_902917.1| malto-oligosyltrehalose synthase [Pelobacter... 27.7 381 gb|ACF09744.1| tRNA methylase Trm12p [uncultured marine group... 27.7 439 ref|XP_002632423.1| C. briggsae CBR-SNF-9 protein [Caenorhabd... 26.9 730
Try to inspect a few of the alignments in details ("+" means similar) - do you find any that look plausible, if we for a moment ignore the length/e-value - (you can qoute a few in your report)?
Svar: Ja se fx. følgende alignment - det ser da overbevisende ud (men det er alt, alt for kort til at være signifikant, og E-værdien er alt, alt for høj):
Score = 27.7 bits (60), Expect = 381, Method: Composition-based stats.
Identities = 10/17 (58%), Positives = 12/17 (70%), Gaps = 0/17 (0%)
Query 1 AYMGSPLSFLSHDHNCI 17
AY GSP + HDHNC+
Sbjct 49 AYPGSPHGYDIHDHNCL 65
If we had used the default E-value cutoff of 10 would any hits have been found?:
Svar: NEJ. Hvis vi også havde kørt søgningerne mod DNA databasen med en e-value cut-off på 1000, ville vi har set sidevis af hits for hver sekvens
Q7b:
If we compare the result from BLAST'ing random DNA sequences to random Peptide sequences - what kind if search has the higher risk of returning false positives (results that appear plausible, maybe even significant, but are truly unrelated)?:
Svar: Risikoen for at samle et "falsk" hit op, er større med DNA - bemærk at vi var nødt til at "snyde" temmeligt meget med cut-off for i det hele taget at se nogle hits mod proteinsekvenserne.
Kommentar ang. word-size og proteinsekvenser: (Vi gennemgik også dette i forelæsningen): BLASTP fungerer noget anderledes end BLASTN mht. at finde et sted i sekvenserne med høj-similaritet: i stedet for at kræve et længere stykke der er 100% det samme (det vil ikke give mening for protein-sekvens), prøver algoritmen at finde to sekvens-stumper (af længde "word-size", default =3) med høj similaritet inden for en kort afstand (default: 40 aa). Så snart et sådant sekvenspar er fundet, fylder BLAST ind med almindeligt lokalt alignment.
QUESTION 8
Do we get any significant hits?
Svar: ja, der er 7 hits med en e-value på "0.0" (dvs så lille at den bliver afrundet til nul) - og de næste er også yderst signifikante. Det første hit (S48754) har desuden en coverage på 100% og en identity på 100%.
Score E Sequences producing significant alignments: (Bits) Value gb|S48754.1| 221 protease=alkaline serine protease [Bacillus ... 2547 0.0 dbj|D13157.1|BACAK221 Bacillus sp. alkaline protease gene 2538 0.0 gb|M65086.1|BACALKPR B.alcalophilus alkaline protease gene, c... 2241 0.0 gb|FJ940727.1| Bacillus alcalophilus strain TCCC11004 alkalin... 2091 0.0 dbj|AP006627.1| Bacillus clausii KSM-K16 DNA, complete genome 2080 0.0 gb|M28537.1|BACYAB Bacillus subtilis (YaB) alkaline elastase ... 677 0.0 dbj|AB005792.1| Bacillus sp. DNA for AprN, complete cds 646 0.0 dbj|D29688.1|BACAPRS Bacillus sp. G-825-6 aprS gene for subti... 477 8e-131 gb|EU391617.1| Expression vector pCL311, complete sequence 403 1e-108 gb|EU130686.1| Bacillus alcalophilus strain PB92 alkaline pro... 255 5e-64 dbj|D29736.1|BACAPRQ Bacillus sp. aprQ gene for subtilisin AL... 123 2e-24 gb|CP001878.1| Bacillus pseudofirmus OF4, complete genome 120 3e-23 dbj|BA000004.3| Bacillus halodurans C-125 DNA, complete genome 118 9e-23 dbj|D26542.1|BACAPRM Bacillus sp. gene for serine protease Ap... 118 9e-23 gb|DQ868543.1| Bacillus clausii strain XJU-5 thermostable alk... 114 1e-21 dbj|D13158.2|BACTAP Bacillus halodurans gene for prepro-therm... 113 4e-21 gb|AY953434.1| Bacillus subtilis clone pGXAA2011 subtilisin A... 111 1e-20
What kind of genes (function) do we find?:
Svar: De er alle basiske ("alkaline") serin-protease fra Bacillus -> vores DNA stammer faktisk fra det første entry ("S48754").
How is the distribution of E-values - a broad spectrum from good to bad, or a sharp division into good/bad scores?:
Svar: Der er en række (næsten) perfekte hits (bit-score: 2100-2700), et par stykker der er signifikante men korte (bit-score: 75-100, e-value: 1e-10 til 1e-18) og en hel masse junk, der består af korte fragmenter med en dårlig e-value. Altså, en rimelig brat overgang fra de gode hits til de virkeligt dårlige.
Som gennemgået i forelæsningen, er DNA alignment-matricen optimeret til at finde hits der matcher godt - bemærk at mismatches straffes hårdt (Bemærk at NCBI også har en 2/-3 matrice - men effekten er basalt set den samme):
a c g t a 1 c -3 1 g -3 -3 1 t -3 -3 -3 1
QUESTION 9
BEMÆRK: HUSK at bruger ORF finderen! Eftersom vi ved at det er et fuld-længde transcript, kan vi godt regne med at både START og STOP codon er med. Vi kan derfor vælge "START-codon = strict".
Note: man kan oversætte med både den standard-genetiske kode (tabel 1) eller alternativt med tabel 11 ("The Bacterial and Plant Plastid Code"). Den eneste forskel er at tabel 11 tillader nogle ekstra yderst sjældne START codons, bla. fra bacterielle phager.
Report your translated protein sequence in FASTA format.:
>Unknown_transcript01_rframe2_ORF MKKPLGKIVASTALLISVAFSSSIASAAEEAKEKYLIGFNEQEAVSEFVEQVEANDEVAI LSEEEEVEIELLHEFETIPVLSVELSPEDVDALELDPAISYIEEDAEVTTMAQSVPWGIS RVQAPAAHNRGLTGSGVKVAVLDTGISTHPDLNIRGGASFVPGEPSTQDGNGHGTHVAGT IAALNNSIGVLGVAPSAELYAVKVLGASGSGSVSSIAQGLEWAGNNGMHVANLSLGSPSP SATLEQAVNSATSRGVLVVAASGNSGAGSISYPARYANAMAVGATDQNNNRASFSQYGAG LDIVAPGVNVQSTYPGSTYASLNGTSMATPHVAGAAALVKQKNPSWSNVQIRNHLKNTAT SLGSTNLYGSGLVNAEAATR
Do we find any conserved protein-domains? (Indicated at the very top of the result page, and during the search).:
Svar: Ja, der findes et "Peptidase S8" domæne. Dette kan se på både "pause-siden" mens der søges, og ved at klikke på "Show conserved domains" øverst på resultat siden.
Do we find any significant hits? (E-value?):
Svar: Ja, ganske mange. (Hvis vi fx. kigger på de første 15 ligger E-value fra 0.0 – 2e-147 – bemærk at det per default kun er de første 100 hits der vises).
Score E Sequences producing significant alignments: (Bits) Value sp|P41362.1|ELYA_BACCS RecName: Full=Alkaline protease; Flags... 754 0.0 sp|P27693.1|ELYA_BACAO RecName: Full=Alkaline protease; Flags... 753 0.0 ref|YP_174261.1| extracellular alkaline serine protease [Baci... 745 0.0 sp|P20724.1|ELYA_BACYA RecName: Full=Alkaline elastase YaB; F... 597 1e-168 dbj|BAA25184.1| AprN [Bacillus sp.] 595 4e-168 dbj|BAA06157.1| prepro-subtilisin Sendai [Bacillus sp. G-825-6] 561 4e-158 sp|P29600.1|SUBS_BACLE RecName: Full=Subtilisin Savinase; Alt... 533 1e-149 pdb|1AH2|A Chain A, Serine Protease Pb92 From Bacillus Alcalo... 531 3e-149 pdb|1GCI|A Chain A, The 0.78 Angstroms Structure Of A Serine ... 530 1e-148 pdb|1IAV|A Chain A, Structure On Native (Asn 87) Subtilisin F... 530 1e-148 pdb|1C9M|A Chain A, Bacillus Lentus Subtilsin (Ser 87) N76dS1... 528 5e-148 pdb|1WSD|A Chain A, Alkaline M-Protease Form I Crystal Strcut... 527 7e-148 pdb|1NDU|A Chain A, Bacillus Lentus Subtilisin Variant S101gV... 526 1e-147 pdb|1C9J|A Chain A, Bacillus Lentus Subtilisin K27rN87SV104YN... 525 3e-147 sp|P29599.1|SUBB_BACLE RecName: Full=Subtilisin BL; AltName: ... 525 4e-147
Are all the best hits the same category of enzymes?:
Svar: Proteaser og Peptidaser (hvilket dækker over det samme). Bemærk at man kan klikke på hits'ne og gå direkte til sekvensen, hvor man kan læse yderligere.
Vi kan slå følgende fast:
- Det er en serin protease.
- Den tilhører S8 familien
- Den er basisk
- Den kommer fra Bacillus.
How does the distribution of E-values looks compared to the DNA search?:
Svar: Fra fuld-længde hits med en e-value på mindre end 1e-168 og gradvis mod højere værdier (per default ser vi kun de 100 bedste hits).
From what you have seen, what is best for identifying intermediate quality hits - DNA or Protein BLAST?:
Svar: BLASTP
QUESTION 10
Svar: Vi vil i sidste ende gerne bruge BLAST til at søge i de store samlede databaser.
Lad os derfor prøve følgende:
- BLASTN
- Oversætte til protein (fx. VirtualRibosome).
- BLASTP
Både med BLASTN og BLASTP ønsker vi at bruge NR databasen - det er den der giver den største bredde. Der er godt at starte med en stor, bred database - så kan man evt. gå i dybten i en af de mere specifikke, når man har skudt sig ind på hvad for en type sekvens man har.
Til at starte med har jeg konverteret DNA sekvensen til FASTA format (brug de jEdit tricks vi lærte første gang):
>CLONE12 AACGGGCACGGGACGCATGTAGCTGGAACAGTGGCAGCCGTAAATAATAATGGTATCGGA GTTGCCGGGGTTGCAGGAGGAAACGGCTCTACCAATAGTGGAGCAAGGTTAATGTCCACA CAAATTTTTAATAGTGATGGGGATTATACAAATAGCGAAACTCTTGTGTACAGAGCCATT GTTTATGGTGCAGATAACGGAGCTGTGATCTCGCAAAATAGCTGGGGTAGTCAGTCTCTG ACTATTAAGGAGTTGCAGAAAGCTGCGATCGACTATTTCATTGATTATGCAGGAATGGAC GAAACAGGAGAAATACAGACAGGCCCTATGAGGGGAGGTATATTTATAGCTGCCGCCGGA AACGATAACGTTTCCACTCCAAATATGCCTTCAGCTTATGAACGGGTTTTAGCTGTGGCC TCAATGGGACCAGATTTTACTAAGGCAAGCTATAGCACTTTTGGAACATGGACTGATATT ACTGCTCCTGGCGGAGATATTGACAAATTTGATTTGTCAGAATACGGAGTTCTCAGCACT TATGCCGATAATTATTATGCTTATGGAGAGGGAACATCCATGGCTTGTCCACATGTCGCC GGCGCCGCC
1) BLASTN - når vi prøver at BLAST'e mod NR får vi slet ikke nogen signifikante hits:
Database: All GenBank+EMBL+DDBJ+PDB sequences (but no EST, STS,
GSS,environmental samples or phase 0, 1 or 2 HTGS sequences)
11,076,294 sequences; 30,408,770,429 total letters
Query= CLONE12
Length=609
Score E
Sequences producing significant alignments: (Bits) Value
gb|CP000159.1| Salinibacter ruber DSM 13855, complete genome 51.8 0.004
ref|XM_002461822.1| Sorghum bicolor hypothetical protein, mRNA 50.0 0.016
gb|CP000813.1| Bacillus pumilus SAFR-032, complete genome 50.0 0.016
ref|XM_002526491.1| Ricinus communis Cucumisin precursor, put... 48.2 0.054
gb|CP001177.1| Bacillus cereus AH187, complete genome 48.2 0.054
gb|CP000485.1| Bacillus thuringiensis str. Al Hakam, complete... 48.2 0.054
gb|CP000001.1| Bacillus cereus E33L, complete genome 48.2 0.054
dbj|AB085752.1| Bacillus sp. KSM-LD1 gene for protease, compl... 48.2 0.054
dbj|AK320060.1| Solanum lycopersicum cDNA, clone: LEFL1004CD0... 44.6 0.66
gb|BT060461.1| Zea mays full-length cDNA clone ZM_BFb0002M22 ... 44.6 0.66
gb|CP000227.1| Bacillus cereus Q1, complete genome 44.6 0.66
ref|NM_001152466.1| Zea mays hypothetical protein LOC10027946... 44.6 0.66
gb|EU957086.1| Zea mays clone 1582547 mRNA sequence 44.6 0.66
ref|XM_001193501.1| PREDICTED: Strongylocentrotus purpuratus ... 44.6 0.66
ref|XM_001196149.1| PREDICTED: Strongylocentrotus purpuratus ... 44.6 0.66
gb|BT014524.1| Lycopersicon esculentum clone 133938F, mRNA se... 44.6 0.66
gb|CP001899.1| Ferroglobus placidus DSM 10642, complete genome 42.8 2.3
gb|BT117303.1| Picea glauca clone GQ03816_N01 mRNA sequence 42.8 2.3
ref|XM_002274806.1| PREDICTED: Vitis vinifera hypothetical pr... 42.8 2.3
ref|XM_002316218.1| Populus trichocarpa predicted protein, mRNA 42.8 2.3
ref|XM_002316216.1| Populus trichocarpa predicted protein, mRNA 42.8 2.3
(TIP: Man kan få dette tekst-baserede output ved at klikke på "Format/reformat" i toppen af siden, og så vælge "Plain text").
Der er simpelt hen ikke noget i databasen, der er ligner vores input tilstrækkeligt - derfor finder vi kun de tilfældige resultater. En søgning på DNA niveau er kun egnet til at finde de sekvenser, der ligner virkeligt meget.
2) Jeg oversætter med Virtual ribosome.
BEMÆRK: vi ved følgende om sekvensen:
- Den er proteinkodende (amplificeret med PCR primere rettet mod en særlig gruppe af enzymer).
- Det er et fragment – kun midten af genet var amplificeret – vi kan derfor ikke forvente at STOP og START codon er med.
Resultat fra ORF finderen:
VIRTUAL RIBOSOME
----------------
Translation table: Standard SGC0
>CLONE12_rframe1_ORF
Reading frame: 1
N G H G T H V A G T V A A V N N N G I G V A G V A G G N G S
5' AACGGGCACGGGACGCATGTAGCTGGAACAGTGGCAGCCGTAAATAATAATGGTATCGGAGTTGCCGGGGTTGCAGGAGGAAACGGCTCT 90
..........................................................................................
T N S G A R L M S T Q I F N S D G D Y T N S E T L V Y R A I
5' ACCAATAGTGGAGCAAGGTTAATGTCCACACAAATTTTTAATAGTGATGGGGATTATACAAATAGCGAAACTCTTGTGTACAGAGCCATT 180
.....................>>>..................................................................
V Y G A D N G A V I S Q N S W G S Q S L T I K E L Q K A A I
5' GTTTATGGTGCAGATAACGGAGCTGTGATCTCGCAAAATAGCTGGGGTAGTCAGTCTCTGACTATTAAGGAGTTGCAGAAAGCTGCGATC 270
.........................................................)))............)))...............
D Y F I D Y A G M D E T G E I Q T G P M R G G I F I A A A G
5' GACTATTTCATTGATTATGCAGGAATGGACGAAACAGGAGAAATACAGACAGGCCCTATGAGGGGAGGTATATTTATAGCTGCCGCCGGA 360
........................>>>..............................>>>..............................
N D N V S T P N M P S A Y E R V L A V A S M G P D F T K A S
5' AACGATAACGTTTCCACTCCAAATATGCCTTCAGCTTATGAACGGGTTTTAGCTGTGGCCTCAATGGGACCAGATTTTACTAAGGCAAGC 450
........................>>>....................................>>>........................
Y S T F G T W T D I T A P G G D I D K F D L S E Y G V L S T
5' TATAGCACTTTTGGAACATGGACTGATATTACTGCTCCTGGCGGAGATATTGACAAATTTGATTTGTCAGAATACGGAGTTCTCAGCACT 540
...............................................................)))........................
Y A D N Y Y A Y G E G T S M A C P H V A G A A
5' TATGCCGATAATTATTATGCTTATGGAGAGGGAACATCCATGGCTTGTCCACATGTCGCCGGCGCCGCC 609
.......................................>>>...........................
(TIP: Husk at du kan få sekvensen i FASTA format ved at klikke på "FASTA" linket på resultat-siden).:
>CLONE12_rframe1_ORF NGHGTHVAGTVAAVNNNGIGVAGVAGGNGSTNSGARLMSTQIFNSDGDYTNSETLVYRAI VYGADNGAVISQNSWGSQSLTIKELQKAAIDYFIDYAGMDETGEIQTGPMRGGIFIAAAG NDNVSTPNMPSAYERVLAVASMGPDFTKASYSTFGTWTDITAPGGDIDKFDLSEYGVLST YADNYYAYGEGTSMACPHVAGAA
3) BLASTP
Det første vi ser, er at at der findes et "Peptidase_S8" domæne. Dette er en meget kraftig indikator på funktionen.
Vi får en række signifikante hits. Hvis vi kigger på de bedste, ser at de alle er serin-proteaser (af S8 familien):
Database: All non-redundant GenBank CDS
translations+PDB+SwissProt+PIR+PRF excluding environmental samples
from WGS projects
10,570,301 sequences; 3,602,205,473 total letters
Query= CLONE12_rframe1_ORF
Length=203
Score E
Sequences producing significant alignments: (Bits) Value
ref|ZP_05857871.1| subtilase family domain protein [Prevotell... 194 6e-48
ref|ZP_06407789.1| subtilase family domain protein [Prevotell... 187 6e-46
ref|ZP_04833939.1| subtilase family domain protein [Prevotell... 187 6e-46
ref|ZP_03300901.1| hypothetical protein BACDOR_02271 [Bactero... 179 1e-43
ref|ZP_05253523.1| protease [Bacteroides sp. 4_3_47FAA] >gb|E... 177 7e-43
ref|ZP_06089720.1| protease [Bacteroides sp. 3_1_33FAA] >gb|E... 177 7e-43
ref|ZP_04555416.1| protease [Bacteroides sp. D4] >gb|EEO46750... 177 7e-43
ref|ZP_04539947.1| protease [Bacteroides sp. 9_1_42FAA] >gb|E... 174 4e-42
ref|ZP_05734780.1| subtilase family domain protein [Prevotell... 171 5e-41
ref|ZP_06251999.1| subtilase family protein [Prevotella copri... 165 2e-39
ref|ZP_03460514.1| hypothetical protein BACEGG_03331 [Bactero... 165 3e-39
ref|ZP_05917065.1| subtilase family domain protein [Prevotell... 164 6e-39
ref|ZP_06422877.1| subtilase family domain protein [Prevotell... 159 1e-37
ref|ZP_06006082.2| conserved hypothetical protein [Prevotella... 156 1e-36
ref|ZP_05416921.1| subtilase family domain protein [Bacteroid... 156 1e-36
ref|NP_809125.1| protease [Bacteroides thetaiotaomicron VPI-5... 154 6e-36
(BEMÆRK: "Hypotetical" proteiner er ikke nær så stærk evidens som et velbeskrevet protein).
Lad os kigge lidt nærmere på det bedste hit:
>ref|ZP_05857871.1| subtilase family domain protein [Prevotella veroralis F0319]
gb|EEX18266.1| subtilase family domain protein [Prevotella veroralis F0319]
Length=921
Score = 194 bits (492), Expect = 6e-48, Method: Compositional matrix adjust.
Identities = 117/214 (54%), Positives = 138/214 (64%), Gaps = 20/214 (9%)
Query 2 GHGTHVAGTVAAVNNNGIGVAGVAGGNGSTNSGARLMSTQIF---NSDGDYTNSETLVYR 58
GHGTHVAGTVAA NNNG GVAG+AGGNG+ +SG RLMS QIF N GD +
Sbjct 283 GHGTHVAGTVAARNNNGKGVAGIAGGNGTADSGVRLMSCQIFRKKNEQGDAAAAIKY--- 339
Query 59 AIVYGADNGAVISQNSWGSQSLT----IKELQKAAIDYFIDYAGMDETGEIQT-GPMRGG 113
ADNGAVI QNSWG S T + +L K A+DYFI AG D G+ + PM+GG
Sbjct 340 ----AADNGAVICQNSWGYSSNTGVTSMPQLLKDAVDYFIKNAGCDANGQQRADSPMKGG 395
Query 114 IFIAAAGNDNVSTPNMPSAYERVLAVASMGPDFTKASYSTFGTWTDITAPGGDIDKFDLS 173
+ I AAGN+N P+ Y ++VASM DFTKASYS + W ITAPGGD D+F +
Sbjct 396 VVIFAAGNENKEFSAYPACYPPAVSVASMAWDFTKASYSNYAKWVTITAPGGDQDRFG-N 454
Query 174 EYGVLSTY----ADNYYAYGEGTSMACPHVAGAA 203
E GVLST A + YAY +GTSMACPHV+G A
Sbjct 455 EAGVLSTVPKSKASSGYAYMQGTSMACPHVSGIA 488
Bemærk at selv om det ikke er et perfekt hit (vores input-sekvens findes jo netop ikke i databasen), så ser det ganske fornuftigt ud: alignment'et er langt med en Identity på 54% og en Similarity på 64%.
Sammenholdt med det faktum, at alle de bedste hits er serin-proteaser, har vi en meget stærk indikation på hvilken type af enzym vi har med at gøre.
Konklusion: Sekvensen koder for en serin-protease, fra S8 familien.
QUESTION 11
(HUSK at det er i SGD og ikke GenBank I skal slå det op.)
Svar: HTA2 og HTA1 er næsten ens ("nearly identical").
Protein sekvens (Vælg “ORF translation -> get sequence):
>YBL003C Chr 2 reverse complement MSGGKGGKAGSAAKASQSRSAKAGLTFPVGRVHRLLRRGNYAQRIGSGAPVYLTAVLEYL AAEILELAGNAARDNKKTRIIPRHLQLAIRNDDELNKLLGNVTIAQGGVLPNIHQNLLPK KSAKTAKASQEL*
QUESTION 12
How many high-confidence hits do we get?:
Svar: 3
Database: Saccharomyces cerevisiae RefSeq protein
5885 sequences; 2,916,286 total letters
Query= YBL003C Chr 2 reverse complement
Length=133
Score E
Sequences producing significant alignments: (Bits) Value
ref|NP_009552.1| Histone H2A, core histone protein required f... 186 8e-49
ref|NP_010511.1| Histone H2A, core histone protein required f... 183 9e-48
ref|NP_014631.1| Histone variant H2AZ, exchanged for histone ... 89.7 1e-19
Does the hits make sense, from what you have read about HTA2 at the SGD webpage?:
Svar: Ja - HTA1 og HTA2 er meget ens (kun 2 aminosyrer er forskellige). Som bonus finder vi også "H2AZ" - en anden variant af HTA2 der noget mere divergeret i sekvens.
QUESTION 13
How many high-confidence hits are found?:
Database: Homo sapiens RefSeq protein
38,716 sequences; 18,793,054 total letters
Query= YBL003C Chr 2 reverse complement
Length=133
Score E
Sequences producing significant alignments: (Bits) Value
ref|NP_003503.1| histone H2A type 1-C [Homo sapiens] 194 1e-50
ref|NP_003507.1| histone H2A type 2-A [Homo sapiens] >ref|NP_... 194 2e-50
ref|NP_003508.1| histone H2A type 2-C [Homo sapiens] 194 2e-50
ref|NP_003500.1| histone H2A type 1 [Homo sapiens] >ref|NP_00... 194 2e-50
ref|NP_734466.1| histone H2A type 1-A [Homo sapiens] 194 2e-50
ref|NP_542163.1| histone H2A type 1-H [Homo sapiens] 194 2e-50
ref|NP_254280.1| histone H2A type 3 [Homo sapiens] 194 3e-50
ref|NP_808760.1| histone H2A.J [Homo sapiens] 193 4e-50
ref|NP_066544.1| histone cluster 1, H2aj [Homo sapiens] 193 4e-50
ref|NP_066409.1| histone H2A type 1-D [Homo sapiens] 193 4e-50
ref|NP_066390.1| histone H2A type 1-B/E [Homo sapiens] >ref|N... 192 5e-50
ref|NP_002096.1| histone H2A.x [Homo sapiens] 192 5e-50
ref|NP_778235.1| histone H2A type 2-B [Homo sapiens] 192 7e-50
ref|NP_061119.1| core histone macro-H2A.2 [Homo sapiens] 153 3e-38
ref|NP_613075.1| core histone macro-H2A.1 isoform 1 [Homo sap... 152 7e-38
ref|NP_004884.1| core histone macro-H2A.1 isoform 2 [Homo sap... 152 8e-38
ref|NP_613258.2| core histone macro-H2A.1 isoform 3 [Homo sap... 152 8e-38
ref|NP_036544.1| histone H2A.V isoform 1 [Homo sapiens] 133 4e-32
ref|NP_002097.1| histone H2A.Z [Homo sapiens] 132 6e-32
ref|NP_619541.1| histone H2A.V isoform 2 [Homo sapiens] 114 2e-26
ref|NP_958844.1| histone H2A.V isoform 3 [Homo sapiens] 112 8e-26
ref|XP_001721252.1| PREDICTED: similar to hCG1643984 [Homo sa... 91.7 2e-19
ref|NP_542451.1| histone H2A-Bbd type 2/3 [Homo sapiens] >ref... 89.0 1e-18
ref|NP_001017990.1| histone H2A-Bbd type 1 [Homo sapiens] 85.9 9e-18
ref|NP_958925.1| histone H2A.V isoform 5 [Homo sapiens] 75.9 8e-15
ref|NP_958924.1| histone H2A.V isoform 4 [Homo sapiens] 62.8 7e-11
ref|NP_005624.2| son of sevenless homolog 1 [Homo sapiens] 49.3 8e-07
ref|NP_001018082.1| ankyrin repeat and BTB/POZ domain-contain... 43.9 4e-05
ref|NP_665803.1| ankyrin repeat and BTB/POZ domain-containing... 38.9 0.001
Lad os ligesom før kræve at hits skal have en e-value på 1e-10 eller bedre (mindre) for at vi regner det for pålideligt. Dette passer også fint med, at disse hits er enige om at der er tale om en histon.
Som udgangspunkt er der derfor 26 gode hits. Man kan dog godt argumentere for at man skal passe på hits der kun har en forudsagt funktion ("PREDICTED") - dem har vi en enkelt af. Som det også ses er der en række af hits'ne der dækker over variationer af samme protein (fx alle dem der hedder: histone H2A type-XYZ). Det er ikke altid at alle detaljerne kommer med i den korte form af overskrifter - nogen gange kan det være nødvendigt mauelt at inspicere et hit (klikke på link'et og læse hvad der står af information i det bagvedliggende database entry).
Lad os arbejde videre med alle 26 gode hits i de næste analyser.
These protein originates from a number of genes - but how many UNIQUE genes?:
Svar: Trick'et er her af en del af hits'ne er fra isoformer af det sammen protein. Et protein med to isofomer stammer stadig kun fra en enkelt gen
Lad os fx kigge på følgende hits:
ref|NP_613075.1| core histone macro-H2A.1 isoform 1 [Homo sap... 152 7e-38 ref|NP_004884.1| core histone macro-H2A.1 isoform 2 [Homo sap... 152 8e-38
Hvis man gå ind og læser sekvens-entry'et for disse hits kan man faktisk direkte se, at de stammer fra sammen gen (/gene="H2AFY", /gene="H2AFY"). BEMÆRK: Her er vi heldige og alt den information vi har burg for står faktisk i overskrifterne.
En hurtig optælling afslører følgende isoformet i vores set af BLAST hits:
ref|NP_613075.1| core histone macro-H2A.1 isoform 1 [Homo sap... 152 7e-38 ref|NP_004884.1| core histone macro-H2A.1 isoform 2 [Homo sap... 152 8e-38 ref|NP_613258.2| core histone macro-H2A.1 isoform 3 [Homo sap... 152 8e-38
ref|NP_036544.1| histone H2A.V isoform 1 [Homo sapiens] 133 4e-32 ref|NP_619541.1| histone H2A.V isoform 2 [Homo sapiens] 114 2e-26 ref|NP_958844.1| histone H2A.V isoform 3 [Homo sapiens] 112 8e-26 ref|NP_958925.1| histone H2A.V isoform 5 [Homo sapiens] 75.9 8e-15 ref|NP_958924.1| histone H2A.V isoform 4 [Homo sapiens] 62.8 7e-11
Dvs at vi i alt har 27 - 2 - 4 = 21 unikke hits til gener.
QUESTION 14
Langt svar: Et lavt komplekst område er et stykke sekvens der ikke indeholder særligt meget information (fx. TTTTAAAA i human – findes milioner af gange i genomet).
BLAST har et indbygget filter der masker disse områder ud i søgningen.
Ud over at kunne slå det til og fra – har NCBI valgt et tredie mulighed (som default) – slå det fra men at vise hvor områderne er (små bogstaver):
>ref|NP_003503.1| histone H2A type 1-C [Homo sapiens]
Length=130
GENE ID: 8334 HIST1H2AC | histone cluster 1, H2ac [Homo sapiens]
(Over 10 PubMed links)
Score = 194 bits (494), Expect = 1e-50, Method: Compositional matrix adjust.
Identities = 96/125 (76%), Positives = 108/125 (86%), Gaps = 0/125 (0%)
Query 4 gkggkagsaakasqsrsakagLTFPVGRVHRLLRRGNYAQRIGSGAPVYLTavleylaae 63
G+G + G A ++SRS++AGL FPVGRVHRLLR+GNYA+R+G+GAPVYL AVLEYL AE
Sbjct 3 GRGKQGGKARAKAKSRSSRAGLQFPVGRVHRLLRKGNYAERVGAGAPVYLAAVLEYLTAE 62
Query 64 ilelaGNAARDNKKTRIIPRHLQLAIRNDDELNKLLGNVTIAQGGVLPNIHQNLLPKKSA 123
ILELAGNAARDNKKTRIIPRHLQLAIRND+ELNKLLG VTIAQGGVLPNI LLPKK+
Sbjct 63 ILELAGNAARDNKKTRIIPRHLQLAIRNDEELNKLLGRVTIAQGGVLPNIQAVLLPKKTE 122
Query 124 KTAKA 128
KA
Sbjct 123 SHHKA 127
Hvis man gentager søgningen med filteret slået til giver det følgende (jeg har også valgt at viser X’er i det filtrerede område – hvilket NORMALT er standard i BLAST):
>ref|NP_003503.1| histone H2A type 1-C [Homo sapiens]
Length=130
GENE ID: 8334 HIST1H2AC | histone cluster 1, H2ac [Homo sapiens]
(Over 10 PubMed links)
Score = 147 bits (372), Expect = 2e-36, Method: Compositional matrix adjust.
Identities = 87/104 (83%), Positives = 93/104 (89%), Gaps = 0/104 (0%)
Query 25 LTFPVGRVHRLLRRGNYAQRIGSGAPVYLTXXXXXXXXXXXXXXGNAARDNKKTRIIPRH 84
L FPVGRVHRLLR+GNYA+R+G+GAPVYL AVLEYL AEILELAGNAARDNKKTRIIPRH
Sbjct 24 LQFPVGRVHRLLRKGNYAERVGAGAPVYLAAVLEYLTAEILELAGNAARDNKKTRIIPRH 83
Query 85 LQLAIRNDDELNKLLGNVTIAQGGVLPNIHQNLLPKKSAKTAKA 128
LQLAIRND+ELNKLLG VTIAQGGVLPNI LLPKK+ KA
Sbjct 84 LQLAIRNDEELNKLLGRVTIAQGGVLPNIQAVLLPKKTESHHKA 127
Bemærk det kortere alignment og den ændrede E-værdi.
Hvis man eksplicit slår filteret fra, ser resultatet således ud:
>ref|NP_003503.1| histone H2A type 1-C [Homo sapiens]
Length=130
GENE ID: 8334 HIST1H2AC | histone cluster 1, H2ac [Homo sapiens]
(Over 10 PubMed links)
Score = 194 bits (494), Expect = 1e-50, Method: Compositional matrix adjust.
Identities = 96/125 (76%), Positives = 108/125 (86%), Gaps = 0/125 (0%)
Query 4 GKGGKAGSAAKASQSRSAKAGLTFPVGRVHRLLRRGNYAQRIGSGAPVYLTAVLEYLAAE 63
G+G + G A ++SRS++AGL FPVGRVHRLLR+GNYA+R+G+GAPVYL AVLEYL AE
Sbjct 3 GRGKQGGKARAKAKSRSSRAGLQFPVGRVHRLLRKGNYAERVGAGAPVYLAAVLEYLTAE 62
Query 64 ILELAGNAARDNKKTRIIPRHLQLAIRNDDELNKLLGNVTIAQGGVLPNIHQNLLPKKSA 123
ILELAGNAARDNKKTRIIPRHLQLAIRND+ELNKLLG VTIAQGGVLPNI LLPKK+
Sbjct 63 ILELAGNAARDNKKTRIIPRHLQLAIRNDEELNKLLGRVTIAQGGVLPNIQAVLLPKKTE 122
Query 124 KTAKA 128
KA
Sbjct 123 SHHKA 127
Præcis samme længde, e-værdi mm. som i det første alignment.
Kort svar: Ja, alignmentet bliver kortere når low-complexity filteret er slået til.
