I have a Swiss-Prot database file that contains several Swiss-Prot Files.
They are copied and pasted underneath each other.
Therefore there is one Swiss-Prot entry after another listed in the same file.
I want to write the ID into another file as the header. Immediately underneath, I want to write the amino acid sequence.
So far I can only read one single Swiss-Prot file and get as output 1ID and 1 amino acid sequence. In other words, I have managed to print out the ID header first and the amino acid sequence second .
How can this code work to read multiple Swiss-Prot file entries from one single file?
How do I do this sequentially for every ID and amino acid sequence from each Swiss-Prot entry listed in the file?
bright_cyan = "\033[0;96m"
bright_yellow = "\033[0;33m"
bright_green = "\033[0;32m"
reset = "\033[0m"
#--------------------------------------------------------------------
import sys
import re
#--------------------------------------------------------------------
def read_data(SPROT_FILE):
''' This function is what is is aint it '''
flag = ''
try:
DNAfile = open(SPROT_FILE , 'r')
except IOError as error:
print(bright_cyan + "double check and see if you entered the correct filename :> ", str(error))
sys.exit(1)
# create a FASTA file to copy the information to and write.
new_outfile = open("first.fsa", 'w')
amino_acid_sequence = ''
for line in DNAfile:
#print(line, end = '')
if re.match(r'ID', line):
ID = line[5:20]
# Stateful Parsing of the amino acid sequence.
if re.match(r'//', line):
flag = False
if flag:
amino_acid_sequence += line
if re.match(r'SQ', line):
flag = True
# Find the modified amino acid residue.
if re.match(r'FT MOD_RES', line):
FT = line
position_switch = ','.join(re.findall(r'\d+',FT))
header_line = '>'+ID.strip()+" phospho:"+position_switch
print(header_line)
#print('>'+ID.strip()+" phospho:"+position_switch, file = new_outfile)
# Print each amino acid sequence outside of the loop.
amino_acid_sequence = amino_acid_sequence.replace(' ', '')
print(amino_acid_sequence)
# Write the amino acid sequence to the file.
print(amino_acid_sequence, file = new_outfile)
DNAfile.close()
new_outfile.close()
# Not sure about this part...
files = input(bright_yellow + 'Type possibly filenames :> ').split()
for filename in files:
read_data(filename)
I hope the question is clear.
Would be great it if you could offer some help.
Thanks in advance
2 answers
If you download older HMMer (say, 2.3.2 version), there is a program called sreformat that will directly convert this format to fasta.
The 2.3.2 version of HMMer doesn't have the ehmmcalibrate.c file - I just checked. That would indicate that you are working with a different version.
Separately, HMMer is one of the best-behaved programs I have ever encountered in terms of compiling. I have tried all the major versions of it, and not once did I have to do anything more than a simple:
./configure ; make
I downloaded and recompiled 2.3.2 as I was writing this, and it took less than 30 seconds. sreformat is part of HMMer's squid library and it will be in the corresponding directory. If you have problems compiling HMMer, chances are that something standard is missing from your system.
You say you aren't interested in this route in part because it isn't 'pythonic'. I'd argue not repeating making software that already exists is very Pythonic as it adheres strong 'DRY' principles.
In fact, most command line software can be run from within a Python script if you are looking to use Python as the backbone of a workflow; often os.system(<command_here>) is easier than subprocess. I demonstrate this near the end of the notebook I'm going to suggest checking out here.
Be that as it may, others may want to follow this route...
This is a case where having another system to do you work on can be handy. It's especially nice because it makes it so there is zero chance of messing up your system trying to install old software. Plus, it makes it more reproducible by eliminating the 'it works on my machine'-issue.
Put https://mybinder.org/v2/gh/jupyterlab/jupyterlab-demo/HEAD?urlpath=lab/tree/demo in your browser's address bar and hit return, or click that URL to launch a JupyterLab session.
We don't need a special environment here, and so it is the same as available for JupyterLab from Try Jupyter.
When the session comes up, open a new Jupyter notebook and paste in the following in a cell and run it with Shift-Enter:
!curl -OL https://gist.githubusercontent.com/fomightez/cb3a7f13a9b1ff74f55ac23835eb28a5/raw/56b098763f50be78f718642ad7e1a956a9df99e3/Guide_to_using_advice_posted_in_Biostars_answer_9500884.ipynb
That will get a notebook Guide_to_using_advice_posted_in_Biostars_answer_9500884.ipynb. When that shows up after a few seconds, in the file navigation panel on the left, double-click to open it and then execute the entire thing by selecting Run > Run All Cells from the toolbar menu. Alternatively, you can just step through running the cells with Shift-Enter to follow along. The 'Preparation' section installs the software and sets up to process your example data. You won't be able to run the lower cells until you install the software on the new session.
The session is ephemeral, and so if you use to do the conversion, make sure to grab anything useful.
All but the use of Python at the end of the linked Jupyter notebook could be done on the command line with the same commands without the exclamation signs or percent symbols. The notebook just makes it easier to share the commands and the result.
Direct link to static version of that notebook:
here in nbviewer which presently renders gists better than github on my system
You say you aren't interested in this route in part because it isn't 'pythonic'. I'd argue not repeating making software that already exists is very Pythonic as it adheres strong 'DRY' principles.
Amen !
If the error occurs at the very end, chances are that sreformat may have been compiled successfully before the compilation stopped. It should be in the squid subdirectory. If you don't have a squid subdirectory, you are not compiling the 2.3.2 version.
Hi, is this what you want to do ?
import os
import sys
import re
def main(input_fname: str, output_fname: str) -> None:
with open(output_fname, 'w') as f_out:
with open(input_fname, 'r') as f_in:
for record in re.split('//', f_in.read())[:-1]:
record_id = re.split('\s+', record[record.index('ID'):])[1]
sequence = ''.join(record[record.index('SQ'):].split('\n')[1:]).replace(' ','')
f_out.write(f'>{record_id}\n{sequence}\n')
if __name__ == '__main__':
main(input_fname=sys.argv[1], output_fname=sys.argv[2])
•
written
by
hugo.avila
I would like all of the amino acid sequences to be directly underneath the ID line, so that the output file is:
ID LINE
- AA sequence
ID LINE
- AA sequence
ID LINE
- AA sequence and so on ....
There is a Value Error of Substring not found in the last line that I am trying to iron out.
enter code here
#!/usr/bin/env python3
import os
import sys
import re.
# Get input file name
if len(sys.argv) == 3:
input_fname = sys.argv[1]
output_fname = sys.argv[2]
else:
sys.stderr.write("Usage: pythonfile.py <input filename> <output filename> \n")
sys.exit(1)
def main(input_fname: str, output_fname: str) -> None:
with open(output_fname, 'w') as f_out:
with open(input_fname, 'r') as f_in:
# Can you write a comment here to explain this line below?
for record in re.split('//', f_in.read())[:-1]:
# What is the record[record.index ?
record_id = re.split('\s+', record[record.index('ID'):])[1]
# Yes, except this only works for one ID and one sequence, as noted in my initial post.
# error here >> sequence = ''.join(record[record.index('SQ'): >>>
split the file from SQ to // so that each AA sequence is seperate i.e. until the next '//' <<<.
].split('\n')[1:]).replace(' ','')
f_out.write(f'>{record_id}\n{sequence}\n')
if __name__ == '__main__':
main(input_fname=sys.argv[1], output_fname=sys.argv[2])
import sys
import re
def main(input_fname: str, output_fname: str) -> None:
with open(output_fname, 'w') as f_out:
with open(input_fname, 'r') as f_in:
# The input file is separated with '//' so if we split the file by these caracters
# it is possible to get a list of records that can by looped:
# "record1//record2//record//" --split('//')--> ["record1", "record2", "record3", ""].
# As you can see above the last item of the splitted string is a "" (empty string)
# so we need to ignore it
# li like this: ["record1", "record2", "record3", ""][:-1] -> ["record1", "record2", "record3"].
for record in re.split('//', f_in.read())[:-1]:
# I cant see way the first code only returned the first id, i did run it and it worked for the sample.
# I think that must be some format error in one of the sequences. This try block is kind of
# ugly but it wil go to the end of your file, write the output and print the unformated records (if any).
try:
# split a record into a list of lines and get only the one that starts with 'ID'.
# I think that maybe u only want the whole line so i did not pull only the id.
id_line = list(filter(lambda x: x.startswith('ID'), record.split('\n')))[0]
sequence = ''.join(record[record.index('SQ'):].split('\n')[1:]).replace(' ','')
f_out.write(f'>{id_line}\n{sequence}\n')
except Exception as e:
print (f'Unformated record\n{record}')
pass
if __name__ == '__main__':
main(input_fname=sys.argv[1], output_fname=sys.argv[2])
Hi, could you post a sample of the input file where the code fails ?
To reliably find the SQ line, you need to look for lines starting with "SQ" and then followed by 3 spaces.
Also, I strongly recommend using the primary accession number (first identifier on the first line starting with AC), and not the ID line. The ID line contains an entry name/mnemonic which cannot be guaranteed to remain stable (https://www.uniprot.org/help/entry_name vs https://www.uniprot.org/help/accession_numbers).
If you have a list of accession numbers, you can post them to https://www.uniprot.org/id-mapping and download in FASTA format.
If not, you could download the Swissknife package from https://swissknife.sourceforge.net/ and run this script over your data:
use strict;
use IO::File;
use SWISS::Entry;
my $inputfile = @ARGV[0];
my $fh = new IO::File $inputfile or
die "Cannot open input file $inputfile: $!";
$/ = "\n\/\/";
while(<$fh>) {
s/\r//g;
(my $entry_txt = $_) =~ s/^\s+//;
next unless $entry_txt;
$entry_txt .= "\n";
my $entry = SWISS::Entry->fromText( $entry_txt );
print $entry->toFasta();
}
Hey, cool implementation !
I added a demo of both 'hugo.avila''s scripts along with the use of sreformat to the Jupyter notebook linked in my post above that was meant to install & demonstrate sreformat.
Hey @Wayne , really nice job with the notebook ! I didn't know about binder, i'll be using it for now on, it seems very practical.
Log in to answer this question.
-
Eukaryotic gene prediction using GeneMark-ES
written by lakhujanivijay
-
Python script error
written by Abdul Rafay Khan
-
GTF to FASTA including stop codon
written by athey.johnc •
-
Beginner in Python- translating DNA given in GenBank file format into its six reading frames as out…
written by oki4 •
-
Retrieving Swissprot FT details for a particular amino acid position using Python
written by a.gardner •
-
Script For Parsing A Biological Sequence From A Public Database In Python
written by Spyros •
Can you give an example of what you mean by swiss-prot file? I think you're describing a fasta file with amino acids. In that case use BioPython to parse the fasta file.
Yes. Here is an example of the file. It is a few thousand lines long so I won't put the whole thing.
Hopefully that is clearer now.
Best