如何统计字符串中n个长度的字符组合出现的次数
How to count the number of occurences of n-length combinations of characters in string
我正在使用下面的一行来列出 ATCG
组合的出现次数,形成长度为 6 的字符串。除了不打印 0 次匹配之外,它工作正常。有没有办法将正则表达式或其他部分更改为打印“0 ATTTAG”之类的内容?
#!/bin/bash
for file in e_coli.fa
do
base=$(basename $file .fa)
cat $file | perl -nE 'say for /(?<=([ATCG]{6}))/g' \
| sort | uniq -c >> ${base}_hexhits_6mer.txt
done
stdout:
465 AAAAAA
607 AAAAAC
661 AAAAAG
581 AAAAAT
563 AAAACA
807 AAAACC
770 AAAACG
373 AAAACT
663 AAAAGA
1213 AAAAGC
由于uniq -c
计算一行出现的次数,所以不可能return0。请求的更改需要完全重写。
perl -e'
while (<>) {
++$counts{$_} for /(?=([ATCG]{6}))/g;
}
for my $seq (glob("{A,C,G,T}" x 6)) {
printf("%7d %s\n", $counts{$seq}, $seq);
}
' "$file" >"${base}_hexhits_6mer.txt"
您要执行的操作要复杂得多。要了解您没有看到的内容,您首先需要了解所有可能的字符组合,然后您可以根据这些组合进行筛选。
在这里,我使用 Perl 中的 substr 的滑动 window 方法来查找 As
字符串中的所有 "seen" ATCG
个字符,散列中的 Ts
、Cs
和 Gs
(从 __DATA__
读取)。然后对这些进行排序,以便首先显示最常见的 6 聚体,然后打印出来。
use strict;
use warnings;
my @bases = qw/ A G C T /;
my %data;
for my $a1(@bases){
for my $a2(@bases){
for my $a3(@bases){
for my $a4(@bases){
for my $a5(@bases){
for my $a6(@bases){
$data{"$a1$a2$a3$a4$a5$a6"} = 0;
}
}
}
}
}
}
my $nucs = <DATA>;
my $len = length($nucs);
for (my $i = 0; $i <= $len - 6; $i++) {
my $kmer = substr($nucs, $i, 6);
next if $kmer =~ tr/ACGT//c;
$data{$kmer}++; # populate hash with "seen" 6-mers
}
# print out sorted hash
foreach my $seq (sort { $data{$b} <=> $data{$a} } keys %data ){
print "$seq,$data{$seq}\n";
}
__DATA__
ATGCCCGTCGTAGTCATGCATGCATCGATCGATGCATGCTACGTGTTGT
显然会有一种 better/prettier 方法来计算字符串中字符的所有排列,而不是我所做的,但它确实有效。
正如 Borodin 所说,这主要打印出 "unseen" 字符串的变体。
最简单的方法是为每个模式构建一个出现次数的散列,然后打印所有可能模式的次数
此程序使用 glob
技巧生成由 A、T、C 和 G 组成的所有可能的六字符字符串列表
use strict;
use warnings 'all';
my @files = qw/ e_coli.fa /;
my %counts;
for my $file ( @files ) {
open my $fh, '<', $file or die qq{Unable to open "$file" for input: $!};
while ( <$fh> ) {
++$counts{} while /(?= ( [ATCG]{6} ) ) /gx;
}
}
for my $pattern ( glob '{A,T,C,G}' x 6 ) {
printf "%4d %s\n", $counts{$pattern} // 0, $pattern;
}
如果您有大量数据并且需要更快的速度,这里有一个 C 解决方案:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
void reader(FILE* in, unsigned long hist[4096]) {
for (unsigned long key=0, count=0;;) {
switch(getc(in)) {
case EOF: return;
case 'A': key <<= 2; break;
case 'C': key <<= 2; key += 1; break;
case 'G': key <<= 2; key += 2; break;
case 'T': key <<= 2; key += 3; break;
default: count=0; continue;
}
if (count == 5) ++hist[key & 0xFFF];
else ++count;
}
}
int putkey(FILE* out, unsigned long key) {
char s[6];
for (int j=6; j--; key >>= 2) s[j] = "ACGT"[key&3];
return fprintf(out, "%.6s", s);
}
void writer(FILE* out, unsigned long hist[4096]) {
for (unsigned long key = 0; key < 4096; ++key) {
fprintf(stdout, "%7lu ", hist[key]);
putkey(out, key);
putchar('\n');
}
}
int main(int argc, char** argv) {
FILE* in = stdin;
if (argc > 1) in = fopen(argv[1], "r");
if (!in) { perror(argv[1]); exit(1); }
unsigned long hist[4096] = {0};
reader(in, hist);
writer(stdout, hist);
return 0;
}
处理一个 31MB 的 fastq 样本(碰巧包括所有 4096 个可能的六字符序列)只花了不到半秒的时间; Perl 解决方案分别用了 12 秒 (fugu) 和 18 秒 (ikegami/borodin)。
我正在使用下面的一行来列出 ATCG
组合的出现次数,形成长度为 6 的字符串。除了不打印 0 次匹配之外,它工作正常。有没有办法将正则表达式或其他部分更改为打印“0 ATTTAG”之类的内容?
#!/bin/bash
for file in e_coli.fa
do
base=$(basename $file .fa)
cat $file | perl -nE 'say for /(?<=([ATCG]{6}))/g' \
| sort | uniq -c >> ${base}_hexhits_6mer.txt
done
stdout:
465 AAAAAA
607 AAAAAC
661 AAAAAG
581 AAAAAT
563 AAAACA
807 AAAACC
770 AAAACG
373 AAAACT
663 AAAAGA
1213 AAAAGC
由于uniq -c
计算一行出现的次数,所以不可能return0。请求的更改需要完全重写。
perl -e'
while (<>) {
++$counts{$_} for /(?=([ATCG]{6}))/g;
}
for my $seq (glob("{A,C,G,T}" x 6)) {
printf("%7d %s\n", $counts{$seq}, $seq);
}
' "$file" >"${base}_hexhits_6mer.txt"
您要执行的操作要复杂得多。要了解您没有看到的内容,您首先需要了解所有可能的字符组合,然后您可以根据这些组合进行筛选。
在这里,我使用 Perl 中的 substr 的滑动 window 方法来查找 As
字符串中的所有 "seen" ATCG
个字符,散列中的 Ts
、Cs
和 Gs
(从 __DATA__
读取)。然后对这些进行排序,以便首先显示最常见的 6 聚体,然后打印出来。
use strict;
use warnings;
my @bases = qw/ A G C T /;
my %data;
for my $a1(@bases){
for my $a2(@bases){
for my $a3(@bases){
for my $a4(@bases){
for my $a5(@bases){
for my $a6(@bases){
$data{"$a1$a2$a3$a4$a5$a6"} = 0;
}
}
}
}
}
}
my $nucs = <DATA>;
my $len = length($nucs);
for (my $i = 0; $i <= $len - 6; $i++) {
my $kmer = substr($nucs, $i, 6);
next if $kmer =~ tr/ACGT//c;
$data{$kmer}++; # populate hash with "seen" 6-mers
}
# print out sorted hash
foreach my $seq (sort { $data{$b} <=> $data{$a} } keys %data ){
print "$seq,$data{$seq}\n";
}
__DATA__
ATGCCCGTCGTAGTCATGCATGCATCGATCGATGCATGCTACGTGTTGT
显然会有一种 better/prettier 方法来计算字符串中字符的所有排列,而不是我所做的,但它确实有效。
正如 Borodin 所说,这主要打印出 "unseen" 字符串的变体。
最简单的方法是为每个模式构建一个出现次数的散列,然后打印所有可能模式的次数
此程序使用 glob
技巧生成由 A、T、C 和 G 组成的所有可能的六字符字符串列表
use strict;
use warnings 'all';
my @files = qw/ e_coli.fa /;
my %counts;
for my $file ( @files ) {
open my $fh, '<', $file or die qq{Unable to open "$file" for input: $!};
while ( <$fh> ) {
++$counts{} while /(?= ( [ATCG]{6} ) ) /gx;
}
}
for my $pattern ( glob '{A,T,C,G}' x 6 ) {
printf "%4d %s\n", $counts{$pattern} // 0, $pattern;
}
如果您有大量数据并且需要更快的速度,这里有一个 C 解决方案:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
void reader(FILE* in, unsigned long hist[4096]) {
for (unsigned long key=0, count=0;;) {
switch(getc(in)) {
case EOF: return;
case 'A': key <<= 2; break;
case 'C': key <<= 2; key += 1; break;
case 'G': key <<= 2; key += 2; break;
case 'T': key <<= 2; key += 3; break;
default: count=0; continue;
}
if (count == 5) ++hist[key & 0xFFF];
else ++count;
}
}
int putkey(FILE* out, unsigned long key) {
char s[6];
for (int j=6; j--; key >>= 2) s[j] = "ACGT"[key&3];
return fprintf(out, "%.6s", s);
}
void writer(FILE* out, unsigned long hist[4096]) {
for (unsigned long key = 0; key < 4096; ++key) {
fprintf(stdout, "%7lu ", hist[key]);
putkey(out, key);
putchar('\n');
}
}
int main(int argc, char** argv) {
FILE* in = stdin;
if (argc > 1) in = fopen(argv[1], "r");
if (!in) { perror(argv[1]); exit(1); }
unsigned long hist[4096] = {0};
reader(in, hist);
writer(stdout, hist);
return 0;
}
处理一个 31MB 的 fastq 样本(碰巧包括所有 4096 个可能的六字符序列)只花了不到半秒的时间; Perl 解决方案分别用了 12 秒 (fugu) 和 18 秒 (ikegami/borodin)。