如何统计字母出现的频率

Question

我有一个这样的数据

>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

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

Answer 1

使用 awk 可以轻松计算字符串中的字符数。为此，您可以使用以下函数gsub：

gsub(ere, repl[, in]) 其行为类似于sub（见下文），但在指定时，它将替换参数中$0或参数中所有出现的正则表达式（如 ed 实用程序全局替换）。in

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

_{来源：Awk Posix 标准}

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

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

以下是基于此 fasta 文件的一些计时结果

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 的 awk ，该 awk 记录在Al Aho、Brian Kernighan 和 Peter Weinberger 所著的“The AWK 编程语言”中（Addison-Wesley，1988，ISBN 0-201-07981-X） 。我不确定这个版本是否与POSIX兼容。