Extracting lines from text files - script with a couple of 'side effects'

M

mstagliamonte

Dear All,

Here I am, with another newbie question. I am trying to extract some lines from a fasta (text) file which match the headers in another file. i.e:
Fasta file:
header1|info1:info2_info3 general text
header2|info1:info2_info3
general text

headers file:
header1|info1:info2_info3
header2|info1:info2_info3

I want to create a third file, similar to the first one, but only containing headers and text of what is listed in the second file. Also, I want to print out how many headers were actually found from the second file to match the first.

I have done a script which seems to work, but with a couple of 'side effects'
Here is my script:
-------------------------------------------------------------------
import re
class Extractor():

def __init__(self,headers_file, fasta_file,output_file):
with open(headers_file,'r') as inp0:
counter0=0
container=''
inp0_bis=inp0.read().split('\n')
for x in inp0_bis:
container+=x.replace(':','_').replace('|','_')
with open(fasta_file,'r') as inp1:
inp1_bis=inp1.read().split('>')
for i in inp1_bis:
i_bis= i.split('\n')
match = re.search(i_bis[0].replace(':','_').replace('|','_'),container)
if match:
counter0+=1
with open(output_file,'at') as out0:
out0.write('>'+i)
print '{} sequences were found'.format(counter0)

-------------------------------------------------------------------
Side effects:
1) The very first header is written as >>header1 rather than >header1
2) the number of sequences found is 1 more than the ones actually found!

Have you got any thoughts about causes/solutions?

Thanks for your time!
P.S.: I think I have removed the double posting... not sure...
Max
 
D

Dave Angel

Dear All,

Here I am, with another newbie question. I am trying to extract some lines from a fasta (text) file which match the headers in another file. i.e:
Fasta file:
header1|info1:info2_info3 general text
header2|info1:info2_info3
general text

headers file:
header1|info1:info2_info3
header2|info1:info2_info3

I want to create a third file, similar to the first one, but only containing headers and text of what is listed in the second file. Also, I want to print out how many headers were actually found from the second file to match the first.

I have done a script which seems to work, but with a couple of 'side effects'
Here is my script:
-------------------------------------------------------------------
import re
class Extractor():

def __init__(self,headers_file, fasta_file,output_file):
with open(headers_file,'r') as inp0:
counter0=0
container=''
inp0_bis=inp0.read().split('\n')
for x in inp0_bis:
container+=x.replace(':','_').replace('|','_')
with open(fasta_file,'r') as inp1:
inp1_bis=inp1.read().split('>')
for i in inp1_bis:
i_bis= i.split('\n')
match = re.search(i_bis[0].replace(':','_').replace('|','_'),container)
if match:
counter0+=1
with open(output_file,'at') as out0:
out0.write('>'+i)
print '{} sequences were found'.format(counter0)

-------------------------------------------------------------------
Side effects:
1) The very first header is written as >>header1 rather than >header1
2) the number of sequences found is 1 more than the ones actually found!

Have you got any thoughts about causes/solutions?

The cause is the same. The first line in the file starts with a "<",
and you're splitting on the same. So the first item of inp1_bis is the
empty string. That string is certainly contained within container, so
it matches, and produces a result of ">"

You can "fix" this by adding a line after the "for i in inp1_bis" line
if not i: continue

But it also seems to me you're doing the search backwards. if the Fasta
file has a line like: >der

it would be considered a match! Seems to me you'd want to only match
lines which contain an entire header.
 
M

MRAB

Dear All,

Here I am, with another newbie question. I am trying to extract some lines from a fasta (text) file which match the headers in another file. i.e:
Fasta file:
header1|info1:info2_info3 general text
header2|info1:info2_info3
general text

headers file:
header1|info1:info2_info3
header2|info1:info2_info3

I want to create a third file, similar to the first one, but only containing headers and text of what is listed in the second file. Also, I want to print out how many headers were actually found from the second file to match the first.

I have done a script which seems to work, but with a couple of 'side effects'
Here is my script:
-------------------------------------------------------------------
import re
class Extractor():

def __init__(self,headers_file, fasta_file,output_file):
with open(headers_file,'r') as inp0:
counter0=0
container=''
inp0_bis=inp0.read().split('\n')
for x in inp0_bis:
container+=x.replace(':','_').replace('|','_')
with open(fasta_file,'r') as inp1:
inp1_bis=inp1.read().split('>')
for i in inp1_bis:
i_bis= i.split('\n')
match = re.search(i_bis[0].replace(':','_').replace('|','_'),container)
if match:
counter0+=1
with open(output_file,'at') as out0:
out0.write('>'+i)
print '{} sequences were found'.format(counter0)

-------------------------------------------------------------------
Side effects:
1) The very first header is written as >>header1 rather than >header1
2) the number of sequences found is 1 more than the ones actually found!

Have you got any thoughts about causes/solutions?

Thanks for your time!
Here's my version:

class Extractor():
def __init__(self, headers_file, fasta_file, output_file):
with open(headers_file) as inp:
headers = set('>' + line for line in inp)

counter = 0
accept = False

with open(fasta_file) as inp, open(output_file, 'w') as out:
for line in inp:
if line.startswith('>'):
accept = line in headers
if accept:
counter += 1

if accept:
out.write(line)

print '{} sequences were found'.format(counter)
 
M

mstagliamonte

Hi,

Thanks for the answers! And Dave, thanks for explaining the cause of the problem I will keep that in mind for the future. You're right, I am doing thesearch backward, it just seemed easier for me to do it in that way. Looks like I need to keep practising...

Both your suggestions work, I will try and learn from them.

Have a nice day
Max
 

Ask a Question

Want to reply to this thread or ask your own question?

You'll need to choose a username for the site, which only take a couple of moments. After that, you can post your question and our members will help you out.

Ask a Question

Members online

No members online now.

Forum statistics

Threads
473,813
Messages
2,569,696
Members
45,483
Latest member
TedDvb6626

Latest Threads

Top