This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Extraction of first sequences from a big fasta file

Dear All,

I would like to ask a question regarding extraction of 100000 sequences in a big fasta file. In the forum, there is a bunch of script handling the sequence extraction based on ID number, but I could not find a script for such a purpose. Basically, is there any script or bash command for extraction first and/or last 100000 sequences in a fasta file?

Many thanks in advance for all your help!

sequence

If not, you could write one easily enough with biopython or bioperl...

Well, the first X records is easier than the last X records, but still.

That is not an answer for the question that was originally asked.

6 answers

This strategy is based on standard nix tools. Get the *first two sequences:

awk -v RS='>' 'NR>1 { gsub("\n", ";", $0); sub(";$", "", $0); print ">"$0 }' seq.fa \
    | head -n 2 \
    | tr ',' '\n'

It assumes the semicolon ; doesn't occur in the sequence names.

Replace head with tail to get the last sequences.

Hi, dariober thank you so much for this script. It works well

For the first 100k sequences:

reformat.sh in=data.fasta out=100k.fasta reads=100000

You can also get a random 100k sequences with Reformat, but not the last 100k.

Hi Brian, thank you for your help and introducing "bbmap" tools for me, it makes my project more easy

Assuming no linebreaks in sequences, i.e every record is exactly two lines. I believe for fastq files the multiplier would be 4:

head -n 200000 input > output

thank you, 5heikki for your help, it is very simple and useful command.

likewise, use tail for the last N sequences.

Alternatives: seqkit head and seqkit range:

(1) Leading 100000 records:

seqkit head -n 100000 input.fa
seqkit range -r 1:100000 input.fa

(2) Last 100000 records:

seqkit range -r -100000:-1 input.fa

(3) Other ranges:

seqkit range -r 100001:200000 input.fa
REQUESTED_LINES=10
awk "/^>/ {n++} n>$REQUESTED_LINES {exit} {print}" input.fasta > output.fasta
hdl = gzip.open(file, 'rt')
records = SeqIO.parse(hdl, 'fastq')
first_read = next(records)
printfirst_read.id)

Hi psschlogl

Though your answer shows some potential it does not really answer the question posted. If you can provide a more complete/correct answer we could keep it in the answers thread. If not, we consider removing it.

Log in to answer this question.