Difference between revisions of "Fabry:Sequence alignments (sequence searches and multiple alignments)/Journal"
From Bioinformatikpedia
Rackersederj (talk | contribs) m |
Rackersederj (talk | contribs) (→Comparison) |
||
Line 115: | Line 115: | ||
=== Comparison === |
=== Comparison === |
||
− | >R |
+ | >R CMD BATCH all_Evalues.R |
− | #This R Script is based on Andrea's Rscript psiBlast.evalueHist.Rscript |
||
− | # Thank you Andrea :) ... Julia |
||
− | |||
− | |||
− | library("animation") |
||
− | |||
− | |||
− | |||
− | cols3 <- hcl(h = seq(30, by=300 / 3, length = 3), l = 65, alpha = 0.5) |
||
− | br <- seq(-180,15,1.25) |
||
− | Blast <- read.table("Blast/blastsearch_default_v700.out_ident_cov.tsv", header = T) |
||
− | Psi <- read.table("Blast/psi_results_10its_eVal_0.000000001.txt_ident_cov.tsv", header = T) |
||
− | HHBlits <- read.table("HHBlits/hhblits_n8_neu.out_ident_cov.tsv", header = T) |
||
− | |||
− | |||
− | saveMovie({ |
||
− | for (i in 1:4) { |
||
− | if(i == 1){ |
||
− | data <- Blast$Evalue |
||
− | } |
||
− | if(i == 2){ |
||
− | data <- Psi$Evalue |
||
− | } |
||
− | if(i == 3){ |
||
− | data <- HHBlits$Evalue |
||
− | } |
||
− | if(i < 4) { |
||
− | hist(data ,border=cols3[i],col=cols3[i],breaks=br,main='E-values',panel.first = grid(),xlab='log_10(E- |
||
− | value)',axes=TRUE,xlim=range(-200,1),ylim=range(0,200)) |
||
− | legend(-200,200,c("BLAST","Psi-BLAST","HHBlits"),col=cols3, pch = c("_","_","_"), lwd = 3) |
||
− | } |
||
− | if(i == 4) { |
||
− | hist(Blast$Evalue,border=cols3[1],col=cols3[1],breaks=br,main='E-values',panel.first = grid(),xlab='log_10(E- |
||
− | value)',axes=TRUE,xlim=range(-200,1),ylim=range(0,200)) |
||
− | hist(Psi$Evalue,border=cols3[2],col=cols3[2],breaks=br,xlim=range(-200,1),ylim=range(0,200),add=TRUE) |
||
− | hist(HHBlits$Evalue,border=cols3[3],col=cols3[3],breaks=br,xlim=range(-200,1),ylim=range(0,200),add=TRUE) |
||
− | legend(-200,200,c("BLAST","Psi-BLAST","HHBlits"),col=cols3, pch = c("_","_","_"), lwd = 3) |
||
− | Sys.sleep(5) |
||
− | } |
||
− | } |
||
− | }, interval = 5, outdir = getwd()) |
Revision as of 07:40, 5 May 2012
Please see Task 2 Results for our results on this topic. Please see also Task 2 Scripts for the used scripts.
Contents
Reference sequence
The reference sequence of α-Galactosidase A that will be used in the following tasks was obtained from Swissprot P06280.
>gi|4504009|ref|NP_000160.1| alpha-galactosidase A precursor [Homo sapiens] MQLRNPELHLGCALALRFLALVSWDIPGARALDNGLARTPTMGWLHWERFMCNLDCQEEPDSCISEKLFM EMAELMVSEGWKDAGYEYLCIDDCWMAPQRDSEGRLQADPQRFPHGIRQLANYVHSKGLKLGIYADVGNK TCAGFPGSFGYYDIDAQTFADWGVDLLKFDGCYCDSLENLADGYKHMSLALNRTGRSIVYSCEWPLYMWP FQKPNYTEIRQYCNHWRNFADIDDSWKSIKSILDWTSFNQERIVDVAGPGGWNDPDMLVIGNFGLSWNQQ VTQMALWAIMAAPLFMSNDLRHISPQAKALLQDKDVIAINQDPLGKQGYQLRQGDNFEVWERPLSGLAWA VAMINRQEIGGPRSYTIAVASLGKGVACNPACFITQLLPVKRKLGFYEWTSRLRSHINPTGTVLLQLENT MQMSLKDLL
Sequence searches
Blast
We searched the "big80" database with Blast with the following command:
blastall -p blastp -d /mnt/project/pracstrucfunc12/data/big/big_80 -i P06280.fasta -m 0 -o blastsearch_default.out -v 700 -b 700 ./extract_ids_blast.sh blastsearch_default.out perl ../download-annotation.pl blastsearch_default_ids.txt perl ../compare_GO_terms.pl P06280 blastsearch_default_ids_GOterms.tsv perl parse_blast.pl blastsearch_default.out
The run took about 2 minutes (see section Time)
Psi-Blast
Iterations: 2 Evalue: 0.002 real 3m30.256s user 2m58.070s sys 0m13.360s Iterations: 2 Evalue: 0.000000001 real 3m8.507s user 3m5.180s sys 0m2.400s Iterations: 2 Evalue: 0.0000000001 real 3m10.271s user 3m7.620s sys 0m2.190s Iterations: 10 Evalue: 0.002 real 15m29.218s user 15m8.910s sys 0m12.730s Iterations: 10 Evalue: 0.000000001 real 16m33.748s user 16m12.500s sys 0m13.080s Iterations: 10 Evalue: 0.0000000001 real 16m20.137s user 15m55.910s sys 0m13.190s
HHblits / HHsearch
We searched the "big80" database with HHblits using the default settings and also with the maximum number of possible iterations (8) with the following commands:
time hhblits -i ../P06280.fasta -d /mnt/project/pracstrucfunc12/data/hhblits/uniprot20_current -e 0.003 -o hhblits_default.out -E 0.003 -z 700 ./extract_ids_hhblits.sh hhblits_default.out perl ../download-annotation.pl hhblits_default_ids.txt perl ../compare_GO_terms.pl P06280 hhblits_default_ids_GOterms.tsv perl parse_hhblits.pl hhblits_default.out time hhblits -i ../P06280.fasta -d /mnt/project/pracstrucfunc12/data/hhblits/uniprot20_current -e 0.003 -o hhblits_n8_neu.out -E 0.003 -n 8 -z 800 -b 800 ./extract_ids_hhblits.sh hhblits_n8_neu.out perl ../download-annotation.pl hhblits_n8_neu_ids.txt perl ../compare_GO_terms.pl P06280 hhblits_n8_neu_ids_GOterms.tsv perl parse_hhblits.pl hhblits_n8_neu.out R CMD BATCH hist_hhblits.R
The first HHblits run took about 2.5 minutes, the second one about 16 minutes (see section Time).
Time
We evaluated the time the programs ran with the command "time"
Method | Parameter | Time |
---|---|---|
Blast v = 700 | b = 700, v = 700 | 1m53.944s |
HHBlits | default | 2m19.519s |
HHBlits | n = 8 | 16m7.754s |
Comparison
>R CMD BATCH all_Evalues.R