ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/owl/trunk/modelit/scripts/align_seq2prof.sh
Revision: 1299
Committed: Tue Jan 11 17:38:33 2011 UTC (10 years, 5 months ago) by hstehr
File size: 2014 byte(s)
Log Message:
adding Model-It server files to repository
Line User Rev File contents
1 hstehr 1299 #!/bin/sh
2     # Convenience script for aligning a sequence against a profile using T-Coffee
3     # Henning Stehr, h.stehr@molgen.mpg.de
4     # 10/Apr/2008
5    
6     if [ -z "$3" ]
7     then
8     echo "Usage: $0 <seq.fasta> <alignm.fasta> <outfile.fasta> <logfile>"
9     exit 1
10     fi
11    
12     seq=$1
13     prof=$2
14     outfile=$3
15     logfile=$4
16    
17     scriptsdir=/project/StruPPi/CASP8/scripts
18     psiblastFile=`dirname $prof`/`basename $prof .templates_aln.fasta`.pdb-psiblast.classic.out
19    
20     numTempl=`grep -c "^>" $prof`
21     if [ "$numTempl" -eq "1" ]
22     then
23     echo "Single template, using alignment from psiblast from file $psiblastFile"
24     if [ -e "$psiblastFile" ]
25     then
26     pdbId=`cat $prof | grep "^>" | sed "s/>//"`
27     pdbCode=`echo $pdbId | cut -c1-4`
28     pdbChainCode=`echo $pdbId | cut -c5`
29     perl $scriptsdir/parse_blast.pl -i $seq -a $psiblastFile -p $pdbCode -c $pdbChainCode > $outfile
30     if [ "$?" -ne "0" ]
31     then
32     echo "ERROR in parsing the psi-blast output, no target to template alignment!"
33     exit 1
34     fi
35     else
36     echo "ERROR: couldn't find $psiblastFile, no target to template alignment!"
37     exit 1
38     fi
39     exit 0
40     fi
41    
42    
43     # t-coffee seems to have a bug in outputting directly to fasta_aln format: it will add
44     # one extra character to all sequences except for the target. This is why we first write
45     # to clustalw format and convert to fasta_aln
46    
47     # yet another bug in t-coffee. This one is an annoyance rather than a bug:
48     # it takes the file name of the profile file as the sequence name, but maximum sequence name length is limited to 100,
49     # so it fails if path of the file name is really long. To avoid this we first copy the $prof to /tmp
50    
51     cp $prof /tmp
52     tmpprof=/tmp/`basename $prof`
53    
54     tmpoutfile=/tmp/`basename $outfile`.tmp
55     tmpoutfile2=/tmp/`basename $outfile`.tmp2
56    
57     t_coffee $seq -profile $tmpprof -profile_comparison=full50 -output=clustalw -outfile=$tmpoutfile -quiet $logfile
58    
59     t_coffee -other_pg seq_reformat -in $tmpoutfile -output fasta_aln > $tmpoutfile2
60    
61     cat $tmpoutfile2 | perl -ne 'if(/^>/){print}else{print uc($_)}' > $outfile
62    
63     rm -f $tmpoutfile $tmpoutfile2 $tmpprof

Properties

Name Value
svn:executable *