How to loop multiple function in shell script?

Question

I need to extract sequences one after another consecutively from a large fasta files (multiple fasta files) and each extracted files to be saved in new fasta file (I mean the first sequence extracted from multiple fasta file to be saved as gene1 and second sequence extracted from multiple fasta file to be saved in another file and so on). For that, I tried a script like this,

for file in *.fasta do cat "$file" \
  | awk '/^>/ {printf("%s%s\t",(N>0?"\n":""),$0);N++;next;} \
          {printf("%s",$0);} \
     END {printf("\n");}' \
  | awk '{if ( NR==1 ) print $0;}' \
  | awk '{print $1"\n"$2}' | uniq > "new_${file}" 

But, its not serving my purpose. Firstly, I could not loop my sequence extraction procedure from multiple fasta file and save as gene1 and gene2 and so on.

For example,

I have fasta files as follows,

org1.fasta
>locus1 
ATGCGTAGAG
>LG12 
TATGCATGAT
>bpl
TAGATGGAT

org2.fasta
>ysr
ATGTAGCGA
>inc
TGATCTATG
>pym
ATGCATAGA

org3.fasta
>siv
TAGTAGTAT
>tdp
ATGAGTCGA
>tem
ATGCTHATG

I need to extract first positioned sequences from all the org1.fasta org2.fasta and org3.fasta files and save in a new file gene1.fasta and second positioned sequences as gene2.fastaand third positioned sequences to be saved as gene3.fasta.

Expected output,

gene1.fasta
>locus1 
ATGCGTAGAG
>ysr
ATGTAGCGA
>siv
TAGTAGTAT

gene2.fasta
>LG12 
TATGCATGAT
>inc
TGATCTATG
>tdp
ATGAGTCGA

gene3.fasta
>bpl
TAGATGGAT
>pym
ATGCATAGA
>tem
ATGCTHATG

Please help me to do the same. Thanks in advance.

Answer 1

You don't need a loop. You can do the whole thing with a simple awk one-liner:

awk -vRS='>' 'FNR>1{ printf ">%s",$0 > "gene"FNR-1".fasta"}' org*fasta

I ran this on your files and got:

$ for file in gene*; do echo "=== File: $file ==="; cat $file;  done
=== File: gene1.fasta ===
>locus1 
ATGCGTAGAG
>ysr
ATGTAGCGA
>siv
TAGTAGTAT
=== File: gene2.fasta ===
>LG12 
TATGCATGAT
>inc
TGATCTATG
>tdp
ATGAGTCGA
=== File: gene3.fasta ===
>bpl
TAGATGGAT
>pym
ATGCATAGA
>tem
ATGCTHATG

Explanation

• -vRS='>' : set the record separator to >. This will make awk consider every fasta sequence as a single "line".
• FNR>1{foo} : FNR is the line number of the current file. Because > is the first character of the entry, the first "line" will always be empty, so it needs t obe skipped.
• printf ">%s",$0 > "gene"FNR-1".fasta" : this does the work. It prints a > followed by the current "line" ($0, this is the entire sequence + the header) and redirects it into a file called gene, then the current line number -1 (-1 is needed because of the issue mentioned above) and a .fasta extension.

=====

Note that org*fasta will be expanded with the default collation order. This means that org10.fasta will be processed before org1.fasta:

$ printf "org1.fasta\norg10.fasta\n" | sort 
org10.fasta
org1.fasta

To avoid this and have the files in numerical order, you can use brace expansion instead of the *. So, if you know your files are all orgN.fasta where N can take values from 1 to 20, do this:

awk -vRS='>' 'FNR>1{ printf ">%s",$0 > "gene"FNR-1".fasta"}' org{1..20}.fasta

Answer 2

The title of your question indicates a shell solution, however you have added R and Python as tags. Here is an R solution making use of Biostrings pacakge (using two input files for simplicity):

library(Biostrings)

# create Biostring objects from each multifasta file
org1 <- readDNAStringSet("org1.fasta")
org2 <- readDNAStringSet("org2.fasta")

# loop over the Biostring objects
for(i in 1:length(org1)){

  output_file <- paste0("gene_", i, ".fasta")

  # create a Biostring object for each index
  genes_at_index_i <- c(org1[i], org2[i])
  assign(gene_name, genes_at_index_i)

  # save Biostring object as a multifasta file
  writeXStringSet(genes_at_index_i, output_file)
}

Output, gene 1:

>locus1 
ATGCGTAGAG
>ysr
ATGTAGCGA

Output, gene 2:

>LG12 
TATGCATGAT
>inc
TGATCTATG

Warning: The code above won't work if you don't have the same number of fasta entries in each multifasta file.