Hi, dariober thank you so much for this script. It works well
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!
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.
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
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)
Log in to answer this question.
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.
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
That is not an answer for the question that was originally asked.