0

I have a fasta file namely test.fasta, pas.fasta, cel.fasta as shown below

test.fasta
>tile
ATGTC
>259
TGAT

pas.fasta
>ta
ATGCT

cel.fasta
>787
TGTAG
>yog
TGTAT
>In
NNTAG

I need to print the file name and the total number of fasta sequences as shown below,

test,2
pas,1
cel,3

I have used the following commands but failed to serve my purpose

grep ">" test.fasta | wc -l && ls test.fasta

Please help me to do the same.

Thanks in advance.

4

3 Answers 3

6

That's what the -c option of grep (to count) is for:

$ grep -ce '^>' -- *.fasta
cel.fasta:3
pas.fasta:1
test.fasta:2

Note that if there's only one matching file, the file name will not be printed. Some grep implementations have a -H option to force the file name to be printed always:

$ grep -Hce '^>' -- *.fasta
cel.fasta:3

To get your exact expected output, you just need to replace .fasta: with ,:

$ grep -Hce '^>' -- *.fasta | sed 's/\.fasta:/,/'
cel,3
pas,1
test,2

(here assuming your file names don't contain other occurrences of .fasta: such as my.fasta:foo.fasta; of course newline or , or " characters and potentially whitespace characters in file names would also be a problem if the output is meant to be in CSV format)

3

Using GNU awk for ENDFILE and gensub():

$ awk '/^>/{c++} ENDFILE{print gensub(/\.[^.]*$/,",",1,FILENAME) c+0; c=0}' *.fasta
cel,3
pas,1
test,2

The above will work no matter if it's 1 or many input files present (as long as it's less than the shell args limit) and no matter which characters are present in your file names as long as it's not = (when awk could, in some cases, treat it as a variable assignment - if that can happen then see https://www.gnu.org/software/gawk/manual/gawk.html#Other-Arguments for how to handle it).

You can do the same with any awk:

$ awk '/^>/{c[FILENAME]++} END{for (i=1; i<ARGC; i++) {f=ARGV[i]; sub(/\.[^.]*$/,",",f); print f c[ARGV[i]]+0} }' *.fasta
cel,3
pas,1
test,2
5
  • 1
    That doesn't work for filenames that contain = characters (and what's on the left side of the first = has the shape of a valid awk variable name). You'd want to use ./*.fasta and strip the ./ in your ENDFILE statement or use -E /dev/null and put the awk code in a -e <expression> Sep 5, 2021 at 15:16
  • @StéphaneChazelas ah, good point. I'll tweak the answer to account for that thx.
    – Ed Morton
    Sep 5, 2021 at 15:18
  • 1
    While you're at using awk, you could also encode the filename in CSV format (wrap in "...", escape "s as ""). Sep 5, 2021 at 15:22
  • All true. I'll add that if the OP asks for it.
    – Ed Morton
    Sep 5, 2021 at 15:24
  • 1
    Thank you @StéphaneChazelas and Ed Morton for your detailed explanation and valuable code. Both codes works well for me.
    – Kumar
    Sep 6, 2021 at 8:37
1

Using Raku (formerly known as Perl_6)

~$ raku -e 'for dir(test => / \.fasta $/ ) -> $fh { put ($fh,  $_.grep(/ ^^ \> /).elems).join(",") given $fh.lines() };' > outfile.txt

#OR

~$ raku -e 'for dir(test => / \.fasta $/ ) -> $fh { put join ",", $fh,  $_.grep(/ ^^ \> /).elems given $fh.lines() };'  > outfile.txt

Sample Input (pwd):

~$ ls *.fasta
cel.fasta   pas.fasta   test.fasta

Sample Output:

cel.fasta,3
pas.fasta,1
test.fasta,2

This solution uses Raku, a member of the Perl-family of programming languages. The code above relies on Raku's dir() and grep() routines, and therefore may be useful on platforms where classic shell-mediated file-globbing is absent or restricted (see SO discussion here).

Briefly, raku is called with the -e option, which tells Raku's compiler (Rakudo) to compile and execute the given one-liner program. The dir() method is called with the for keyword which tells Raku to loop through the filename values obtained using the filter test => / \.fasta $/. Filenames thus found are loaded individually into $fh variable and analyzed in a {…} block.

Inside the block, reading from right-to-left, lines are extracted from the $fh filehandle variable, which automatically loads into the $_ topic variable. These lines (now held within $_) are grep-ped through to find lines with a ^^ start-of-line >, signifying the start of a new fasta sequence. Lines thus found are counted using elems, and the $fh and count are join-ed with a , comma and returned linewise using put.

The OP originally asks for the "file name" to be returned, so technically the code above suffices. However because $fh are IO::Path objects, it's easy enough to remove the .fasta extension: just append the extension routine to $fh to recognize the .fasta extension, and replace it with nothing (""):

~$ raku -e 'for dir(test => / \.fasta $/ ) -> $fh { put ($fh.extension(""), $_.grep(/ ^^ \> /).elems).join(",") given $fh.lines() };'
cel,3
pas,1
test,2

https://docs.raku.org/routine/dir
https://docs.raku.org/routine/grep
https://docs.raku.org/routine/extension
https://raku.org

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .