如何计算字母的频率答案

【问题标题】：How can I count the frequency of letters如何计算字母的频率
【发布时间】：2019-07-09 18:45:49
【问题描述】：

我有这样的数据

>sp|Q96A73|P33MX_HUMAN Putative monooxygenase p33MONOX OS=Homo sapiens OX=9606 GN=KIAA1191 PE=1 SV=1
RNDDDDTSVCLGTRQCSWFAGCTNRTWNSSAVPLIGLPNTQDYKWVDRNSGLTWSGNDTCLYSCQNQTKGLLYQLFRNLFCSYGLTEAHGKWRCADASITNDKGHDGHRTPTWWLTGSNLTLSVNNSGLFFLCGNGVYKGFPPKWSGRCGLGYLVPSLTRYLTLNASQITNLRSFIHKVTPHR
>sp|P13674|P4HA1_HUMAN Prolyl 4-hydroxylase subunit alpha-1 OS=Homo sapiens OX=9606 GN=P4HA1 PE=1 SV=2
VECCPNCRGTGMQIRIHQIGPGMVQQIQSVCMECQGHGERISPKDRCKSCNGRKIVREKKILEVHIDKGMKDGQKITFHGEGDQEPGLEPGDIIIVLDQKDHAVFTRRGEDLFMCMDIQLVEALCGFQKPISTLDNRTIVITSHPGQIVKHGDIKCVLNEGMPIYRRPYEKGRLIIEFKVNFPENGFLSPDKLSLLEKLLPERKEVEE
>sp|Q7Z4N8|P4HA3_HUMAN Prolyl 4-hydroxylase subunit alpha-3 OS=Homo sapiens OX=9606 GN=P4HA3 PE=1 SV=1
MTEQMTLRGTLKGHNGWVTQIATTPQFPDMILSASRDKTIIMWKLTRDETNYGIPQRALRGHSHFVSDVVISSDGQFALSGSWDGTLRLWDLTTGTTTRRFVGHTKDVLSVAFSSDNRQIVSGSRDKTIKLWNTLGVCKYTVQDESHSEWVSCVRFSPNSSNPIIVSCGWDKLVKVWNLANCKLK
>sp|P04637|P53_HUMAN Cellular tumor antigen p53 OS=Homo sapiens OX=9606 GN=TP53 PE=1 SV=4
IQVVSRCRLRHTEVLPAEEENDSLGADGTHGAGAMESAAGVLIKLFCVHTKALQDVQIRFQPQL
>sp|P10144|GRAB_HUMAN Granzyme B OS=Homo sapiens OX=9606 GN=GZMB PE=1 SV=2
MQPILLLLAFLLLPRADAGEIIGGHEAKPHSRPYMAYLMIWDQKSLKRCGGFLIRDDFVLTAAHCWGSSINVTLGAHNIKEQEPTQQFIPVKRPIPHPAYNPKNFSNDIMLLQLERKAKRTRAVQPLRLPSNKAQVKPGQTCSVAGWGQTAPLGKHSHTLQEVKMTVQEDRKCES
>sp|Q9UHX1|PUF60_HUMAN Poly(U)-binding-splicing factor PUF60 OS=Homo sapiens OX=9606 GN=PUF60 PE=1 SV=1
MGKDYYQTLGLARGASDEEIKRAYRRQALRYHPDKNKEPGAEEKFKEIAEAYDVLSDPRKREIFDRYGEEGLKGSGPSGGSGGGANGTSFSYTFHGDPHAMFAEFFGGRNPFDTFFGQRNGEEGMDIDDPFSGFPMGMGGFTNVNFGRSRSAQEPARKKQDPPVTHDLRVSLEEIYSGCTKKMKISHK
>sp|Q06416|P5F1B_HUMAN Putative POU domain, class 5, transcription factor 1B OS=Homo sapiens OX=9606 GN=POU5F1B PE=5 SV=2
IVVKGHSTCLSEGALSPDGTVLATASHDGYVKFWQIYIEGQDEPRCLHEWKPHDGRPLSCLLFCDNHKKQDPDVPFWRFLITGADQNRELKMWCTVSWTCLQTIRFSPDIFSSVSVPPSLKVCLDLSAEYLILSDVQRKVLYVMELLQNQEEGHACFSSISEFLLTHPVLSFGIQVVSRCRLRHTEVLPAEEENDSLGADGTHGAGAMESAAGVLIKLFCVHTKALQDVQIRFQPQLNPDVVAPLPTHTAHEDFTFGESRPELGSEGLGSAAHGSQPDLRRIVELPAPADFLSLSSETKPKLMTPDAFMTPSASLQQITASPSSSSSGSSSSSSSSSSSLTAVSAMSSTSAVDPSLTRPPEELTLSPKLQLDGSLTMSSSGSLQASPRGLLPGLLPAPADKLTPKGPGQVPTATSALSLELQEVEP
>sp|O14683|P5I11_HUMAN Tumor protein p53-inducible protein 11 OS=Homo sapiens OX=9606 GN=TP53I11 PE=1 SV=2
MIHNYMEHLERTKLHQLSGSDQLESTAHSRIRKERPISLGIFPLPAGDGLLTPDAQKGGETPGSEQWKFQELSQPRSHTSLKVSNSPEPQKAVEQEDELSDVSQGGSKATTPASTANSDVATIPTDTPLKEENEGFVKVTDAPNKSEISKHIEVQVAQETRNVSTGSAENEEKSEVQAIIESTPELDMDKDLSGYKGSSTPTKGIENKAFDRNTESLFEELSSAGSGLIGDVDEGADLLGMGREVENLILENTQLLETKNALNIVKNDLIAKVDELTCEKDVLQGELEAVKQAKLKLEEKNRELEEELRKARAEAEDARQKAKDDDDSDIPTAQRKRFTRVEMARVLMERNQYKERLMELQEAVRWTEMIRASRENPAMQEKKRSSIWQFFSRLFSSSSNTTKKPEPPVNLKYNAPTSHVTPSVK

我想计算每个字母有多少，所以如果我有一个，我会这样数

>sp|Q96A73|P33MX_HUMAN Putative monooxygenase p33MONOX OS=Homo sapiens OX=9606 GN=KIAA1191 PE=1 SV=1
RNDDDDTSVCLGTRQCSWFAGCTNRTWNSSAVPLIGLPNTQDYKWVDRNSGLTWSGNDTCLYSCQNQTKGLLYQLFRNLFCSYGLTEAHGKWRCADASITNDKGHDGHRTPTWWLTGSNLTLSVNNSGLFFLCGNGVYKGFPPKWSGRCGLGYLVPSLTRYLTLNASQITNLRSFIHKVTPHR

cat input.txt | grep -v ">" | fold -w1 | sort | uniq -c



   6 A
   9 C
  10 D
   1 E
   7 F
  18 G
   5 H
   4 I
   7 K
  21 L
  15 N
   7 P
   6 Q
  11 R
  16 S
  18 T
   7 V
   8 W
   7 Y

但是，我想以更好的方式和更高效的方式为所有人进行计算，尤其是当数据量很大时

【问题讨论】：

您是逐行报告计数还是为整个文件报告一次？
@karakfa 获取整个文件
Bash script to find the frequency of every letter in a file的可能重复

标签： awk sed bioinformatics fasta

【解决方案1】：

用 awk 可以很容易地计算字符串中的字符。为此，您可以使用函数gsub：

gsub(ere, repl[, in]) 表现得像sub（见下文），除了它应替换$0 或in 参数中所有出现的正则表达式（如ed 实用程序全局替换）。

sub(ere, repl[, in ]) 用字符串repl 代替字符串中扩展正则表达式ERE 的第一个实例并返回替换次数。如果in 被省略，awk 将使用当前记录 ($0) 代替它。

_{来源：Awk Posix Standard}

以下两个函数以这种方式执行计数：

function countCharacters(str) {
    while(str != "") { c=substr(str,1,1); a[toupper[c]]+=gsub(c,"",str) }
}

或者如果可能出现很多相等的连续字符，以下解决方案可能会缩短几秒钟。

function countCharacters2(str) {
    n=length(str)
    while(str != "") { c=substr(str,1,1); gsub(c"+","",str);
       m=length(str); a[toupper[c]]+=n-m; n=m
    }
}

您可以在下面找到基于第一个函数的 4 个实现。前两个在标准 awk 上运行，后两个在 fasta 文件的优化版本上运行：

1.读取序列并逐行处理：

awk '!/^>/{s=$0; while(s!="") { c=substr(s,1,1); a[c]+=gsub(c,"",s) } }
     END {for(c in a) print c,a[c]}' file

2。连接所有序列并在最后进行处理：

awk '!/^>/{s=s $0 }
     END {while(s!="") { c=substr(s,1,1); a[c]+=gsub(c,"",s) }
         for(c in a) print c,a[c]}' file

3.与 1 相同，但使用 bioawk:

bioawk -c fastx '{while ($seq!=""){ c=substr($seq,1,1);a[c]+=gsub(c,"",$seq) } }
                 END{ for(c in a) print c,a[c] }' file

4.与 2 相同，但使用 bioawk:

bioawk -c fastx '{s=s $seq}
                 END { while(s!="") { c=substr(s,1,1); a[c]+=gsub(c,"",s) }
                       for(c in a) print c,a[c]}' file

这里有一些基于this fasta-file的计时结果

OP            : grep,sort,uniq : 47.548 s
EdMorton 1    : awk            : 39.992 s
EdMorton 2    : awk,sort,uniq  : 53.965 s
kvantour 1    : awk            : 18.661 s
kvantour 2    : awk            :  9.309 s
kvantour 3    : bioawk         :  1.838 s
kvantour 4    : bioawk         :  1.838 s
karafka       : awk            : 38.139 s
stack0114106 1: perl           : 22.754 s
stack0114106 2: perl           : 13.648 s
stack0114106 3: perl (zdim)    :  7.759 s

注意： BioAwk 基于Brian Kernighan's awk，记录在"The AWK Programming Language", by Al Aho, Brian Kernighan, and Peter Weinberger (Addison-Wesley, 1988, ISBN 0-201-07981-X) 中。我不确定这个版本是否兼容POSIX。

【讨论】：

谢谢。不幸的是，这并不是真正的最佳解决方案，这比 O(n) 更像是一个 O(n^2) 解决方案（由于 gsub）
gsub 计数是聪明的主意。
哇，这个bioawk 真的很高效。感谢所有这些:)
我突然想到，1 将 $0 保存在字符串中，然后 match()/gsub()-ing 可能比直接在 $0 上操作更有效，因为每次更改时$0 您可能会再次调用字段拆分（取决于您的 awk 实现 - 有些仅在脚本中明确提及字段时才进行字段拆分）。
@EdMorton 我刚刚进行了测试。我没有看到时间上的差异，但你是对的。通过一些 awk 实现，这将有所不同。（更新了实现）

【解决方案2】：

试试这个 Perl 解决方案以获得更好的性能。

$ perl -lne ' 
            if( ! /^>/ ) { while(/./g) { $kv{$&}++} }  
        END { for ( sort keys %kv) { print "$_ $kv{$_}" }} 
' learner.txt
A 107
C 41
D 102
E 132
F 65
G 140
H 52
I 84
K 114
L 174
M 39
N 67
P 107
Q 88
R 101
S 168
T 115
V 101
W 27
Y 30

$

另一个使用 Perl 的解决方案，针对性能进行了优化。

$ time perl -lne ' 
     if( ! /^>/ ) { for($i=0;$i<length($_);$i++) 
     { $x=substr($_,$i,1); $kv{$x}++ } }  
    END { for ( sort keys %kv) { print "$_ $kv{$_}" }} 
' chrY.fa
A 2994088
C 1876822
G 1889305
N 30812232
T 3002884
a 4892104
c 3408967
g 3397589
n 140
t 4953284

real    0m15.922s
user    0m15.750s
sys     0m0.108s

$

编辑以进一步优化性能

下面报告的所有时间都是在桌面上运行 3-5 次的平均值，大约在同一时间完成，但为了避免明显的缓存效应而交换了时间。

将 C 样式的 for 循环更改为 for my $i (0..length($_)) 将第二个解决方案的速度从 9.2 秒缩短到 6.8 秒。

然后，在每个操作中删除一个标量 ($x)，使用

if (not /^>/) { for $i (0..length($_)) { ++$kv{ substr($_,$i,1) } } }

将此速度提高到 5.3 秒。

通过复制$_ 进一步减少变量的使用，从而释放循环以使用$_

if (not /^>/) { $l=$_; ++$kv{ substr($l,$_,1) } for 0..length($l) }

只有一点帮助，运行时间 5.2 秒。

这与awk 解决方案（在kvantour answer 中给出的kvantour 2 进行了很好的比较，在6.5 秒（在此系统上）。

当然，这一切都无法与优化的bioawk（C 代码？）程序相提并论。为此，我们需要用 C 语言编写它（使用 Inline C 并不难）。

请注意，使用删除每个字符的子调用（到substr）

if (not /^>/) { ++$kv{$_} for split //; }

平均结果“只有”6.4 秒，不如上述调整好；这是一个惊喜。

这些时间是在 v5.16 的桌面上。在 v5.24 上，在同一台机器上，最佳情况（substr 循环中没有额外变量）时间是 4.8 秒，而没有substr（但有split）的时间是 5.8 秒。很高兴看到新版本的 Perl 性能更好，至少在这些情况下是这样。

供其他人参考和方便计时，完整代码以获得最佳运行

time perl -lne'
    if (not /^>/) { $l=$_; ++$kv{ substr($l,$_,1) } for 0..length($l) }
    END { for ( sort keys %kv) { print "$_ $kv{$_}" }}
' chrY.fa

【讨论】：

将您的解决方案添加到时间表中：22.745 秒
谢谢@kvantour..我通过优化添加了一个解决方案..在我的机器上你的大文件花了 15.92 秒..bioawk 太棒了！
添加了您的第二个解决方案（13 秒）
@kvantour 注意进一步的改进。
@stack0114106 我冒昧地在您的答案中添加了一个部分。如果这不符合您的喜好，我深表歉意，请告诉我，我会将其拆分为单独的答案。我认为它最适合这里，因为它只是调整您的代码。

【解决方案3】：

不确定这会快多少，但如果您尝试，请发布您的时间

$ awk '!/^>/ {n=split($0,a,""); for(i=1;i<=n;i++) c[a[i]]++} 
       END   {for(k in c) print k,c[k]}' file | sort

A 6
C 9
D 10
E 1
F 7
G 18
H 5
I 4
K 7
L 21
N 15
P 7
Q 6
R 11
S 16
T 18
V 7
W 8
Y 7

这会报告文件的计数，而不是逐行报告。如下所述，并非所有awk 都支持空字符串拆分。

以下是三种方法的时间安排：

$ time grep -v ">" filey | fold -w1 | sort | uniq -c >/dev/null

real    0m11.470s
user    0m11.746s
sys     0m0.260s

$ time awk '{n=split($0,a,""); for(i=1;i<=n;i++) c[a[i]++]} END{for(k in c) print k,c[k]}' filey >/dev/null 

real    0m7.441s
user    0m7.334s
sys     0m0.060s

$ time awk '{n=length($0); for(i=1;i<=n;i++) c[substr($0,i,1)]++} END{for(k in c) print k,c[k]}' filey >/dev/null

real    0m5.055s
user    0m4.979s
sys     0m0.047s

用于测试文件

$ wc filey

  118098   649539 16828965 filey

substr 比 split 快让我感到惊讶。可能是由于数组分配。

【讨论】：

您应该提到，根据 POSIX，在空字符上拆分是未定义的行为，因此它仅适用于某些 awk，例如GNU awk。
@Tiw 你误解了我的评论。在split($0,a,"") 中，问题不在于$0 是一个空字符串，而是"" 是一个空字符串。使用空字符串作为“正则表达式”进行拆分是每个 POSIX 未定义的行为。无论是 FS 的值还是 split() 的第三个参数都是如此。
哦，你的意思是……很抱歉我误会了。在 Freebsd 的 awk 上尝试过，但效果很好。 @埃德莫顿
@Tiw 是的，有些 awks 会做你想做的事，但其他人不会。请参阅 the POSIX awk spec - 在讨论 FS 时 If FS is a null string, the behavior is unspecified. 和稍后在 split(...,fs) 概要中 The effect of a null string as the value of fs is unspecified.。
还在我的回答中添加了一些时间。

【解决方案4】：

在任何 UNIX 机器上的任何 shell 中使用任何 awk：

$ cat tst.awk
!/^>/ {
    for (i=1; i<=length($0); i++) {
        cnt[substr($0,i,1)]++
    }
}
END {
    for (char in cnt) {
        print char, cnt[char]
    }
}

$ awk -f tst.awk file
A 107
N 67
P 107
C 41
Q 88
D 102
E 132
R 101
F 65
S 168
G 140
T 115
H 52
I 84
V 101
W 27
K 114
Y 30
L 174
M 39

或者如果您愿意：

$ awk -v ORS= '!/^>/{gsub(/./,"&\n"); print}' file | sort | uniq -c
    107 A
     41 C
    102 D
    132 E
     65 F
    140 G
     52 H
     84 I
    114 K
    174 L
     39 M
     67 N
    107 P
     88 Q
    101 R
    168 S
    115 T
    101 V
     27 W
     30 Y

【讨论】：

@Ed Morton awk: syntax error at source line 1 source file tst.awk context is cat >>> tst. <<< awk awk: bailing out at source line 11
@Learner 再看看我的回答，你没有做它所显示的。似乎您可能已将 cat 行复制到 awk 脚本中。
当我看到这个问题时，我正在考虑这个解决方案，但我认为多个substr 将比只拆分一次效率低得多。但是我不确定这一点，这次我想对了吗？
@Tiw 在空字符上拆分不可移植 - 它会在某些 awk 中拆分为字符，在其他行中保持原样，并且可以做任何其他事情并且仍然符合 POSIX因为那是未定义的行为。我真的不知道单个 split() 是否比多个 substr() 更有效，但恕我直言，在这种情况下可移植性更重要。