extracting contigs

Question

[hollyannj7]

Hiya,

I have ordered my contigs by sequence length and wish to extract all of the contigs up to 310 (L50) or >000310F into a separate file.

Can this be done with grep?

Thanks!

[Carambakaracho]

I'm not sure I understand whether you want the shorter or the longer than L50 contigs.

Anyhow, provided your contigs are ordered by contig length as you wrote, you can simply get the last line of the contig at L50 and head/tail the file, for example

head -n $LINE_L50 <contig_all.fna >contig_above_L50.fna

[iraun]

Hi! Give a try to reformat.sh using maxlenght option (see BBMap reformat from BBMap toolkit).

Answer 1

[Mark]

The amazing seqkit can do this very easily:

seqkit seq seqs.fasta -m 310 > seqs.310.fasta

You can install seqkit using conda:

conda install -c bioconda seqkit

The docs are very good: https://bioinf.shenwei.me/seqkit/usage/#seq