I have a FASTA file of bacteria which has multiple strains like below. I would like to go through the FASTA file and extract/separate the largest strain from the FASTA file. Does anyone help me how i can do that?
>NZ_PKKG01000001.1 Lactobacillus crispatus strain UMB1398 .21837_8_80.1, whole genome shotgun sequence
TCTAATTTCACGGCTGATAGTTGATTTATGACAGCCTACTTCTTGAGCGATAGCAGTGCAGGAGGTTATGCCAGAATCAA
GCAGAGCTTGAATTTGTCCACGTTGTTCGCTGTTCAATTGATGATAATGCTTGGAAATGCTAGAATTTGAGTTGGTCATG
AAGATCTTCCTTTCTTGATTTTTGGTCACTTCAAGTTTAGGTCTTCATGGCCTTTTTGTTTAACAATTAGTGTTGCACTT
>NZ_PKKG01000101.1 Lactobacillus crispatus strain UMB1398 .21837_8_80.101, whole genome shotgun sequence
AAGGGCAAGAGCATGAAAGACAAGTCAAAAGCTTATGTTCAACAGGCAACTGATGCCATTAACCACAAATATCGCCGAAT
CCTCCAATATCACACAGCAGAGGAACTCTTCAAGCAATATATCTCTTCATAACCTAACTGTTGCACTTAATTTGACAATT
CAGGCAACTTGTAAATTGAACAAAAAATAGGGAAGATGAAACTAATCATTTCCCTATTTGTTTTTAAGAACAGTATTCAA
...
•
written
by
Shelle •
2 answers
Get the length of each sequence using the following command -
cat file.fa | awk '$0 ~ ">" {print c; c=0;printf substr($0,2,100) "\t"; } $0 !~ ">" {c+=length($0);} END { print c; }' reference
Try this
longest=$(seqkit sort --quiet -lr your.fasta.file.fa | seqkit fx2tab -nl | head -1 | cut -f1); seqkit grep -n -r -p "$longest" your.fasta.file.fa > longest.fa
Explanation
A. Fetching the header of the longest sequence
longest=$(seqkit sort --quiet -lr your.fasta.file.fa | seqkit fx2tab -nl | head -1 | cut -f1);
seqkit sort --quiet -lr your.fasta.file.f : sort by length
seqkit fx2tab -nl: take the lengths of sorted fasta file
head -1 | cut -f1 : first line would have the header and length of the longest sequence. head -1 is for taking first line and cut -f1 is to take just the header. This is saved in the variable called longest.
B. Fetching the sequence
seqkit grep -n -r -p "$longest" your.fasta.file.fa > longest.fa
grep : fetch header stored in longest variable
-n: by full name
-r: using regex (regular expression)
-p: pattern
@cpad : did that, longest is more appropriate. I liked the shortened version ;)
Hello shelle,
Please use the formatting bar (especially the code option) to present your post better. I've done it for you this time.
This question has been asked many times. A solution using seqkit
seqkit fx2tab -nl your.fasta.file.fa | sort -nr -k2 | head -1
Thank you!
Log in to answer this question.
-
How to create a tax_id_map file for getting scientific names in BLASTN output
written by Maria •
-
How do I extract a fasta sequence using the specific keyword (**Psychro**) ?
written by sunnykevin97
-
need to remove some part of header fasta file
written by Princy •
-
How to change headers of fasta files with the help of another file (nucl_wgs.accession2taxid.gz) co…
written by Sebastian •
-
Parse multifasta files based on a values in two columns in a metadata file
written by jmwhitha •
-
Better method to find specific enzymes in whole genome sequences?
written by Mel161803 •
-
How to separate the largest strain when it has the genome of its chromosome
written by Shelle •
-
Change NCBI fasta file headers to makeblastdb format
written by chland •
-
accession to taxid for fasta files downloaded from refseq
written by badribio
-
how to locate a gene sequence among fastq files containing short reads
written by jerrybug109 •
Outputs the sequences sorted by length and largest one on the top. You can use
seqkit sort --quiet -2lr sequence.fa |seqkit head -n 1for the largest sequence (provided there are no two equal length largest sequences) . Download seqkit from here. If you are a programmer, you can also use index of your fasta file (.fai from samtools faidx) to getting largest sequence.