perl - 随机选择一个区域并进行多次处理

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

我想从中随机选取一个包含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 

Note, we should not pick up the same region twice

Answer 1

我必须在这里假设一些事情，并保留一些限制

随机选择的区域不依赖于行；他们只是从文字

中挑选出来的

顺序无关紧要；应该通过文件

散布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;
}

请展开错误检查。例如，查看 this post和详细信息链接。

最简单的用法: driver_random.pl -n 4，但是您可以提供所有主程序的参数。

注意:被调用程序(上面的 random_regions.pl)必须是可执行的。

†从简单到功能更强大的一些: IPC::System::Simple， Capture::Tiny， IPC::Run3。 (然后在这里使用 IPC::Run。)另请参阅 String::ShellQuote，以准备命令而不引用问题，shell注入(inject)错误和其他问题。例如，请参阅 this post中组装的链接(示例)。