我怎样才能计算字母的频率
How can I count the frequency of letters
我有这样的数据
>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
但是,我想以更好的方式和更高效的方式为所有计算,尤其是当数据量很大时
不确定这会快多少,但如果您尝试,请post您的时间安排
$ awk '!/^>/ {n=split([=10=],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([=11=],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([=11=]); for(i=1;i<=n;i++) c[substr([=11=],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
快。可能是由于数组分配。
在任何 UNIX 机器上的任何 shell 中使用任何 awk:
$ cat tst.awk
!/^>/ {
for (i=1; i<=length([=10=]); i++) {
cnt[substr([=10=],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
使用 awk 可以轻松地对字符串中的字符进行计数。为此,您可以使用函数 gsub
:
gsub(ere, repl[, in])
Behave like sub
(see below), except that it shall replace all occurrences of the regular expression (like the ed utility global substitute) in [=20=]
or in the in
argument when specified.
sub(ere, repl[, in ])
Substitute the string repl
in place of the first instance of the extended regular expression ERE
in string in and return the number of substitutions. <snip> If in
is omitted, awk shall use the current record ([=20=]
) in its place.
source: 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=[=12=]; 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 [=13=] }
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 is based on Brian Kernighan's awk 记录在 "The AWK Programming Language",
作者:Al Aho、Brian Kernighan 和 Peter Weinberger
(Addison-Wesley, 1988, ISBN 0-201-07981-X)
。我不确定这个版本是否与 POSIX.
兼容
试试这个 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 2
给出了很好的比较,在 6.5 秒(在此系统上).
当然 none 这可以与优化的 bioawk
(C 代码?)程序进行比较。为此,我们需要用 C 语言编写(使用 Inline C
并不难)。
请注意,使用
为每个字符移除子调用(对 substr
)
if (not /^>/) { ++$kv{$_} for split //; }
结果 "only" 平均 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
我有这样的数据
>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
但是,我想以更好的方式和更高效的方式为所有计算,尤其是当数据量很大时
不确定这会快多少,但如果您尝试,请post您的时间安排
$ awk '!/^>/ {n=split([=10=],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([=11=],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([=11=]); for(i=1;i<=n;i++) c[substr([=11=],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
快。可能是由于数组分配。
在任何 UNIX 机器上的任何 shell 中使用任何 awk:
$ cat tst.awk
!/^>/ {
for (i=1; i<=length([=10=]); i++) {
cnt[substr([=10=],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
使用 awk 可以轻松地对字符串中的字符进行计数。为此,您可以使用函数 gsub
:
gsub(ere, repl[, in])
Behave likesub
(see below), except that it shall replace all occurrences of the regular expression (like the ed utility global substitute) in[=20=]
or in thein
argument when specified.
sub(ere, repl[, in ])
Substitute the stringrepl
in place of the first instance of the extended regular expressionERE
in string in and return the number of substitutions. <snip> Ifin
is omitted, awk shall use the current record ([=20=]
) in its place.source: 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=[=12=]; 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 [=13=] }
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 is based on Brian Kernighan's awk 记录在 "The AWK Programming Language", 作者:Al Aho、Brian Kernighan 和 Peter Weinberger (Addison-Wesley, 1988, ISBN 0-201-07981-X) 。我不确定这个版本是否与 POSIX.
兼容试试这个 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 2
给出了很好的比较,在 6.5 秒(在此系统上).
当然 none 这可以与优化的 bioawk
(C 代码?)程序进行比较。为此,我们需要用 C 语言编写(使用 Inline C
并不难)。
请注意,使用
为每个字符移除子调用(对substr
)
if (not /^>/) { ++$kv{$_} for split //; }
结果 "only" 平均 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