6
$\begingroup$

I have the following mwe for filtering a Swissprot file based on a certain feature, in this case, transmembrane regions.

from Bio import SeqIO
records=[]
for record in SeqIO.parse("Input.txt", "swiss"):
    transmembrane_protein=False
    print record.id
    for i, feature in enumerate(record.features):
        if feature.type == "TRANSMEM":
            transmembrane_protein=True
    if transmembrane_protein==True:
        records.append(record)
SeqIO.write(records, "Output.txt", "swiss")

The script works when SeqIO.write(records, "Output.txt", "swiss") becomes SeqIO.write(records, "Output.txt", "fasta")

However, such a method is not yet supported.

ValueError: Reading format 'swiss' is supported, but not writing

From the docs, I see that writing is not supported for swiss:

Note that while Bio.SeqIO can read all the above file formats, it cannot write to all of them.

Is there any unofficial way of using biopython/python to write swissprot files from parsed swissprot files?

$\endgroup$
  • $\begingroup$ I read that "SeqIO cannot output the file format (e.g. the plain text SwissProt format, "swiss") ". But you could paste the proteins without using Biopython once you know the id of the features you want $\endgroup$ – llrs Jun 29 '17 at 9:12
  • $\begingroup$ It must be in some of the docs somewhere I'm sure, I just haven't seen it explicitly stated, but it is pointed out in the error handling so that's good enough for me to know there is no easy way. I'm dealing with several proteomes. Manual pasting wouldn't be fun ;) I could really do with a programmatic solution. $\endgroup$ – James Jun 29 '17 at 9:20
  • $\begingroup$ You could use grep, or awk to select the records and output them in another file using a pipe. $\endgroup$ – llrs Jun 29 '17 at 9:22
  • $\begingroup$ @Llopis That would certainly be one way forward. However I'd like to keep the script in a single python file because it's cross platform and transportable. $\endgroup$ – James Jun 29 '17 at 9:59
6
$\begingroup$

Using SeqIO.index rather than SeqIO.parse lets you read all the records into a dict, from which you can then extract the raw entry:

from Bio import SeqIO

record_dict = SeqIO.index('Input.txt', 'swiss')
for key in record_dict:
    print(record_dict.get_raw(key).decode())

Now you should be able to apply your test for a transmembrane protein to each entry, and write out only the ones you want to keep.

from Bio import SeqIO

record_dict = SeqIO.index('Input.txt', 'swiss')
out = open('Output.txt', 'wb')
for key in record_dict:
    record = record_dict[key]
    for i, feature in enumerate(record.features):
    if feature.type == "TRANSMEM":
        out.write(record_dict.get_raw(key).decode())
| improve this answer | |
$\endgroup$
  • $\begingroup$ Perfect solution, worked like a charm. $\endgroup$ – James Jun 30 '17 at 4:12
2
$\begingroup$

I don't know how to do this with biopython, but it's simple enough in Perl using Swisskinfe:

#!/usr/bin/env perl
use strict;
use warnings;
use SWISS::Entry;

my $flat=$ARGV[0]||die("Need a flat file\n");

## Change the line termination string so we read an 
## entire entry at a time
local $/ = "\n//\n";

my $fh;
## Deal with compressed or uncompressed files
if($flat=~/\.gz$/){
    open($fh,"zcat $flat |")|| die("cannot open flat file $flat : $! \n");
}
else{
    open($fh,"$flat")|| die("cannot open flat file $flat : $!: $@\n");
}
## Parse the file. 
entry: while (my $current_entry = <$fh>) {
    ## Read the current entry as an "entry" object
    my $entry = SWISS::Entry->fromText($current_entry);
    ## Iterate over the entry's features
    foreach my $feature ($entry->FTs->elements()) {
        ## If any of the features are 'TRANSMEM', print the entry
        if ($feature->[0] eq 'TRANSMEM') {
            print "$current_entry\n";
            ## Skip to the next entry
            next entry;
        }
    }
}

You need to install the Swissknife module first:

  1. Download Swissknife_1.73.tar.gz from here (SourceForge link).

  2. Install

    tar xvzf Swissknife.tar.gz
    cd Swissknife_1.73/
    perl Makefile.PL
    sudo make install
    

Now, make the script above executable and run with:

foo.pl file.flat
| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.