This is a test version of Biostars. For the public version, visit https://www.biostars.org.
extract sequences from fasta starting with a specific nucleotide

Hi, have a fasta file of 500 sequences, I want to extract the sequences starting with a specific nucleotide in the example A. Can anyone help?

E.g.

>seq_1
ACACACCGCTTCTAGCTG
>seq_2
ACAGGCAGAATTCTACAAGGA
>seq_3
CAAATATAATGACTATGGAATACC
>seq_4
CAATCGCCCGCTCACCTAGGTCT
>seq_5-493
TAACAGGCAGAATTCTACAAGGAC

Desired output:

>seq_1
ACACACCGCTTCTAGCTG
>seq_2
ACAGGCAGAATTCTACAAGGA

thank you in advance for your help

next-gen rna-seq sequence
• 4,124 views • link
written by Chris • Chris

Thank you all for your answers. Both methods are working!

• 1 views • link
written by Chris • Chris

2 answers

Something like following should work

grep '^A' -B 1 file.fa | sed '/--/d' > new_file.fa

Update: (Credits - Pierre)

grep '^A' -B 1 --no-group-separator file.fa > new_file.fa
• 1 views • link
written by venu venu

A yes, grep -B1 ! :-) , you know there is a secret option in grep to remove the double hyphen: --no-group-separator

• 1 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Aww..this is cool. Thank you Pierre.

• 1 views • link
written by venu venu

--no-group-separator is not available in all implementations of grep. So having a sed/grep -v to exclude separators is a safe bet.

• 1 views • link
written by GenoMax GenoMax

It is not there in man page also. (grep (GNU grep) 2.16).

• 1 views • link
written by venu venu

assuming there are only 2 lines per record:

cat input.fa |paste - - | awk  '($2 ~ /^A/)' | tr "\t" "\n"
• 1 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

What if I want to extract sequences with T at the 10th position? Thanks.

• 0 views • link
written by Seq225 Seq225

Log in to answer this question.

More posts like this