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 (8 years, 8 months ago) by hstehr
File size: 2014 byte(s)
Log Message:
adding Model-It server files to repository
Line File contents
1 #!/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 *