How to extract fasta headers that contains specific text and it's sequence

Question

[Optimist]

Greetings!!!

I have barrnap output of 100 Pseudomonas aeruginosa genomes.

The output looks like this (sequences have been trimmed to avoid huge lines in biostars)

>16S_rRNA::Pseudomonas_aeruginosa_PAOC_Seq_1:6516148-6517679(-)
TGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACACATGCAAGTCGAGCGGATGAAGGGAGCTTGCTCCTGGATTCAGCGGCGGACGGGTGAG
>23S_rRNA::Pseudomonas_aeruginosa_PAOC_Seq_1:6512786-6515674(-)
TCAAGTGAAGAAGCGCATACGGTGGATGCCTTGGCAGTCAGAGGCGATGAAAGACGTGGTAGCCTGCGAAAAGCT
>5S_rRNA::Pseudomonas_aeruginosa_PAOC_Seq_1:6512529-6512639(-)
TGACGATCATAGAGCGTTGGAACCACCTGATCCCTTCCCGAACTCAGAAGTGA

I want to extract only 16s rRNA headers and sequences from all the outputs.

Result output should look like this

>16S_rRNA::Pseudomonas_aeruginosa_PAOC_Seq_1:6516148-6517679(-) 
TGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACACATGCAAGTCGAGCGGATGAAGGGAGCTTGCTCCTGGATTCAGCGGCGGACGGGTGAG

How can I get this output

Thank you all

[lieven.sterck]

What have you tried so far?

did you had a look at for instance seqtk or SeqKit ?

[Optimist]

filter_fasta.py has worked.

[cpad0112]

with awk and flattend fasta:

awk '/^>16S/{print;getline;print}' seq.fa

[Neel]

Hi optimist,

can you share your script? how did you used barrnap on 100 P.aeruginosa genome ?

Thank you so much!

Answer 1

[Optimist]

I have found a solution to this.

seqkit rmdup -n -i -j <threads> <infilename> > outfilename

Thank you all for the responses

Cheers Have a great time!!!

[cpad0112]

This would remove dups by sequence ID, not by sequence. Try using seqkit rmdup -s -i -j <threads> <infilename> > outfilename

Answer 2

[5heikki]

If there are no linebreaks in the sequences I would do it like this:

paste - - <file.fasta | awk 'BEGIN{FS="\t";OFS="\n"}{if($1~/16S/){print $1,$2}}'

If you really have duplicate headers, i.e. identical strings then you can (note if the sequences under identical headers differ you will lose that information here):

paste - - <file.fasta | awk 'BEGIN{OFS=FS="\t"}{if($1~/16S/){print $0}}' | sort -t $'\t'-uk1,1 | awk 'BEGIN{FS="\t";OFS="\n"}{print $1,$2}'

Answer 3

[Fatima]

If each sequence is one and only one line:

grep -A 1 "16S_rRNA:" filename

Should do the trick!

[Optimist]

This has worked too. Thank You

is there a way to remove duplicate headers along with the associated fasta sequence

Looks like most of the genomes have more than 1 16SrRNA sequences (as expected).

[Fatima]

You can use cd-hit with cutoff value of 1 (-c 1), it will remove the redundant sequences, but will also remove them if their header is different (as long as sequences are identical). So, for example when you have three identical sequences, the output file will only contain one of them ( the representative).

[Optimist]

Thank you

This might not help because most of the 16S rRNA seqs would be conserved (identical).

[Fatima]

Are the headers exactly identical?

grep "16S_rRNA" filename > headers 

#This is not very efficient but should do the work
cat headers | sort | uniq | while read line ; do grep -A 1 "${line}" filename >> output ; done