Extracting specific FASTA sequences from a complex FASTA file

Question

[boileau]

I have been trying to troubleshoot this for the past couple days but I can't seem to find the issue, any help you could provide is greatly appreciated. I am a noob but will be writing similar scripts to this in the future- so any added resources are welcomed.

I have a large file with a bunch of FASTA sequences and I am trying to extract specific ones based on a list of IDs. While I would like the output to be a collection of whatever ID == FASTA sequences it won't give me an output with results. I know that if I change t to a definite string (ie. 'NM_xxxxx') then it will give a desired output so I suspect the issue is in reading from the csv file.

import re
import sys
import csv
from Bio import SeqIO

output = open("matches.fasta", "w")
x = SeqIO.parse(open(sys.argv[1], 'rU'),'fasta')
y = csv.reader(open(sys.argv[2], 'rU'))

for row in y:
    print str(row[1])
    t = str(row[1])
    for fasta in x:
            name, sequence = fasta.id, fasta.seq.tostring()
            searchObj = re.search( t, fasta.id, re.M|re.I)
            if searchObj:
                SeqIO.write(fasta, output, "fasta")

output.close()

[Ram]

I hope you read through the dozens of ID matching posts here.

[boileau]

I did, but the fact my experience is almost exclusively python and very limited at that means I am unable to take advantage of related posts.

[Ram]

From your description, there could be one problem with the data that is a bit undetectable unless you look for it. For a few sequences, after reading them with Bio.SeqIO, try printing their ID, name and description attributes. Check if your expectations of what should be printed are in fact what are printed.

If the FASTA header has white spaces, the SeqIO parser can be a bit crazy. I'd usually say that the part before spaces goes into id and the part after goes into description, but this assumption was recently proven wrong.

Check out if the data conforms to expectations, and then dig into your code for changes.

[boileau]

The main issue is that the FASTA header has a wierd header (ie. >NM_xxxxxxx_2000_from_chromosome_spongebob) which is why I was trying to use RE to match with my candidate list containing just the RefSeq Ids (NM_xxxxxx). Else, Bio.SeqIO is behaving like it should given the control tests I've run.

[Ram]

Ah, I see. Have you tried printing str(row[1]), as plain text as well as an ASCII sequence? Maybe a new line or a black space at the end (or a weird unicode character) is messing up your RE search.

[boileau]

You were right, looks like the export .txt left some weird ascii characters. Thanks for the help Ram!

[barvazduck]

In linux

#exract fasta sequences with a reference list ALL SEQUENCE LENGTH

#LIST is your list as text and IN is your fasta
>awk 'BEGIN{while((getline<"LIST")>0)l[">"$1]=1}/^>/{f=l[$1]?1:0}f' IN.fa > OUT.fa

[Ram]

Couple of notes:

a. LIST is a file, not text. substituting that can be a bit non-trivial

b. Heng Li's bioawk is awk for biological formats. It will obviate the need for multi-line checks and header line awareness in your awk script.

Answer 1

[dago]

Try pyfaidx, it is really easy to use.

pip install pyfaidx

then to extract your seq

xargs faidx -d "<here the separator used in the IDs of your fasta>" <Your fastafile.fasta> < <Here the file containing the IDs>

Example

xargs faidx -d "|" My_seq.fasta < My_Ids.txt

PS: this is a clear example of the great "teaching power" of Biostars. extract sequences from multi fasta using partial ID

Answer 2

[GouthamAtla]

What you are missing is seek() in biopython. Assuming the seq_ids in a file called "names"

How to use seek() with Biopython SeqIO object

from Bio import SeqIO

fasta=SeqIO.parse("fasta.fa","fasta")
seq_dict={}

for record in fasta:
        seq_dict[record.id]=record.seq
# If you feel its memory intensive, try next solution

for line in open("names","r"):
        line=line.strip()
        if line in seq_dict.keys():
                print ">"+line+"\n"+seq_dict[line]

Otherwise write a small function like:

from Bio import SeqIO

def getSeqRec(seq_id):
        fasta=SeqIO.parse("fasta.fa","fasta")
        for record in fasta:
                if seq_id==record.id:
                        return record

for line in open("names","r"):
        line=line.strip()
        seq_rec=getSeqRec(line)
        print ">"+seq_rec.id
        print seq_rec.seq

Answer 3

[mark.ziemann]

Does it need to be Python? samtools faidx does it.

samtools faidx Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa chr22

Answer 4

[barvazduck]

In linux

#exract fasta sequences with a reference list ALL SEQUENCE LENGTH
#LIST is your list as text and IN is your fasta

awk 'BEGIN{while((getline<"LIST")>0)l[">"$1]=1}/^>/{f=l[$1]?1:0}f' IN.fa > OUT.fa

Answer 5

[shenwei356]

Here is a script I wrote. Before using it, Perl and perl module BioUtil should be install.

fasta_extract_by_pattern.pl could extract FASTA sequences by header or sequence, exactly matching or regular expression matching are both supported. The query pattern could read from files. And negation of the result is also easy to get. What's the most important, it could read from STDIN.

Combining fasta2tab and tab2fasta with cvs_grep could also have the same function.

Examples

1. sequences WITH "bacteria" in header

fasta_extract_by_pattern.pl -r -p Bacteria *.fa > result.fa

2. sequences WITHOUT "bacteria" in header

fasta_extract_by_pattern.pl -r -n -p Bacteria seq1.fa seq2.fa > result.fa

3. sequences with TTSAA (AgsI digest site) in SEQUENCE. Base S stands for C or G.

fasta_extract_by_pattern.pl -r -s -p 'TT[C|G]AA' seq.fa > result.fa

4. sequences (read from STDIN ) with header that matches any patterns in list file

zcat seq.fa.gz | fasta_extract_by_pattern.pl -pf name_list.txt > result.fa

See more: Manipulation on FASTA/Q format file

Answer 6

[Brian Bushnell]

To extract sequences with specific names from a large file efficiently, I suggest you use my filterbyname tool:

A: Extract sequence with header from a fasta file with specific ID given in another