EXTRACT MAQ sequences?

[johnsequence]

how to simply extract the sequences of a gene list (~1000) in FASTA format from a sequence database (~400MB) in FASTA format generated by MAQ?

[maubp]

Are you starting with a ~400MB FASTA file, containing ~1000 sequences, and you just want the list of sequence identifiers ("gene names")?

Try something like this at the Unix command line:

grep "^>" my_database.fasta

That string "^>" is a regular expression meaning look for any lines starting ("^") with the greater than symbol.

[johnsequence]

Thanks for reply. I actually need the IDs (headers) and sequences in FASTA format.

[maubp]

Thanks for reply. I actually need the IDs (headers) and sequences in FASTA format.

So you have a large FASTA file, and a list of entries you want to extract from it, to give a smaller subset as a new FASTA file? This is a very general problem, and not specific to MAQ at all.

How is your list of identifiers stored? e.g. a text file with one id per line?

I would suggest you write a simple script, e.g. using Perl (perhaps with BioPerl) or Python (perhaps with Biopython), or your preferred script language.

http://bioperl.org/wiki/HOWTO:SeqIO

http://biopython.org/wiki/SeqIO

Or, if you are happier just working at the command line, you can probably do this with EMBOSS seqret.

http://emboss.sourceforge.net/apps/release/6.1/emboss/apps/seqret.html

[kmcarr]

You can use a couple of the utilities in the BLAST package from NCBI. Take your large FASTA file and create a BLAST database from it using formatdb. Then retrieve just the sequences you want from the BLASTdb using the fastacmd tool.

%> formatdb -i <your.FASTA.file> -p F -n <your.blast.db>

%> fastacmd -d <your.blast.db> -i <your.ID.file> > <output.file>

The first command takes your FASTA file and creates your.blast.db. The -p F tells formatdb that this is a nucleotide database. The second command extracts the FASTA formatted reads from the BLAST db based on the list of of IDs in your.ID.file.