Difference between revisions of "ARSA search protocol"
From Bioinformatikpedia
m (→PSI-Blast) |
m (→PSI-Blast) |
||
| Line 42: | Line 42: | ||
2040.800u 54.190s 36:10.13 96.5% 0+0k 19511152+12824io 2375pf+0w |
2040.800u 54.190s 36:10.13 96.5% 0+0k 19511152+12824io 2375pf+0w |
||
| − | # To evaluate the final results, we have to dissect the output into separate files: Look for the first hit, find at what line numbers it occurs and then cut the files accordingly. |
+ | # To evaluate the final results for big_80, we have to dissect the output into separate files: Look for the first hit, find at what line numbers it occurs and then cut the files accordingly. |
grep -n G3IH84 ARSA.psiBlast.j2.hDefault |
grep -n G3IH84 ARSA.psiBlast.j2.hDefault |
||
tail -n +4679 ARSA.psiBlast.j2.hDefault > ARSA.psiBlast.j2.hDefault.lastIter |
tail -n +4679 ARSA.psiBlast.j2.hDefault > ARSA.psiBlast.j2.hDefault.lastIter |
||
| − | |||
# Extract unique matches as above for Blast |
# Extract unique matches as above for Blast |
||
uniq -w 44 ARSA.psiBlast.j2.hDefault.lastIter > ARSA.psiBlast.j2.hDefault.lastIter.uniq |
uniq -w 44 ARSA.psiBlast.j2.hDefault.lastIter > ARSA.psiBlast.j2.hDefault.lastIter.uniq |
||
| + | |||
| + | # Or directly work on the run with the complete big, in order to compare to hhblits |
||
| + | blastpgp -d /mnt/project/pracstrucfunc12/data/big/big -m 8 -b 100000 -j 1 -i ARSA.fas -o ARSA.psiBlast.j2.hDefault.big -R ARSA.psiBlast.j2.hDefault.chk |
||
| + | ... |
||
| + | uniq -w 44 ARSA.psiBlast.j2.hDefault.big > ARSA.psiBlast.j2.hDefault.lastIter.big |
||
| + | |||
| + | |||
| + | # generate histogram of e-value distributions in R: (see psiBlast.evalueHist.Rscript) |
||
| + | Rscript psiBlast.evalueHist.Rscript |
||
| + | |||
| + | |||
| + | |||
| + | |||
| + | # find the pdb ids of the different structure groups in COPS: |
||
| + | head -1 ARSA.psiBlast.j10.h1e-10.pdb.ids.uniq |
||
| + | grep -i 1n2l /mnt/project/pracstrucfunc12/data/COPS/COPS-ChainHierarchy.txt |
||
| + | grep 3nkm /mnt/project/pracstrucfunc12/data/COPS/COPS-ChainHierarchy.txt | cut -c 1-4 | sort | uniq > COPS-ARSA-L30.pdbIds |
||
| + | grep 1p49 /mnt/project/pracstrucfunc12/data/COPS/COPS-ChainHierarchy.txt | cut -c 1-4 | sort | uniq > COPS-ARSA-L40.pdbIds |
||
| + | grep "1p49,A 1p49,A" /mnt/project/pracstrucfunc12/data/COPS/COPS-ChainHierarchy.txt | cut -c 1-4 | sort | uniq > COPS-ARSA-L60.pdbIds |
||
# To evaluate the true/false positives use the checkpoint files to search against pdb |
# To evaluate the true/false positives use the checkpoint files to search against pdb |
||
| Line 54: | Line 72: | ||
cut -f 2 ARSA.psiBlast.j2.hDefault.pdb | cut -c 9-14 | uniq > ARSA.psiBlast.j2.hDefault.pdb.ids.uniq |
cut -f 2 ARSA.psiBlast.j2.hDefault.pdb | cut -c 9-14 | uniq > ARSA.psiBlast.j2.hDefault.pdb.ids.uniq |
||
... |
... |
||
| + | cut -c 1-4 ARSA.psiBlast.j2.hDefault.pdb.ids.uniq | sort | uniq > ARSA.psiBlast.j2.hDefault.pdb.ids.uniq2 |
||
| + | ... |
||
| + | |||
| + | # Now we can just do a diff to find out which / how many pdb ids are missing / too much |
||
| + | |||
| + | |||
| + | # |
||
| + | |||
| + | |||
| + | |||
| + | |||
| + | # |
||
| − | # generate histogram of e-value distributions in R: |
||
| − | #colors |
||
| − | cols4 <- hcl(h = seq(30, by=360 / 4, length = 4), l = 65, alpha = 0.5) |
||
| − | # read data |
||
| − | j2_hDefault = read.table("ARSA.psiBlast.j2.hDefault.big.uniq") |
||
| − | j10_hDefault = read.table("ARSA.psiBlast.j10.hDefault.big.uniq") |
||
| − | j2_h10 = read.table("ARSA.psiBlast.j2.h1e-10.big.uniq") |
||
| − | j10_h10 = read.table("ARSA.psiBlast.j10.h1e-10.big.uniq") |
||
| − | # hist1 |
||
| − | hist(log10(j2_hDefault$V11),border=cols4[1],col=cols4[1],breaks=151,main='PSI-BLAST logarithm of E-values',panel.first = grid(),xlab='lg(E-value)',axes=FALSE,xlim=range(-150,1),ylim=range(0,1000)) |
||
| − | # hist2-4 |
||
| − | hist(log10(j10_hDefault$V11),border=cols4[2],col=cols4[2],breaks=151,add=TRUE) |
||
| − | hist(log10(j2_h10$V11),border=cols4[2],col=cols4[2],breaks=151,add=TRUE) |
||
| − | hist(log10(j10_h10$V11),border=cols4[2],col=cols4[2],breaks=151,add=TRUE) |
||
| − | #axes |
||
| − | axis(1, las=1, labels=seq(-150,0,50), at=seq(-150,0,50), tcl=0.5) |
||
| − | axis(2, las=1, labels=seq(0,1000,100), at=seq(0,1000,100), tcl=0.5) |
||
| − | #legend |
||
| − | legend(0.8,100,c('j 2, h 0.002','j 10, h 0.002','j 2, h 1e-10','j 10, h 1e-10'),col=cols4,fill=cols4,border=cols4) |
||
Revision as of 12:34, 30 April 2012
Blast
# Run blast: blastall -p blastp -d /mnt/project/pracstrucfunc12/data/big/big_80 -m 8 -b 10000 -i ARSA.fas > ARSA.blast
# find unique matches cut -f 2 ARSA.blast | uniq > ARSA.blast.matchIDs.uniq uniq -w 44 ARSA.blast > ARSA.blast.uniq
# use R to plot:
# within R: histogram of log10 Evalue for blast run
blastTable = read.table("ARSA.blast.uniq")
pdf("ARSA_blastHistogram_eVal.pdf")
logeval.hist <- hist(log10(blastTable$V11), breaks=100, plot=FALSE)
plot(logeval.hist, col="blue", main="BLAST logarithm of E-values", xlab="lg(E-Value)")
dev.off()
# within R: histogram of %identity for blast run
blastTable = read.table("ARSA.blast.uniq")
pdf("ARSA_blastHistogram_id.pdf")
id.hist <- hist(blastTable$V3, breaks=5, plot=FALSE)
plot(id.hist, col="blue", main="BLAST sequence identity distribution", xlab="percent identity")
dev.off()
# create lists of Uniprot IDs for analysis for GO annotations grep -n "1e-10" ARSA.blast.uniq tail -n +1721 ARSA.blast.uniq > ARSA.blast.uniq.worseE10 cut -f 2 ARSA.blast.uniq.worseE10 > ARSA.blast.uniq.worseE10.matchIDs cut -d "|" -f 2 ARSA.blast.uniq.worseE10.matchIDs > ARSA.blast.uniq.worseE10.matchIDs.clean
PSI-Blast
# run PSI-Blast with different combinations of parameters (generate checkpoint files for scanning PDB) > blastpgp -d /mnt/project/pracstrucfunc12/data/big/big_80 -m 8 -b 10000 -j 2 -i ARSA.fas -o ARSA.psiBlast.j2.hDefault -C ARSA.psiBlast.j2.hDefault.chk 280.470u 30.420s 5:55.56 87.4% 0+0k 18005712+1880io 792pf+0w > blastpgp -d /mnt/project/pracstrucfunc12/data/big/big_80 -m 8 -b 10000 -j 10 -i ARSA.fas -o ARSA.psiBlast.j10.hDefault -C ARSA.psiBlast.j10.hDefault.chk [blastpgp] ERROR: ncbiapi [000.000] ObjMgrNextAvailEntityID failed with idx 2048 2111.750u 44.740s 36:55.25 97.3% 0+0k 12951072+11856io 1271pf+0w > blastpgp -d /mnt/project/pracstrucfunc12/data/big/big_80 -m 8 -b 10000 -j 2 -h 1e-10 -i ARSA.fas -o ARSA.psiBlast.j2.h1e-10 -C ARSA.psiBlast.j2.h1e-10.chk 280.140u 6.580s 4:55.89 96.9% 0+0k 2536896+2000io 258pf+0w > blastpgp -d /mnt/project/pracstrucfunc12/data/big/big_80 -m 8 -b 10000 -j 10 -h 1e-10 -i ARSA.fas -o ARSA.psiBlast.j10.h1e-10 -C ARSA.psiBlast.j10.h1e-10.chk [blastpgp] ERROR: ncbiapi [000.000] ObjMgrNextAvailEntityID failed with idx 2048 2040.800u 54.190s 36:10.13 96.5% 0+0k 19511152+12824io 2375pf+0w
# To evaluate the final results for big_80, we have to dissect the output into separate files: Look for the first hit, find at what line numbers it occurs and then cut the files accordingly. grep -n G3IH84 ARSA.psiBlast.j2.hDefault tail -n +4679 ARSA.psiBlast.j2.hDefault > ARSA.psiBlast.j2.hDefault.lastIter # Extract unique matches as above for Blast uniq -w 44 ARSA.psiBlast.j2.hDefault.lastIter > ARSA.psiBlast.j2.hDefault.lastIter.uniq
# Or directly work on the run with the complete big, in order to compare to hhblits blastpgp -d /mnt/project/pracstrucfunc12/data/big/big -m 8 -b 100000 -j 1 -i ARSA.fas -o ARSA.psiBlast.j2.hDefault.big -R ARSA.psiBlast.j2.hDefault.chk ... uniq -w 44 ARSA.psiBlast.j2.hDefault.big > ARSA.psiBlast.j2.hDefault.lastIter.big
# generate histogram of e-value distributions in R: (see psiBlast.evalueHist.Rscript) Rscript psiBlast.evalueHist.Rscript
# find the pdb ids of the different structure groups in COPS: head -1 ARSA.psiBlast.j10.h1e-10.pdb.ids.uniq grep -i 1n2l /mnt/project/pracstrucfunc12/data/COPS/COPS-ChainHierarchy.txt grep 3nkm /mnt/project/pracstrucfunc12/data/COPS/COPS-ChainHierarchy.txt | cut -c 1-4 | sort | uniq > COPS-ARSA-L30.pdbIds grep 1p49 /mnt/project/pracstrucfunc12/data/COPS/COPS-ChainHierarchy.txt | cut -c 1-4 | sort | uniq > COPS-ARSA-L40.pdbIds grep "1p49,A 1p49,A" /mnt/project/pracstrucfunc12/data/COPS/COPS-ChainHierarchy.txt | cut -c 1-4 | sort | uniq > COPS-ARSA-L60.pdbIds
# To evaluate the true/false positives use the checkpoint files to search against pdb blastpgp -d pdb_seqres.fasta -m 8 -b 10000 -j 1 -i ARSA.fas -o ARSA.psiBlast.j2.hDefault.pdb -R ARSA.psiBlast.j2.hDefault.chk ... cut -f 2 ARSA.psiBlast.j2.hDefault.pdb | cut -c 9-14 | uniq > ARSA.psiBlast.j2.hDefault.pdb.ids.uniq ... cut -c 1-4 ARSA.psiBlast.j2.hDefault.pdb.ids.uniq | sort | uniq > ARSA.psiBlast.j2.hDefault.pdb.ids.uniq2 ...
# Now we can just do a diff to find out which / how many pdb ids are missing / too much
#
#
