0
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

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
$\endgroup$
  • $\begingroup$ @K.Dineshkumar you're welcome! If one of the answers here solved your issue, please take a moment to accept it by clicking on the checkmark on the left. That will mark the question as answered and is the way that thanks are conveyed on the Stack Exchange sites. $\endgroup$ – terdon Nov 11 at 13:05
  • $\begingroup$ Thank you so much @terdon. However, I am facing a problem, this script perfectly fine for up to 9 number of fasta files. If it exceed 9, it takes 10th fasta file and printing first instead of 1st fasta file sequence. Same way, it prints 20th fasta file second instead of 2nd fasta file. I do not know how to solve this. $\endgroup$ – K. Dineshkumar Nov 11 at 13:43
  • 1
    $\begingroup$ @K.Dineshkumar ah yes, that is because of the order of collation. When you do echo org*fasta, it will sort alphabetically, so org10.fasta comes before org1.fasta (try with printf "org1.fasta\norg10.fasta\n" | sort). Try the updated answer. $\endgroup$ – terdon Nov 11 at 14:14
  • $\begingroup$ Many thanks @terdon for your detailed explanation. With your explanation, I learned how to play around the script, in order to serve our purpose. $\endgroup$ – K. Dineshkumar Nov 11 at 14:24
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thank you @haci for your help. However, I could not use this approach for 100 fasta files or above. Creating variable for each fasta file is tedious job. $\endgroup$ – K. Dineshkumar Nov 11 at 13:25
  • $\begingroup$ You can do so in a loop, i.e. creating a list of Biostring objects with for or lapply and then looping over that list. $\endgroup$ – haci Nov 11 at 13:30
  • $\begingroup$ Thank you @haci, I would try to do the same. $\endgroup$ – K. Dineshkumar Nov 11 at 13:43

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.