This is a test version of Biostars. For the public version, visit https://www.biostars.org.
There are any way to use a script for feed emboss with fasta sequence?

Hi there
I have a compressed fasta file (really big) and I was thinking if it is possible to use the SimpleFastaParser(biopython), gzip.open to feed the sequeces into emboss compseq through cml?

Thank you

sequence
• 4,300 views • link
written by schlogl 18

What is cml? Have you tried simple shell scripting - why do you need to use biopython?

• 0 views • link
written by Ram 4.5K

to read the compressed file together with gzip and then pass/stream to emboss each sequence at time.

And because a don't know much about bash scripting.

• 0 views • link
written by schlogl 18

You can do everything you need to do inside (bio)python:

http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc54

You haven't said exactly how large your file is, and python is not famed for its speed, so even with the best code in the world I still can't guarantee that this will be fast enough for you.

It seems EMBOSS will not accept piped output however, so without much hacking, genomax's advice of finding a tool that natively supports compressed amino acid data would be preferable. Through python you could make many iterative subprocess calls though once biopython has parsed the entry out for you (not the same as truly piping the data though).

• 0 views • link
written by Joe 2.2K

If you are thinking about feeding as in streaming in then most likely not.

• 0 views • link
written by GenoMax 16K

You probably can use something like this seqtk subseq in.fq name.lst > out.fq, but pipe the results into the other program rather than to a file.

• 0 views • link
written by onestop_data 33

I will checkit out if seqtk works with compressed files. Thank you

• 0 views • link
written by schlogl 18

You can use other options like jellyfish, kmercountexact.sh from BBMap if you are only looking get unique k-mer counts from your file. These programs will work with compressed files directly (ones from BBMap will)..

• 0 views • link
written by GenoMax 16K

Hi genomax

jeelyfish works with proteins too? I don't think so, sorry I forgot to say that the file is fasta.aa...

Thanks for your time.

• 0 views • link
written by schlogl 18

You could try this program. It was important to note that you have protein sequences.

• 0 views • link
written by GenoMax 16K

I don't see any reason that your described approach wouldn't work. Have you already tried it?

• 0 views • link
written by Joe 2.2K

How big is really big?

• 0 views • link
written by Joe 2.2K

Hi Joe it is the same file we already discussed in another question... 67GB! I am counting 5mers using python, but I would like to use another tool as a confirmation.

I have some code to read the file as a generator, but I reading the emboss manual and about seqcomp.

I tried to pass the sequences like this:

Calculate the composition of unique words in sequences
Input sequence(s): FASTA::hiv
Error: Failed to open filename 'hiv'
Error: Unable to read sequence 'FASTA::hiv'
Input sequence(s): fasta::hiv.txt
Error: Failed to open filename 'hiv.txt'
Error: Unable to read sequence 'fasta::hiv.txt'
Died: compseq terminated: Bad value for '-sequence' and no more retries

As specified in here: http://www.csd.hku.hk/bruhk/emboss.html

The test file in this case has only one sequence, but i always got this error.

What I am doing wrong?

PS- You right Joe it seems EMBOSS it is not a good tool in my case.
I will try genomax suggestion.
Thanks anyway.

• 0 views • link
written by schlogl 18

I'm not sure I understand what you're showing me (what is :: supposed to be?)

In order to pass to compseq you'll need to reconstitute a file to pass to -sequence. If you really want to do this with python still, take a look at the tempfile module.

The steps at a high level would be:

  1. Read compressed file
  2. Use biopython to get 1 record at a time
  3. Create a temporary file as a fasta or similar via biopython+tempfile
  4. Use tempfile with subprocess to dispatch compseq jobs to the shell
  5. Optionally purge the tempfile once the output file is generated.

This will be reasonably memory efficient I should think (though I'm not familiar with how python's gzip module works), but will be slowwww (if you have any high-speed storage, such as that which /tmp/ is usually mounted to, I'd strongly advise using that )

I used a similar approach here if you want some reference material: https://github.com/jrjhealey/Oread/blob/master/oread/__main__.py#L208-L246

• 0 views • link
written by Joe 2.2K

Thats what I was thinking to do some similar like that...

I never worked with subprocess and tempfile modules, but I will check it out.

Thanks Joe.

• 0 views • link
written by schlogl 18

That's a scary code 8)

I don't know if I am skilled enough to write some like that, but I will try.

Paulo

• 0 views • link
written by schlogl 18

You don't need everything thats in the code I linked. If you just look at the documentation/examples for each of the modules, they'll show you what you need to as a minimum.

It will be tricky, as thats a lot of steps to glue together, but this is the sort of thing that python is actually particularly good at so this will be good practice.

• 0 views • link
written by Joe 2.2K

Thanks Joe.

I appreciate you effort to help.

Thank you very much.

• 0 views • link
written by schlogl 18

You can try parallel. It can send in chunks of lines or bytes, of input file. But binary file I am not sure.

• 0 views • link
written by cpad0112 2.1K

0 answers

No answers yet.

Log in to answer this question.

More posts like this