This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Sort the fasta sequences based on the desired amino acid content

Dear all, I have a fasta file having multiple protein sequences. I want to sort those protein sequence based on amount of desired amino acids. How can I do that?

For example: I want the following sequences in descending order with D and G amino acid content:

>FastaA
ASDFGHILMNV
>FastaB
SKSYGLKQAPPDTITLIAAKSNS
>FastaC
FQRRYVVWILAVSRHIVFLEN
>FastaD
LAPKDYKLELDDGSDVMK

Output file:

>FastaD
LAPKDYKLELDDGSDVMK
>FastaB
SKSYGLKQAPPDTITLIAAKSNS
>FastaA
ASDFGHILMNV
>FastaC
FQRRYVVWILAVSRHIVFLEN
fasta protein seqio biophython
• 1,356 views • link
updated by cpad0112 cpad0112 • written by Pradeep • Pradeep

What exactly do you mean by "amount of desired amino acids"? In any case, if you can code in python, you can easily sort strings (sequences) by a custom function using the sort(key=your_function) syntax.

• 0 views • link
written by liorglic liorglic

Assuming that sequences are in single line

$ awk -v OFS="\t" '/^>/ {getline seq} {print $0, seq, length(seq),gsub("[DG]","",seq)}' test.fa   |  sort -k4rV,4 -k3rV,3 | awk '{print $1"\n"$2}'    

>FastaD
LAPKDYKLELDDGSDVMK
>FastaB
SKSYGLKQAPPDTITLIAAKSNS
>FastaA
ASDFGHILMNV
>FastaC
FQRRYVVWILAVSRHIVFLEN

If there is a tie, longer sequence gets printed first

• 0 views • link
written by cpad0112 cpad0112

1 answer

linearize, create a new column with the number of character, sort on that column, restore the fasta:

awk '/^>/ {printf("%s%s\t",(N>0?"\n":""),$0);N++;next;} {printf("%s",$0);} END {printf("\n");}' in.fa |\
awk -F '\t' '{S=$2;gsub(/[^GC]/,"",S);printf("%d\t%s\t%s\n",length(S),$1,$2);}' |\
sort -t $'\t' -k1,1n |\
cut -f 2- |\
tr "\t" "\n"

>FastaC
FQRRYVVWILAVSRHIVFLEN
>FastaA
ASDFGHILMNV
>FastaB
SKSYGLKQAPPDTITLIAAKSNS
>FastaD
LAPKDYKLELDDGSDVMK
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Log in to answer this question.

More posts like this