for loops to grep multiple texts from parent file to multiple files in single command

manpreet Best Answer 2 years ago

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?

manpreet 2 years ago

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.