Extraction of first sequences from a big fasta file

Question

[vahapel]

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!

[Devon Ryan]

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.

[Tanvir Ahamed]

You can have a look at this. May be helpful for you

Extract sequence with header from a fasta file with specific ID given in another file

[GenoMax]

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

Answer 1

[dariober]

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.

[vahapel]

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

Answer 2

[Brian Bushnell]

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.

[vahapel]

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

[Brian Bushnell]

You're welcome!

Answer 3

[5heikki]

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

[vahapel]

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

[st.ph.n]

likewise, use tail for the last N sequences.

Answer 4

[AK]

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

Answer 5

[Simply Bioinformatics]

REQUESTED_LINES=10
awk "/^>/ {n++} n>$REQUESTED_LINES {exit} {print}" input.fasta > output.fasta

Answer 6

[psschlogl]

hdl = gzip.open(file, 'rt')
records = SeqIO.parse(hdl, 'fastq')
first_read = next(records)
printfirst_read.id)

[lieven.sterck]

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.