FASTA sequence extraction using list of sequence name that are not identical with the headers?

Question

[seta]

Hi all,

I 'm trying to extract some sequences from a fasta file that their header are as following:

>c0_g1_i1 len=526 path=[27:0-525]
TCTAGCTTGAAACCTGACATTAAAATGAAAAGGAGATCACTTGTAGTCCAATAAACCAAA
CTAGTTTTTAGAGAAATTTATTACTGTTTTATCTTCCCAACACTTGCAAATTTAGTAAGG
TATACATAAAAGTATGACATAATATTTAGATAAATAGATCTCTATTTTCATATGTTCAAG
AAAAAATCATCCTCATTTCTCAAAATTCTGAATAGAAAAGCTAAAATTCATGTGCTCTGA
TCATGTTATCCCTATGTCTAGCATTATTTGGTTCAGCCTGCGGTTAGTACAACTGAAGGG
TCACTGCCACGAAACTCCAACCATCGCAAATTCAAACATTTTACAAAAGAGCATGCGCTA
AATGACAATAATAGACTCCAGGTGATGTGTTCCATTGTTTCTAACACTTTGGCTATGGCA
AATTACTTGATAGCAAAATAATACAAGATGAGAAACCACACTTGAAATTACTTAAACCAC
TCTGCTCCAATTAGTCAGTTTGAACTACTAAACTAATAATTTCCTG
>c1_g1_i1 len=472 path=[27:0-471]
GAACTGGATTTAAATCAAGGATGTGGTTATGGATAGACTGCATGCGCATCTCACACATGA
AAACTAAGATCCCCTAGCACTCCCACCGTGCATCCCGCCATTCACAACAATCGTGATGTC
GCCCGGAACGTTTTCCATATGTCCAGGAGGTATCTTTGTCGCAATAACCTGGTTCTCCCG
CACAAACGATACCCGCAAGCAGTCCTGTGCACAGTCGAATACGGTGATGGATCCTGCTTT
GAACTCGTGGAAGGCCTTCTTGTCTTTCTCAGTGAAAGTCGCCTCGGTCCCCGTCACCCA
TAGGGAGTCAAGGTCCGATTTCGTAGTCATGTTTATTACCAAGGCTTGAGCTGCTTGCGT
TTCTAACTCTGCGGAATATGCTGCTGTATCTAGTACTTTGATGCCCAAAGGCAGCGGAAC
GGCCACGTCACAAATAGTTTAAGTTGGGCAGCAAATTAGGCGTGTGTCATGC

But, I have a text file containing sequence name that doesn't exactly match with fasta sequences header, it's like

>c0_g1_i1
>c1_g1_i1

Could you please help me out to extract my sequences of interest? Thanks in advance.

[5heikki]

grep -w -f file.txt file.fasta

should return the complete headers..

[seta]

If there is one space among name in the text file, the command return us the sequence with complete header. Now, I don't know how to create space among name in the text file?. It can be best command if it takes no long time.

Answer 1

[Anima Mundi]

In Python (for a FASTA file named foo.fa and a header list named header_list.txt):

def search(header):
    found = False
    for line in open('foo.fa'):
        if header.replace('\n', '') in line:
            print line,
            found = True
        elif '>' in line and found:
            break
        elif '>' not in line and found:
            print line,
            

for j in open('header_list.txt'):
        if len(j) > 2:
            search(j)

[seta]

Many thanks for all replies. I try them

[seta]

The command 's Anima Mundi sounds great. But, it takes long time to finish for a fasta file with about 170 MB. Is it usual?

[Anima Mundi]

The previous script was not designed to be particularly fast. How much did it take to run it on your machine? Did you manage? This one anyhow should be faster (though it could be made even faster!), hope it helps:

headers = []
for j in open('header_list.txt'):
    if len(j) > 2:
        headers.append(j.replace('\n', ''))

found = False
for line in open('foo.fa'):
    if '>' in line and not found:
        for item in headers:
            if item in line:
                print line,
                found = True
                break
    elif '>' in line and found:
        for item in headers:
            if item in line:
                print line,
                found = True
                break
            else:
                found = False
    elif '>' not in line and found:
        print line,

[seta]

Actually, I tested it on just 4 sequences to see if works or not. it run on my original file for 2-3 hours and it has not still finished. it seems unusual. I try new script, thanks

[seta]

Unfortunately, this script did not work at all. I executed it directly on my original file and after about 40-50 min, an empty output file was appeared, without any error!

[Anima Mundi]

It is working for me for the sequences and headers that you posted (now I added a few fake ones and it still works). What happens if you try to run it on a few sequences? What do you get with grep -c '^>' foo.fa?

[seta]

It works fine for two example sequences that I posted, but it did not work on my original file composed of about 130000 sequences. output of the command was 0 for the original file, strange!

[Anima Mundi]

Ok, I see. But next time please provide more information at the beginning (on the input files but also about your machine), as the size of your files is a fundamental part of the problem here ;).

Maybe this solution will work for you:

a) add a > after the last line of your foo.fa

b) grep '^>' header_list.txt > header_list_clean.txt

c) launch this Python script and output on a file called out.fa:

output = {}
curr_head = ''
curr_seq = ''
for line in open('foo.fa'):
    if '>' in line:
        output[curr_head] = curr_seq
        curr_head = ''
        curr_seq = ''
        curr_head += line.replace('\n', '')
    elif '\n' != line:
        curr_seq += line.replace('\n', '')

for keys,values in output.items():
    print(keys) + '\t' + (values)

d) as 5heikki suggests, grep -w -f header_list_clean.txt out.fa > out_desired.fa

e) refine your out_desired.fa with this Python script:

for line in open ('out_desired.fa'):
    print line.replace('\t', '\n'),

[seta]

Thanks friend. I finally clean fasta header so only first column of header reminded using awk '{print $1}' file.fa > output.fa, then extract desired sequence using the following python script:

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""Extract sequences from a fasta file if their name is in a 'wanted' file.

Wanted file contains one sequence name per line.

Usage:
    %program <input_file> <wanted_file> <output_file> """

import sys
import re

try:
    from Bio import SeqIO
except:
    print "This program requires the Biopython library"
    sys.exit(0)

try:
    fasta_file = sys.argv[1]  # Input fasta file
    wanted_file = sys.argv[2] # Input wanted file, one gene name per line
    result_file = sys.argv[3] # Output fasta file
except:
    print __doc__
    sys.exit(0)

wanted = set()
with open(wanted_file) as f:
    for line in f:
        line = line.strip()
        if line != "":
            wanted.add(line)

fasta_sequences = SeqIO.parse(open(fasta_file),'fasta')

with open(result_file, "w") as f:
    for seq in fasta_sequences:
        name = seq.id
        if name in wanted and len(seq.seq.tostring()) > 0:
            wanted.remove(name) # Output only the first appearance for a name
            SeqIO.write([seq], f, "fasta")

It was so easy and fast, just took a few seconds.

[Anima Mundi]

All's well that ends well, bye :)

[glich]

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""Extract sequences from a fasta file if their name is in a 'wanted' file.

Wanted file contains one sequence name per line.

Usage:
    %program <input_file> <wanted_file> <output_file>"""

import sys
import re

try:
    from Bio import SeqIO
except:
    print "This program requires the Biopython library"
    sys.exit(0)

try:
    fasta_file = sys.argv[1]  # Input fasta file
    wanted_file = sys.argv[2] # Input wanted file, one gene name per line
    result_file = sys.argv[3] # Output fasta file
except:
    print __doc__
    sys.exit(0)

wanted = set()
with open(wanted_file) as f:
    for line in f:
        line = line.strip()
        if line != "":
            wanted.add(line)

fasta_sequences = SeqIO.parse(open(fasta_file),'fasta')

with open(result_file, "w") as f:
    for seq in fasta_sequences:
        name = seq.id
        if name in wanted and len(str(name)) > 0:
            wanted.remove(name) # Output only the first appearance for a name
            SeqIO.write([seq], f, "fasta")

I changed this line:

if name in wanted and len(seq.seq.tostring()) > 0:

to:

if name in wanted and len(str(name)) > 0:

because

/usr/local/lib/python2.7/dist-packages/Bio/Seq.py:341: BiopythonDeprecationWarning: This method is obsolete; please use str(my_seq) instead of my_seq.tostring().

Answer 2

[venu]

grep '^>' original.faa > headers.txt
grep -Fwf your_headers.txt headers.txt > final_headers.txt
./faSomeRecords original.faa final_headers.txt output.faa

Download faSomeRecords here. Change mode and use it.

NOTE: Remove '>' from headers file while using faSomeRecords

Answer 3

[h.mon]

Several answers from this post should work with little or no modification for you.

Answer 4

[nterhoeven]

One of my colleagues wrote this tool, it does what you need. Just run

./SeqFilter --help

to get an overview of the command line options

Answer 5

[Tanvir Ahamed]

To extract sequence form a big list of file you can use faSomeRecords. Detailed instruction here.