Hi all, I need a trained python eye for this :)
I need to remove 100's of genes from a proteome file contains 1000s genes. Obviously I do not want to do it manually. I have pulled the python code pasted below from somewhere, which is a few years old. It is supposed to do what I want, but it does not. It just copies all the files from the original file to the output file, ignoring the remove.file. This code requires 3 files which I supplied. File 1; "123.fasta" - the file with my original unedited proteome, file 2; "remove.txt" - the file with the list of gene ID's to be removed. File 3. "new.fasta" - the output file with the edited proteome minus the genes listed in the remove.txt file. Ideally, I would like the code to identify the genes in "123.fasta" by the fasta format sequence ID (eg. >sequence1, >sequence2 etc).
This is the code:
import Bio
from Bio import SeqIO
import sys
fasta_file = ("123.fasta")
remove_file = ("remove.txt")
result_file = ("new.fasta")
remove = set (">")
with open(remove_file) as f:
for line in f:
line = line.strip()
if line != "":
remove.add(line)
fasta_sequences = SeqIO.parse(open(fasta_file), "fasta")
with open(result_file, "w") as f:
for seq in fasta_sequences:
nam = str()
nam = nam.stripseq.id)
nuc = str(seq.seq)
SeqIO.write([seq], f, "fasta")
As I said, no matter what I tweak, it just copies and pastes all of the 123.fasta file into the output file, no deletions. Any of the python people see what may be the problem? I am not a trained python operator , just using it for my work.
•
written
by
Bacteria_forever •
On closer inspection there is a lot wrong with the code in your opening post, so its easier just to start from scratch:
from Bio import SeqIO
import sys
with open(sys.argv[1], 'r') as rfh:
to_remove = [line.strip() for line in rfh]
to_keep = []
for record in SeqIO.parse(sys.argv[2], "fasta"):
if record.id not in to_remove:
to_keep.append(record)
SeqIO.write(to_keep, sys.argv[3], "fasta")
python scriptname.py list_to_remove.txt input.fasta output.fastaNot the worlds most efficient script, but works well enough.
If your IDs aren't simple strings that correspond directly, (e.g. they contain the > character etc) you'll need to do a bit more work to generate the to_remove list.
Hmmmm, thank you, But your amendment returns this:
Traceback (most recent call last): File "biostarans.py", line 19, in <module> nam = nam.stripseq.id() AttributeError: 'str' object has no attribute 'stripseq' ??
Yeah seems weird, the code is not coming out right when I submit in message above, hence your amendment didn't work. Sorry. There are brackets missing, which will not reproduce even I post them here now. But still, when I fix it in my terminal with your addition it still just straight up copying and pasting without editing out the remove list. Damn frustrating!
Log in to answer this question.
Since this question is not about python code you wrote consider
faSomeRecordsutility from Jim Kent (LINK). After downloading the file add execute permissions (chmod u+x faSomeRecords). Use as follows