looking for a simple script to pull a subset of contigs from an assembly

[lac302]

i'm sure this is a simple enough task but i'm just an end user no scripting experience at all.

looking for a script to pull contigs listed in a .txt file from assembly.fa and output the results to a new .fa file

any help would be much appreciated. thanks.

[Monika_bioinf]

Use seqtk subseq, https://github.com/lh3/seqtk

[GenoMax]

Just in case subseq does not work with fasta files here are other alternatives:

http://www.biostars.org/p/2822/

A perl one liner:

http://edwards.sdsu.edu/labsite/index.php/robert/381-perl-one-liner-to-extract-sequences-by-their-identifer-from-a-fasta-file

[JackieBadger]

Out of curiosity, GenoMax, how does second perl function handle the searching of the sample name file?
I wrote something similar under linux bash and it went bonkers with the RAM. I think the issue was that when ever it found an ID in the fasta file, it would then not remove this ID from the file containing the query IDs, and then start the search again from the start. Either way, my scripting abilities produced something that seemed unfeasible to execute on a large query ID file and large fasta file.

Cheers,

J

[lac302]

Thanks all!

[JohnN]

I wrote this a few years back.

https://code.google.com/p/nash-bioinformatics-codelets/downloads/detail?name=subset_fasta.pl&can=2&q=

[Richard Finney]

awk 'BEGIN{while((getline x<ARGV[1])>0){a[i++]=x;}while((getline y<ARGV[2])>0){if(substr(y,0,1)==">"){m=0;for(j=0;j<i;j++){if(y==a[j])m=1;}}if(m==1)print y;}}' $1 $2


$1 is match file
$2 is fasta file

[GenoMax]

Out of curiosity, GenoMax, how does second perl function handle the searching of the sample name file?
I wrote something similar under linux bash and it went bonkers with the RAM. I think the issue was that when ever it found an ID in the fasta file, it would then not remove this ID from the file containing the query IDs, and then start the search again from the start. Either way, my scripting abilities produced something that seemed unfeasible to execute on a large query ID file and large fasta file.

Cheers,

J

@JackieBadger: Second perl function is using -n and -e switches. -n wraps a while loop around the program while -p feeds the program value of $_ each time.

A nice example that illustrates this (equivalent to unix 'cat' command)

$ perl -ne 'print $_' filename

or

$ perl -ne 'print' filename

[Birdman]

This little BioPython script will nicely do the job:

from Bio import SeqIO
import sys

#Usage: filter_fasta_per_ids.py input.fasta filter_ids.txt output.fasta

input_file =sys.argv[1]
id_file =sys.argv[2]
output_file =sys.argv[3]
wanted = set(line.rstrip("\n").split(None,1)[0] for line in open(id_file))
print("Found %i unique identifiers in %s" % (len(wanted), id_file))
records = (r for r in SeqIO.parse(input_file, "fasta") if r.id in wanted)
count = SeqIO.write(records, output_file, "fasta")
print("Saved %i records from %s to %s" % (count, input_file, output_file))
if count < len(wanted):
    print("Warning %i IDs not found in %s" % (len(wanted)-count, input_file))