This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to extract fasta headers that contains specific text and it's sequence

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

barrnap 16srrna fasta
• 4,285 views • link
updated by Neel • Neel • written by Optimist Optimist

What have you tried so far?

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

• 0 views • link
written by lieven.sterck lieven.sterck

filter_fasta.py has worked.

• 0 views • link
written by Optimist Optimist

with awk and flattend fasta:

awk '/^>16S/{print;getline;print}' seq.fa
• 0 views • link
written by cpad0112 cpad0112

Hi optimist,

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

Thank you so much!

• 0 views • link
written by Neel • Neel

3 answers

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!!!

• 0 views • link
updated by Ram Ram • written by Optimist Optimist

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

• 0 views • link
written by cpad0112 cpad0112

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}'
• 0 views • link
written by 5heikki 5heikki

If each sequence is one and only one line: 

grep -A 1 "16S_rRNA:" filename

Should do the trick!

• 0 views • link
written by Fatima Fatima

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).

• 0 views • link
written by Optimist Optimist

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).

• 0 views • link
written by Fatima Fatima

Thank you

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

• 0 views • link
written by Optimist Optimist

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
• 0 views • link
written by Fatima Fatima

Log in to answer this question.

More posts like this