Compare fasta file with txt file full of sequence id

user_14529630 0
Viewed 92

I need help cause I´m stuck. I have one txt file with sequence ids, it looks like this -->

tr|K9RTD0|K9RTD0_SYNP3

tr|K9RSV3|K9RSV3_SYNP3

tr|K9RRE8|K9RRE8_SYNP3

tr|K9RMU9|K9RMU9_SYNP3

Then I have a typical fasta file.

>sp|P00115|CYC6_SYNP3 Cytochrome c6 OS=Synechococcus sp. (strain ATCC 27167 / PCC 6312) OX=195253 GN=petJ PE=1 SV=2
MKTLLTILALTLVTLTTWLSTPAFAADIADGAKVFSANCAACHMGGGNVVMANKTLKKEA
LEQFGMNSADAIMYQVQNGKNAMPAFGGRLSEAQIENVAAYVLDQSSKNWAG
>tr|K9RTH7|K9RTH7_SYNP3 N-acyl-D-glucosamine 2-epimerase OS=Synechococcus sp. (strain ATCC 27167 / PCC 6312) OX=195253 GN=Syn6312_2130 PE=4 SV=1
MAPQINFPFSDLIAGYVTSYDTETDIFGLKTSDGREFPVKLSPMAYAKVIQNFDEGYPDA
TSTMRAWLTPGRFLFVYGVFYPDTDVFDAKQVVFAGKKEDDYVFEKQDWWIQQINALGKF
YVKAQFGQEEIDYRNYRTDLSVSGERSSVKFRQETDTISRLVYGFATAFMMTGDEVFLEA
AEKGTEYLRDHMRFVDRDEDIIYWYHGIDVQGEKELKIFASEFGDDYDAIPAYEQIYALA
GPIQTYRCTGDPRILSDAEQTIKLFDKFFLDQSEYGGYFSHIDPLMLDPRSDSLGRNKAR
KNWNSVGDHAPAYLINLWLATGEQKYADMLEYTFDTIEKYFPDYENSPFVQERFYEDWSH
DTTWGWQQNRAVVGHNLKIAWNLMRMQSLKPKEQYVGLAQKIADLMPSVGSDQQRGGWSD
TVERLLTNNSKFHQFVWHDRKAWWQQEQAILAYLILGGILEHDDYHRLGREAAAFYNAWF
LDLEDGGVYFNVLANGISYLARGNERAKGSHSMSGYHSFELCYLAAVYTNFLITKHPMDF
YFKPLPNGFPDRILRVSPDILPPGSILLESVEIDGKAYTDFDSQALTVKLPETKERVKVK
VRLAPKS
>tr|K9RXQ9|K9RXQ9_SYNP3 Uncharacterized protein OS=Synechococcus sp. (strain ATCC 27167 / PCC 6312) OX=195253 GN=Syn6312_3008 PE=4 SV=1
MKVEILKKRLNKECPMTTTRMPEDVIQELKQIASLLVFWGYQPLIGADIGQGLRTDLEQL
EDDKVSALVASLKRHRVSDEVLQTALMETTIN

I need to compare this two files, and find description of sequence based on the id and print it. My code:

from Bio import SeqIO
from Bio.SeqRecord import SeqRecord
import sys

p = "proteome.fasta"
file = "reference.txt"
out = "jopik.txt"


with open(out, "w") as o:
    sys.stdout = o
    for seq_record in SeqIO.parse(open(p, mode = "r"),"fasta"):
        seq_record.description=' '.join(seq_record.description.split()[1:])
        with open(file,"r") as f:
            line = f.readlines()
            print(line)
            if (seq_record.id == line):
                    i = seq_record.description
                    print(i)
biopython python
1 Answers

You're just missing some sort of loop for x in y: additionally, file handlers are iterable in Python (iterating by-lines for non-binary mode), which will save you from having to load the entire file into memory before starting iteration (as .readlines() does)

# load first file and create a helpful structure
compare_dict = {}
with open("reference.txt") as fh:
    for line in fh:
        if line:  # throw out empty lines, could do a stricter compare
            compare_dict[line.strip()] = None

# form a tuple of possible prefixes
compare_tuple = tuple(">" + a for a in compare_dict.keys())

with open("proteome.fasta") as fh:
    for line_no, line in enumerate(fh, 1):  # lines start at 1, not 0
        if line.startswith(compare_tuple)
            key, value = line.split(" ", 1)
            key = key[1:]  # strip ">" from prefix
            compare_dict[key] = value
            print("found {} on L{}: {}".format(key, line_no, value))

# optionally display keys which were not in your .fasta file
for key, value in compare_dict.items():
    if value is None:
        print("failed to find a definition for {}".format(key))
user_45410450
Related