我有这样的数据
>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
我想从中随机选取一个包含10个字母的区域，然后计算F的数量，我想这样做一定次数，例如1000次甚至更多次
例如，我随机选择
LVPSLTRYLT 0
然后
ITNLRSFIHK 1
然后再随机去接连续10个字母
AHSRIRKERP 0
这一直持续到达到要求的运行次数为止。我想将所有随机选择的值与它们的值一起存储，因为那样我想计算出看到F的次数
所以我做以下
# first I remove the header
grep -v ">" data.txt > out.txt
然后随机获得一个我尝试使用的10个字母的区域shuf，但没有成功，
shuf -n1000 data.txt
然后我尝试使用awk但也没有成功
awk 'BEGIN {srand()} !/^$/ { if (rand() == 10) print $0}'
然后计算F的数量并将其保存在文件中
grep -i -e [F] |wc -l
请注意，我们不应两次拾取相同的区域
我必须在这里假设一些事情，并保留一些限制
选择的随机区域不取决于线；他们只是从文本中挑选出来的
顺序无关紧要；文件中应该只分布N个区域
文件的大小可以达到千兆字节，因此无法先完整读取（会容易得多！）
有未处理（边缘或不太可能）的情况，在代码后进行讨论
首先建立一个排序的随机数列表；这些是文件中区域开始的位置。然后，在读取每一行时，计算其在文件中的字符范围，并检查我们的数字是否在其中。如果有的话，它们会标记每个随机区域的开头：从这些字符开始，选择所需长度的子字符串。检查子字符串是否适合该行。
use warnings;
use strict;
use feature 'say';
use Getopt::Long;
use List::MoreUtils qw(uniq);
my ($region_len, $num_regions) = (10, 10);
my $count_freq_for = 'F';
#srand(10);
GetOptions(
'num-regions|n=i' => \$num_regions,
'region-len|l=i' => \$region_len,
'char|c=s' => \$count_freq_for,
) or usage();
my $file = shift || usage();
# List of (up to) $num_regions random numbers, spanning the file size
# However, we skip all '>sp' lines so take more numbers (estimate)
open my $fh, '<', $file or die "Can't open $file: $!";
$num_regions += int $num_regions * fraction_skipped($fh);
my @rand = uniq sort { $a <=> $b }
map { int(rand (-s $file)-$region_len) } 1..$num_regions;
say "Starting positions for regions: @rand";
my ($nchars_prev, $nchars, $chars_left) = (0, 0, 0);
my $region;
while (my $line = <$fh>) {
chomp $line;
# Total number of characters so far, up to this line and with this line
$nchars_prev = $nchars;
$nchars += length $line;
next if $line =~ /^\s*>sp/;
# Complete the region if there wasn't enough chars on the previous line
if ($chars_left > 0) {
$region .= substr $line, 0, $chars_left;
my $cnt = () = $region =~ /$count_freq_for/g;
say "$region $cnt";
$chars_left = -1;
};
# Random positions that happen to be on this line
my @pos = grep { $_ > $nchars_prev and $_ < $nchars } @rand;
# say "\tPositions on ($nchars_prev -- $nchars) line: @pos" if @pos;
for (@pos) {
my $pos_in_line = $_ - $nchars_prev;
$region = substr $line, $pos_in_line, $region_len;
# Don't print if there aren't enough chars left on this line
last if ( $chars_left =
($region_len - (length($line) - $pos_in_line)) ) > 0;
my $cnt = () = $region =~ /$count_freq_for/g;
say "$region $cnt";
}
}
sub fraction_skipped {
my ($fh) = @_;
my ($skip_len, $data_len);
my $curr_pos = tell $fh;
seek $fh, 0, 0 if $curr_pos != 0;
while (<$fh>) {
chomp;
if (/^\s*>sp/) { $skip_len += length }
else { $data_len += length }
}
seek $fh, $curr_pos, 0; # leave it as we found it
return $skip_len / ($skip_len+$data_len);
}
sub usage {
say STDERR "Usage: $0 [options] file", "\n\toptions: ...";
exit;
}
取消注释该srand行，以便始终进行相同的运行以进行测试。注意事项如下。
一些极端情况
如果10个长的窗口从其随机位置开始不适合该行，则在下一行中完成它-但该行上的任何（可能）其他随机位置都将被忽略。因此，如果我们的随机列表有1120和1122，而一条线在1125结束，则跳过从1122开始的窗口。不太可能，没有任何结果（除了减少一个区域）。
当一个不完整的区域中的下一行（第一填充if在while循环中），则可能的是该行是比其余所需的字符短（$chars_left）。这是极不可能的，因此需要在此处进行额外的检查。
随机数被删除。这使顺序变得歪斜，但是在此刻不重要的是什么。我们可能留下的数量少于要求的数量，但数量很少
处理有关随机性的问题
这里的“随机性”是很基本的，似乎很合适。我们还需要考虑以下内容。
在文件大小int(rand -s $file)（减去区域大小）的间隔内绘制随机数。但是>sp会跳过行，并且不会使用可能落在这些行内的任何数字，因此最终得出的区域可能少于绘制的数字。这些行较短，因此在其上具有数字的机会较小，因此丢失了很多数字，但是在某些运行中，我看到甚至跳过了十分之三的数字，最终得到了所需样本大小的70％。
如果这很麻烦，则有一些方法可以解决。为了不进一步分散发行版，它们都应涉及预处理文件。
上面的代码对文件进行了初始运行，以计算将被跳过的字符比例。然后将其用于增加绘制的随机点的数量。这当然是“平均”的措施，但仍应产生接近足够大文件所需区域的数量。
更详细的度量将需要查看（更大）分布的哪些随机点将丢失到跳过的行中，然后重新采样以解决该问题。这可能仍然会与发行版混淆，这在这里可能不是问题，但更重要的是根本不需要。
在所有这些中，您两次读取了大文件。额外的处理时间应仅以秒为单位，但是如果不可接受，则将功能更改为fraction_skipped仅读取文件的10-20％。对于大文件，这仍应提供合理的估计。
注意特定的测试用例
使用srand(10)（在开头附近注释掉的行），我们得到随机数，使得在一行上该区域在行尾之前开始8个字符！因此，这种情况确实会测试代码以完成下一行的区域。
一个简单的驱动程序可以运行上述给定次数，以进行统计。
仅使用内置工具（system，qx）进行操作会很麻烦而且很麻烦；实际上，一个需要模块。所以我在这里使用IPC :: Run。还有很多其他选择。†
调整并添加代码以根据需要进行统计；输出在文件中。
use warnings;
use strict;
use feature 'say';
use Getopt::Long;
use IPC::Run qw(run);
my $outdir = 'rr_output'; # pick a directory name
mkdir $outdir if not -d $outdir;
my $prog = 'random_regions.pl'; # your name for the program
my $input = 'data_file.txt'; # your name for input file
my $ch = 'F';
my ($runs, $regions, $len) = (10, 10, 10);
GetOptions(
'runs|n=i' => \$runs,
'regions=i' => \$regions,
'length=i' => \$len,
'char=s' => \$ch,
'input=s' => \$input
) or usage();
my @cmd = ( $prog, $input,
'--num-regions', $regions,
'--region-len', $len,
'--char', $ch
);
say "Run: @cmd, $runs times.";
for my $n (1..$runs) {
my $outfile = "$outdir/regions_r$n.txt";
say "Run #$n, output in: $outdir/$outfile";
run \@cmd, '>', $outfile or die "Error with @cmd: $!";
}
sub usage {
say STDERR "Usage: $0 [options]", "\n\toptions: ...";
exit;
}
请展开错误检查。例如，请参阅这篇文章和详细信息链接。
最简单的用法：driver_random.pl -n 4，但是您可以提供所有主程序的参数。
注意：被调用的程序（random_regions.pl上面）必须是可执行的。
† 从简单到更强大：IPC :: System :: Simple，Capture :: Tiny，IPC :: Run3。（然后IPC::Run在这里使用。）另请参见String :: ShellQuote，以准备命令而不引用问题，shell注入错误以及其他问题。见链接（例子）在组装此篇，例如。