Lab Journal of Task 3 (MSUD)

For task 3 we have used the reference sequence of BCKDHA and other given example proteins.

Secondary structure

Prediction and assignment

• PSSMs were created with Psi-Blast:

blastpgp -d /mnt/project/pracstrucfunc13/data/big/big_80 -i P10775.fasta -j 2 -h 10e-10 -Q P10775_big80.blastPsiMat

blastpgp -d /mnt/project/pracstrucfunc13/data/swissprot/uniprot_sprot -i P10775.fasta -j 2 -h 10e-10 -Q P10775_SwissProt.blastPsiMat

• ReProf was run for P10775 with a simple fasta file and with a PSSM (generated with big_80 and SwissProt, respectively) as input:

reprof -i P10775.fasta

reprof -i P10775_big80.blastPsiMat

reprof -i P10775_SwissProt.blastPsiMat

• For the remaining proteins (P12694, Q08209 and Q9X0E6) ReProf was run with a big_80 PSSM analogously as shown above.

• PsiPred and DSSP were run on the following servers: PsiPred, DSSP_server.

• The PDB files used as input for DSSP are located at /mnt/home/student/schillerl/MasterPractical/task3/pdb_structures/. If there were more than one PDB structures available those which the highest coverage over the sequence and with the highest resolution were taken preferentially. These structures were used: 2BFD (P12694), 2BNH (P10775), 1AUI (Q08209) and 1KR4 (Q9X0E6).

• To parse the output of ReProf and PsiPred, we used Sonja's script. For parsing the DSSP output, we used the following Python script (which is located at /mnt/home/student/schillerl/MasterPractical/task3/dssp/parse_dssp_output.py):

<source lang=python> Parse DSSP output to a sequence of H, E and L (secondary structure). If there is more than one sequence in the DSSP output, only the first sequence is parsed.

Call with python parse_dssp_putput.py <dssp output file> <sequence length> <output file>

@author: Laura Schiller

import sys

seq_len = int(sys.argv[2]) dssp_file = open(sys.argv[1]) out_file = open(sys.argv[3], "w")

1. translation of DSSP secondary structure to H, E and L

secstr = { 'H': "H",

           'G': "H", 
           'I': "H", 
           'E': "E", 
           'B': "E", 
           'S': "L", 
           'T': "L", 
           ' ': "L" }

line = dssp_file.readline() while not line.startswith(" # RESIDUE"):

   line = dssp_file.readline()

line = dssp_file.readline()

current_position = 1 while (line):

   try:
       position = int(line[7:10])
   except ValueError:
       if(line[13:15] == "!*"): # new amino acid chain
           break
       else:
           line = dssp_file.readline()
           continue
   structure = secstr[line[16:17]]
   if(current_position != position):
       for i in range(position - current_position):
           out_file.write("-") # no DSSP output for that residue
           current_position += 1
   out_file.write(structure)
   current_position += 1    
   line = dssp_file.readline()

for i in range(seq_len - current_position + 1):

   out_file.write("-")

dssp_file.close() out_file.close() </source>

Note: for the proteins P12694 and Q9X0E6 the residue numbering in the DSSP output did not agree with the fasta sequences, so the output of the above script for these proteins was manually modified to fit the fasta sequences (for P12694 the beginning of the sequence was missing in the DSSP output, so there was added sufficiently many '-' for this part; for Q9X0E6 there were some additional residues in the beginning and end in the DSSP output, which are missing in the fasta sequence, so these parts were deleted). The reason for these discrepancies may be, that the PDB files that were used as input for DSSP don't exactly match the fasta sequences from UniProt.

Evaluation of prediction approaches

• The ReProf predictions for P10775 (derived with different approaches, see above) were compared with the DSSP assignment by using the following Python script (located at /mnt/home/student/schillerl/MasterPractical/task3/evaluate_secstr_reprof.py), which calculates the recall, precision and f-measure of the predictions. Positions that lack a DSSP assignment (parsed as '-' by the above script) were ignored for the calculation.

Recall and Precision are defined as follows:

• recall = TP / (TP + FN)

• precision = TP / (TP + FP)

• f-measure = 2 * recall * precision / (recall + precision)

where TP means true positive, FP false positive and FN false negative.

<source lang="python"> dssp_file = open("./dssp/P10775_secstr.txt") dssp = dssp_file.readline() dssp_file.close()

for reprof_run in ["./reprof/P10775_secstr.txt", "./reprof/P10775_big80_secstr.txt", "./reprof/P10775_SwissProt_secstr.txt"]: reprof_file = open(reprof_run) reprof = reprof_file.readline() reprof_file.close()

assert len(dssp) == len(reprof)

sum1 = {'E': 0, 'H': 0, 'L': 0, 'all': 0} sum2 = {'E': 0, 'H': 0, 'L': 0, 'all': 0} found = {'E': 0, 'H': 0, 'L': 0, 'all': 0} right = {'E': 0, 'H': 0, 'L': 0, 'all': 0}

for i in range(len(dssp)): for secstr in ['E', 'H', 'L']: if dssp[i] == secstr: sum1[secstr] += 1 if reprof[i] == secstr: found[secstr] += 1 if reprof[i] == secstr: sum2[secstr] += 1 if dssp[i] == secstr: right[secstr] += 1

for sum in [sum1, sum2, found, right]: sum['all'] = sum['E'] + sum['H'] + sum['L']

recall = {'E': 0.0, 'H': 0.0, 'L': 0.0, 'all': 0} precision = {'E': 0.0, 'H': 0.0, 'L': 0.0, 'all': 0}

print "-----------" print "%s:" % reprof_run print "-----------"

for secstr in ['H', 'E', 'L', 'all']: recall[secstr] = (float(found[secstr]) / sum1[secstr]) print "Recall for %s: %f" % (secstr, recall[secstr]) precision[secstr] = (float(right[secstr]) / sum2[secstr]) print "Precision for %s: %f" % (secstr, precision[secstr]) print "F-measure for %s: %f" % (secstr, (2 * precision[secstr] * recall[secstr] / (precision[secstr] + recall[secstr]))) </source>

Comparison of ReProf to PsiPred and DSSP

The Reprof results for the four proteins were compared with the PsiPred predictions and DSSP assignments with the following Python script (located at/mnt/home/student/schillerl/MasterPractical/task3/cmp_secstr_reprof_psipred_dssp.py), which calculates the percent identities (number of matching residues devided by all residues; for single secondary structure elements the matching number was divided by the maximum number of appearances of this element in one of the results; positions with no DSSP assignment were excluded from the comparison with DSSP):

<source lang=python> for protein in ["P12694", "P10775", "Q08209", "Q9X0E6"]:

reprof_file = open("./reprof/" + protein + "_big80_secstr.txt") reprof = reprof_file.readline() reprof_file.close()

psipred_file = open("./psipred/" + protein + "_secstr.txt") psipred = psipred_file.readline() psipred_file.close()

dssp_file = open("./dssp/" + protein + "_secstr.txt") dssp = dssp_file.readline() dssp_file.close()

assert len(reprof) == len(dssp) == len(psipred)

ident_psipred = {'E': 0, 'H': 0, 'L': 0} ident_dssp = {'E': 0, 'H': 0, 'L': 0}

for i in range(len(reprof)): if reprof[i] == psipred[i]: if reprof[i] == "E": ident_psipred['E'] += 1 elif reprof[i] == "H": ident_psipred['H'] += 1 elif reprof[i] == "L": ident_psipred['L'] += 1 if reprof[i] == dssp[i]: if reprof[i] == "E": ident_dssp['E'] += 1 elif reprof[i] == "H": ident_dssp['H'] += 1 elif reprof[i] == "L": ident_dssp['L'] += 1

print "-----------" print "%s:" % protein print "-----------" print "percent identity to psipred:" print "H %.3f" % (float(ident_psipred['H']) / max(sum(c == 'H' for c in reprof), sum(c == 'H' for c in psipred))) print "E %.3f" % (float(ident_psipred['E']) / max(sum(c == 'E' for c in reprof), sum(c == 'E' for c in psipred))) print "L %.3f" % (float(ident_psipred['L']) / max(sum(c == 'L' for c in reprof), sum(c == 'L' for c in psipred))) print "all %.3f" % (float(ident_psipred['H'] + ident_psipred['E'] + ident_psipred['L']) / len(reprof)) print "percent identity to dssp:" print "H %.3f" % (float(ident_dssp['H']) / max(sum((reprof[i] == 'H') and (dssp[i] != '-') for i in range(len(reprof))), sum(c == 'H' for c in dssp))) print "E %.3f" % (float(ident_dssp['E']) / max(sum((reprof[i] == 'E') and (dssp[i] != '-') for i in range(len(reprof))), sum(c == 'E' for c in dssp))) print "L %.3f" % (float(ident_dssp['L']) / max(sum((reprof[i] == 'L') and (dssp[i] != '-') for i in range(len(reprof))), sum(c == 'L' for c in dssp))) print "all %.3f" % (float(ident_dssp['H'] + ident_dssp['E'] + ident_dssp['L']) / sum(c != '-' for c in dssp)) </source>

Disordered protein

IUPred

• Predictions were performed through the web server of IUPred. Graphical profiles of the results were downloaded.
• Output of IUPred are stored in the directory /mnt/home/student/weish/master-practical-2013/task03/02-disordered-protein/iupred
• We have also performed the prediction from command-line, following is the bash script:

<source lang="bash">

1. !/bin/sh -e

INPUT=$HOME/master-practical-2013/task03 OUTPUT=$HOME/master-practical-2013/task03/02-disordered-protein/iupred PARAMS="long short glob"

if [ ! -d $OUTPUT ]; then

       mkdir $OUTPUT

fi

for seq in $INPUT/*.fasta do

       filename=`basename $seq`
       for param in $PARAMS
       do
               iupred $seq $param > $OUTPUT/iupred_${filename}_$param.tsv
       done

done </source>

MetaDisorder(MD)

• As the man page of metadisorder describes, the prediction of disordered region is based on the results of other programs such as NORSnet, PROFbval etc. Rather than directly call metadisorder we have used the wrapper program predictprotein as is described on the exercise page.
• Comparison to DisProt database:
  • 3 sequences do not have entries in DisProt
  • Local alignments were performed to compare the assignment of disordered regions between DisProt and metadisorder.

Following script was called for the task: <source lang="bash">

1. !/bin/sh -e

INPUT=$HOME/master-practical-2013/task03 OUTPUT=$HOME/master-practical-2013/task03/metadisorder EXE=predictprotein

1. make output directory

if [ ! -d $OUTPUT ]; then

       mkdir $OUTPUT

fi

1. call metadisorder for all query sequences

for seq in $INPUT/*.fasta do

       filename=`basename $seq`
       $EXE --seqfile $seq --target metadisorder -p metadisorder_$filename \
       -o $OUTPUT

done echo Done! </source>

Transmembrane helices

Polyphobius

• Preparation of input data of Polyphobius:
  • Pipeline for data generation:
    • search for homologous sequences (blastget)
    • moving query sequence to top of sequence list
    • perform multiple sequence alignment (kalign)
    • run polyphobius (jphobius)
• For each of the 4 sequences, we have also used the web-tool of Polyphobius to perform predictions and get resulted polts.

Following Perl script was used for the automatic running of Polyphobius: polyphobius.pl

For each given example sequences and the reference sequence of BCKDHA, following shell commands were used: <source lang="bash"> cd /mnt/home/student/weish/master-practical-2013/task03/03-transmembrane-helices for q in seqs/*.fasta; do

   ../polyphobius.pl -q $q -o polyphobius

done </source>

Signal peptides

SignalP version 3.0 was run as follows:

signalp -t euk P02768.fasta > signalp/P02768.txt

signalp -t euk P47863.fasta > signalp/P47863.txt

signalp -t euk P11279.fasta > signalp/P11279.txt

Additionally version 4.1 was run on the SignalP server with default options.