1

I have 29 fasta files (.fa as extension) named and stored sequences according to their genes.

(Example: ribosomal protein L1, ribosomal protein L6P/L9E,...)

There were a total of 722 species existed among all those 29 fasta files. Each sequence has their genes and species name labelled in the first line and second line was filled with its sequence.

There will be more than 1 gene sequence for 1 species.

I want to transfer the 722 species from 29 fasta files sorted according to genes into separate 722 files (sort them under species instead of genes).

The name of the species in the parent file were enclosed by square brakets [ ].

How can I use for loops to extract the 722 files and name the files according to its sequence name?

Example from Ribosomal Protein L1.fa:

>gi|103486926|ref|YP_616487.1| 50S ribosomal protein L1 [Sphingopyxis alaskensis RB2256]
MAKLTKKQKALEGKVDAQKLHGVDEAIKLVRELATAKFDETLEIAMNLGVDPRHADQMVRGVVTLPAGTGKDVKVAVFAR

Example from Ribosomal Protein L6PL9E.fa:

>gi|410479108|ref|YP_006766745.1| ribosomal protein L6P/L9E [Leptospirillum ferriphilum ML-04]
MGFTHTVEFTLPSLIKASIEKQTIITLSSPDKELLGQFAADVRSIRPPEPYKGKGIKYSGEKILRKEGKTGKK

For the first example,

Species name: Sphingopyxis alaskensis RB2256

Gene sequence: MAKLTKKQKALEGKVDAQKLHGVDEAIKLVRELATAKFDETLEIAMNLGVDPRHADQMVRGVVTLPAGTGKDVKVAVFA

I wanted to name the file as Sphingopyxis alaskensis RB2256.fa and insert all sequences with this species name into this file.

I am using bash shell to do this. I can use grep to have things done:

grep -A+1 "Sphingopyxis alaskensis RB2256" *.fa >> Sphingopyxis alaskensis RB2256.fa

But I will need to do it 722 times to get my sequences sorted according to species.

Is grep in for loops can be used to simplify the work? Or there are alternative methods to do so?

  • 1
    @muru You had no way of knowing this, but the > is part of the file. It's a requirement for the fasta format. – terdon Apr 5 at 10:30
  • why are you assuming your sequences will always be just a single line? Also, while this is perfectly on topic here, you might want to use our sister site, Bioinformatics in the future to avoid this sort of confusion we had here with the > and so your question will be asked of people who are familiar with the file formats you deal with. – terdon Apr 5 at 10:34
  • @terdon I am new to bioinformatics as well as scripting. I appreciate the site recommendation! Also, thanks for notifying about error on my sequences. I made some mistakes in extracting the sequences. They are supposed to be more than one line. – web Apr 8 at 0:51
2

The Fasta format doesn't require that all sequences be on a single line. In fact, that isn't even common, since most biological sequences are long. So your grep will fail in any cases where there's more than one line of sequence for the ID. Also, your grep command will create a file called Sphingopyxis and not a file called Sphingopyxis alaskensis RB2256.fa.

In any case, you can do something like this to get each sequence into a file names after the species:

awk -F'[][]' '/>/{n=$2}; {print >> n".fa"}' *.fa 

However, I strongly urge you not to use spaces in your file names since that will only make your life harder. A safer approach would be:

awk -F'[][]' '/>/{n=$2; gsub(/ /,"_",n)}; {print >> n".fa"}' *.fa 

The gsub replaces all spaces in the species name with _, resulting in these files:

Leptospirillum_ferriphilum_ML-04.fa  Sphingopyxis_alaskensis_RB2256.fa

Note that both approaches above can deal with multi-line sequences.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.