This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How To Remove Certain Sequences From A Fasta File

Hi everyone

Is there a fast way to do this filter?

I have a huge Fasta file (sequences are short reads coming from an Illumina instrument). I have also a list of nucleotide sequences (not Fasta, just the sequences) and I want to remove from the big Fasta file all entries identical to those in the list.

My idea was simply to go down through the Fasta file and then, for every read, check all the sequences of the list. If the read matches one of the sequences then do nothing, otherwise print the read into a new file. I made this with perl but it takes ages!

The list is made up of nucleotide sequences, not IDs. It's something like this:

AACGACTACTTATCGATC

TCGGCGATATACGTAC

CCAGTTTCGGGGCTAT ....

Thanks!

fasta parsing short
• 35,064 views • link
updated by Eric Normandeau Eric Normandeau • written by Emmecola Emmecola

can you make a better example... you have a file with a lot of sequences, and a list of sequence ids to remove from it. Is it correct?

• 0 views • link
written by Giovanni M Dall'Olio Giovanni M Dall'Olio

Yes, except that the list is not made up of IDs. They are nucleotide sequences.

• 0 views • link
written by Emmecola Emmecola

Do I understand correctly that you have a list of nucleotide sequences, and you want to go through a FASTA file and remove all entries that contain an exact, full-length hit to one of the nucleotide sequences in your list? Are the nucleotide sequences in your list by any chance all of the same length?

• 0 views • link
written by Lars Juhl Jensen Lars Juhl Jensen

I'm fairly sure this has been asked/answered before - check the "Related" box or search the archives.

• 0 views • link
written by Neilfws Neilfws

3 answers

  • linearize your fastq file into 3 columns: sequence, ID, qualities (using 'awk', or 'perl')
  • sort it using the sequence as the key ( unix 'sort')
  • sort your 2nd file of sequences ( unix 'sort')
  • use unix 'join' to filter out the first file and re-transform it to FASTQ using awk or perl
• 2 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Pierre's solution is the way to go for sure, assuming you don't need to account for mismatches, INDELs, etc. Use the -v option in join to exclude the FASTQ sequences that match your filter set.

The beauty of this recipe is that is applies to many other scenarios and is not constrained by memory. I probably do something close to this recipe (in other situations) 100 times a year.

• 0 views • link
updated by Ram Ram • written by Aaronquinlan Aaronquinlan

Hi,

EDIT: I had misunderstood that you actually had the nucleotide sequences in the exclusion file. I modified the code accordingly.

You can use the following script. Put it in a file named fasta_remove.py, make the file executable:

chmod +x fasta_remove.py

and use it by typing:

./fasta_remove.py <input_fasta> <remove_file> <output_fasta>

The remove_file file contains one sequence name per line that needs to be removed from the input_fasta file. The input fasta file is not limited in size. You will need to have Biopython installed and this works under Linux, obviously.

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

import sys
from Bio import SeqIO

fasta_file = sys.argv[1]  # Input fasta file
remove_file = sys.argv[2] # Input wanted file, one gene name per line
result_file = sys.argv[3] # Output fasta file

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:
        nuc = seq.seq.tostring()
        if nuc not in remove and len(nuc) > 0:
            SeqIO.write([seq], f, "fasta")

Cheers!

• 2 views • link
updated by Ram Ram • written by Eric Normandeau Eric Normandeau

Hi

I am working with old sanger reads, and I have the fasta and qual files separated. I will like to remove some sequences from the qual file using a list of names of those sequences but I am having problems with these. ¿Do yo have a solution for this?

Thankyou very much

• 0 views • link
written by shinken123 shinken123

For me, I had to change the code to this:

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

import sys
from Bio import SeqIO

fasta_file = sys.argv[1]  # Input fasta file
remove_file = sys.argv[2] # Input wanted file, one gene name per line
result_file = sys.argv[3] # Output fasta file

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 = seq.id
        **nuc = str(seq.seq)
        **if nam not in remove and len(nuc) > 0:
            SeqIO.write([seq], f, "fasta")

(I changed the lines beginning with '**' )

Versions: Python 2.7.11 and Biopython 1.67

• 1 views • link
written by tlaguna • tlaguna

Hi Eric. Thanks for the script. Do you have a suggestion so as how the 'wanted' file should look exactly? I tried with text file like this: 

>TRINITY_DN10009_c0_g1_i1   
>TRINITY_DN103518_c0_g1_i1
>TRINITY_DN104798_c0_g1_i1

or just to be sure also like this

TRINITY_DN10009_c0_g1_i1    
TRINITY_DN103518_c0_g1_i1
TRINITY_DN104798_c0_g1_i1

but the script didn't work. The output file was the same as input.

Any suggestions for the troubleshooting? or other script? MilleThanks.

• 0 views • link
updated by Ram Ram • written by marie.pazoutova • marie.pazoutova

Hi I tried to use the script but it needs to be updated. I am not an expert so just changed my_seq.tostring() to str() but it is not working. Also does the fasta file needs to have a single line per sample?

Message:

'BiopythonDeprecationWarning: This method is obsolete; please use str(my_seq) instead of my_seq.tostring().
  BiopythonDeprecationWarning)'

Thanks

• 0 views • link
updated by Ram Ram • written by catagui • catagui

An alternative to Pierres suggestion: use

grep -v SeqToBeDeleted inputfile > outputfile

and iterate on all your Sequences to be removed.

• 0 views • link
updated by Ram Ram • written by Grünschnabel • Grünschnabel

Unless the number of sequences that you want to remove is very small, this will be prohibitively slow. In either case, if you want to use grep for it, here are two simple optimations: 1) use the -f option to give grep a file with all sequences that you want to search with instead of running grep multiple times, and 2) use fgrep instead of grep since you only need to do exact string matching and not full regular expression matching.

• 0 views • link
written by Lars Juhl Jensen Lars Juhl Jensen

Unless the number of sequences that you want to remove is very small, this will be prohibitively slow. In any case, if you want to use grep for it, here are two simple optimizations:

  1. use the -f option to give grep a file with all sequences that you want to search with instead of running grep multiple times, and
  2. use fgrep instead of grep since you only need to do exact string matching and not full regular expression matching.
• 0 views • link
updated by Ram Ram • written by Lars Juhl Jensen Lars Juhl Jensen

Plus it will only remove the sequence part of the FASTQ entry and will result in a mangled file consisting of proper FASTQ entries that should remain and incomplete (i.e. head, separator, qual, yet no sequence) entries for the matched sequences.

• 0 views • link
written by Aaronquinlan Aaronquinlan

Log in to answer this question.

More posts like this