☰

如何打印文件名和fasta序列总数？

tag-icon • tag-icon shell-script python bioinformatics
如何打印文件名和fasta序列总数？

我有一个 fasta 文件，即 test.fasta、pas.fasta、cel.fasta 如下所示

test.fasta
>tile
ATGTC
>259
TGAT

pas.fasta
>ta
ATGCT

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

我需要打印文件名和 fasta 序列总数，如下所示，

test,2
pas,1
cel,3

我使用了以下命令但未能达到我的目的

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

请帮助我做同样的事情。

提前致谢。

答案1

这就是(to ount)-c选项的用途：grepc

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

请注意，如果只有一个匹配文件，则不会打印文件名。某些grep实现可以-H选择强制始终打印文件名：

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

要获得确切的预期输出，您只需替换.fasta:为,：

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

（这里假设您的文件名不包含其他出现的情况，.fasta:例如my.fasta:foo.fasta; 当然，如果输出采用 CSV 格式，文件名中的换行符或,或字符以及潜在的空白字符也会成为问题）"

答案2

使用 GNU awk 来表示ENDFILE和gensub()：

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

无论存在 1 个还是多个输入文件（只要它小于 shell args 限制），并且无论文件名中存在哪些字符，只要不存在=（当 awk 可以时，在某些情况下），上述内容都将起作用，将其视为变量赋值 - 如果会发生这种情况，请参阅https://www.gnu.org/software/gawk/manual/gawk.html#Other-Arguments了解如何处理）。

您可以对任何 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

答案3

使用 Raku（以前称为 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

输入示例（密码）：

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

示例输出：

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

该解决方案使用 Raku，Perl 编程语言家族的成员。上面的代码依赖于 Raku 的dir()和grep()例程，因此在经典 shell 介导的文件通配不存在或受到限制的平台上可能很有用（请参阅 SO 讨论这里）。

简而言之，raku使用-e选项调用，它告诉 Raku 的编译器 (Rakudo) 编译并执行给定的单行程序。dir()使用关键字调用该方法for，该关键字告诉 Raku 循环使用过滤器获得的文件名值test => / \.fasta $/。由此找到的文件名被单独加载到$fh变量中并在{…}块中进行分析。

在块内部，从右到左读取，lines是从$fhfilehandle 变量中提取的，该变量会自动加载到$_topic 变量中。这些行（现在保存在 中$_）grep通过 -ped 来查找具有^^start-of-line的行>，表示新序列的开始fasta。这样找到的行使用 进行计数elems，并且$fh和 计数使用逗号进行join分隔,并使用 逐行返回put。

OP最初要求返回“文件名”，因此从技术上讲，上面的代码就足够了。然而，因为$fh是IO::Path对象，所以删除扩展名很容易.fasta：只需附加extension例程即可$fh识别.fasta扩展名，然后将其替换为空（""）：

~$ 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

相关内容