Printing specific sequence and ID from combined fasta file using bash commands

Question

[SaltedPork]

I have a combined fasta file with all my sequences. I want to print the ID lines and their DNA sequence that end in .4

So far I have awk /'>*.4'/ {getline;print} combined.fasta Which prints the sequences that I want, how do I get the ID lines as well?

[Joe]

Can you show an example of your input data and ideal output?

[SaltedPork]

Thanks for responding, I have ID's that end in .1, .2 and so on up until .8. I just want the .4's and their sequence.

Input:

>16U035667-26_S26_L001.1
AGCTACGT
>16U035667-26_S26_L001.2
AGCTAACGTAC
>16U035667-26_S26_L001.4
ACGTACGTACTGAC

Output:

>16U035667-26_S26_L001.4
ACGTACGTACTGAC

Answer 1

[Joe]

Since you requested bash only:

#!/bin/bash

string="$2"

while read line ; do
    if [[ ${line:0:1} == ">" ]] ; then
        header="$line"
    else
        seq="$line"
        if [[ "$header" == *"$string" ]]; then
            echo -e "$header""\n""$seq"
        fi
    fi
done < $1

Put it in a script file and run it like:

$ bash parseheaders.sh seqs.fasta .4

Some of my own test data:

$ cat seqs.fasta
>tpg|Magnaporthiopsis_incrustans|JF414846
ACTGTAGTAGCTACGATCGATCAGATGATCACGTAGCATCGATCGATCATCGACTAGTAGATCACTCGACATAGATCCACATCAATAGATCATCATCATCATAATCGATCACTAGCAGC
>tpg|Pyricularia_pennisetigena|AB818016
GCAAGNTTCATGACGATGTAGAATGGCTTATCGAAGGGAGCAGGCCAGGGATTGAGGTCCGTCTCACGGGTTGGCTTCACTCCCCCACTGCCAGCCCTCTTGCTGCAACTCCACCAGAA
>tpg|Inocybe_sororia|EU525947
AACCANGCCGCGACGGCGGTGCGATCGGGAAACGCGGCGGTGGCGGAGGAATCGGCCATCCTTCACCATATCGGCCAAGGATTGTGGTTCCTGTAGGGCTCGCGCAGCCCAGGACGCGC

$ bash parseheaders.sh seqs.fasta 947
>tpg|Inocybe_sororia|EU525947
AACCANGCCGCGACGGCGGTGCGATCGGGAAACGCGGCGGTGGCGGAGGAATCGGCCATCCTTCACCATATCGGCCAAGGATTGTGGTTCCTGTAGGGCTCGCGCAGCCCAGGACGCGC

[SaltedPork]

This works very well, many thanks! PS. also didn't need to put quotes around the .4.

[GenoMax]

You can accept multiple answers (green check mark). So you should accept this one one as well.

[Joe]

Just as a final comment though, I'd advise following some of the other suggestions here and use proper parsers like bioperl, biopython, bioawk and so on. My own personal go-to script for pulling out sequences is here, using Biopython (though it finds the key of interest anywhere in the header, not just at the end, so it wasn't directly applicable to this case.

[Joe]

It shouldn't, but in case the . doesn't play nicely with your terminal, just enclose the string to look for in quotes: ".4"

Answer 2

[Pierre Lindenbaum]

grep '^>.*4$' -A 1 --no-group-separator in.fa

[SaltedPork]

Fantastic one-liner, many thanks!

[Joe]

I believe this will only work so long as your sequences are single lines FYI

Answer 3

[Sean Davis]

How about:

awk '/>*.4/ {print $0; getline; print}' combined.fasta

[SaltedPork]

hmmmm, I just get everything from combined.fasta printed to terminal.

[Sean Davis]

I'm on Mac OS. You may need to adjust the syntax a bit. Sorry....

[GenoMax]

Odd. Works fine on CentOS.

Answer 4

[Ram]

Check out bioawk - that should address your use case very well.