请教一个关于范围操作符（..）的问题

飞越地平线

回复 10# zhlong8

@bed里的文件（也就是目录B下的文件）有20来个，每个接近100M。之前我怕我的电脑内存不够就没有这样处理了。
用12G的电脑是不是就快多了？
直接操作数组是不是指这样？

1. #用@bedfile直接全部存入split后的bed文件内容，然后类似下面这样统计？
   
2. foreach my $line(@bedfile){
   
3.         my $sNum=$line->[1];
   
4.         # stop condition
   
5.         last if($sNum > $sMax);
   
6.         # start condition
   
7.         if($sNum >= $sMin){
   
8.                 $num++;
   
9.                 next;
   
10.         }
    
11. }

复制代码

zhlong8

回复 11# 飞越地平线


    你这么写12G是远远不够的，只保存 $line->[1] 即数据中的数字部分才勉强够，而用C只用1G多点就够了，这就是我前面说的 Perl 的限制。

jason680

回复 8# 飞越地平线

Sorry, it's hard to know what you want
Note: you didn't give us the example for result


and hard to know what you said ...
4楼
目录B有多个文件（文件名为chr1, chr2, chr3)，比如chr1文件内的内容像这样：
...

8楼
目录B下的文件名称是像bed_chr_1.bed, bed_chr_2.bed这样的   
...

xiumu2280

文件B的bed文件看样子应该是分染色体的。你循环A文件的时候多次重复打开文件句柄。
所以你要把1号染色体的位点放到一起，然后循环，只打开一次文件句柄。
文件B是SNP信息。
你在计算位点数量的时候，可以把位点存到数组，计算数量，要简单一些。
另外你用的map实在太多了。
threads这是个多线程的模块。

飞越地平线

回复 13# jason680

不好意思。4楼我说目录B里的文件（chr1, chr2, chr3）其实就是 8楼说的（bed_chr_1.bed, bed_chr_2.bed ), 为了简单化才这么说的。

如果你懂生物的话应该大概知道我想做的是什么，其实我就是想计算基因的SNP（单核苷酸多态）数。

目录A里有多个文件，每个文件保存着一类的基因（遍布各染色体），每行都是一个基因信息。比如（chr1 3 7 aa）指的就是（基因名称为aa， 在1号染色体（chr就是chromosome）的起始3，终止7的位置上）；

目录B也有多个文件（bed_chr_1.bed, bed_chr_2.bed, bed_chr_3.bed )。比如文件bed_chr_1.bed的内容就是保存着1号染色体所有位置的SNP信息。比如一行里（chr1 7）就表示在1号染色体的7号位有一个SNP。

目的是统计目录A里所有的文件里的基因的SNP信息。目录A里有多少文件，就同样地生成多少文件。这些文件的结果里每行是（一个名称，一个SNP数）像是这样：
aa 10
bb 3
cc 5
dd 89

大概这样解释了。其中有些地方（例如bed文件的格式）不好讲，代码里统计SNP数的地方就不用太细看了。

飞越地平线

回复 12# zhlong8

嗯，我再试试。

   

飞越地平线

回复 14# xiumu2280

谢谢关注！
如果把1号染色体的位点放到一起，然后循环，每次处理一行的时候，不是都需要重复的打开bed_ch_1.bed文件吗？
“你在计算位点数量的时候，可以把位点存到数组，计算数量，要简单一些。”这句话不懂。直接用$num记录？
我看了下，14行和15行的map可以并在一起的。 另外想问下，map和foreach效率怎么样？

想问一下，从哪里看到threads的模块说明？
   

rubyish

是不是指这样？

1. #!/usr/bin/perl
   
2. use 5.018;
   
3. @ARGV = glob "/tmp/A/*";
   
4. my %A;
   
5. while (<>) {
   
6.     my ( $C, $N, $X, $G ) = split;
   
7.     push @{ $A{$C} }, { min => $N, max => $X, G => $G, F => $ARGV };
   
8. }
   
9. 
   
10. while ( my ( $k, $v ) = each %A ) {
    
11.     my ($N) = $k =~ /(\d+)$/;
    
12.     open my $B, "/tmp/B/bed_chr_$N.bed" or next;
    
13.     my @CHR = map { (split)[1] } <$B>;
    
14.     for my $S (@$v) {
    
15.         open my $NSP, '>>', $S->{F} . '.nsp';
    
16.         say $NSP $S->{G}, "\t", SUM( \@CHR, $S->{min}, $S->{max} );
    
17.         close $NSP;
    
18.     }
    
19. }

复制代码

sub SUM:

1. sub SUM {
   
2.     my ( $A, $N, $X, $S ) = @_;
   
3.     return 0 if $N > $A->[ $#{$A} ] or $X < $A->[0];
   
4.     for (@$A) {
   
5.       last if $_ > $X;
   
6.       next if $_ < $N;
   
7.       $S++;
   
8.     }
   
9.     $S;
   
10. }

复制代码

OR sub SUM:

1. sub SUM {
   
2.     my ( $A, $N, $X, $C ) = @_;
   
3.     return 0 if $N > $A->[ $#{$A} ] or $X < $A->[0];
   
4.     my $bs = sub {
   
5.         my ( $W, $l, $r, $m ) = ( shift, shift, $#{$A} );
   
6.         return $r if $W >= $A->[$r];
   
7.         return 0  if $W <= $A->[$l];
   
8.         goto R    if $C;
   
9.         while ( $l <= $r ) {
   
10.             $m = int 0.5 * ( $l + $r );
    
11.             $W > $A->[ $m - 1 ]
    
12.               ? $W <= $A->[$m] ? return $m : ( $l = $m + 1 )
    
13.               : ( $r = $m - 1 );
    
14.         }
    
15.       R: while ( $l <= $r ) {
    
16.             $m = int 0.5 * ( $l + $r );
    
17.             $W < $A->[ $m + 1 ]
    
18.               ? $W >= $A->[$m] ? return $m : ( $r = $m - 1 )
    
19.               : ( $l = $m + 1 );
    
20.         }
    
21.     };
    
22.     my $L = $bs->( $N, 0 );
    
23.     $bs->( $X, $L, $C++ ) - $L + 1;
    
24. }

复制代码

xiumu2280

1.你可以看下  programing perl中的 哈希的数组 这一章。文件A用chr1和基因名做成key,然后把后面碱基范围作为VALUE存到一个数组。可以把A里的所有文件cat到一起。
2.因为B类文件，应该是按染色体分类的吧？ 把每个文件中的数据位点存成hash的key，这样打开一次就可以了。然后循环上面hash中的value。把存在的结果push到一个数组里。算数量的时候就$num=@data。
3.map和foreach效率我还真不知道···map挺好，用了之后括号就没那么乱了。我的意思是你map很长，；的时候可以换一行。
4.去CPAN上看模块啊··
回复 17# 飞越地平线


   

xiumu2280

1. #!/usr/bin/perl
   
2. use strict;
   
3. use warnings;
   
4. ###默认你把A目录下文件cat到一起了
   
5. open (IN,"A");
   
6. my %hash_A;
   
7. while (<IN>) {
   
8.         chomp;
   
9.         my ($chr,$pos_start,$pos_end,$gene) = split;
   
10.         $hash_A{$chr}{$gene}=[($pos_start..$pos_end)];
    
11. }
    
12. close IN;
    
13. my @bed_B = glob "*bed";
    
14. foreach  my $bedB (@bed_B) {
    
15.                 my @bed_name = split /\/|_|./,$bedB;
    
16.                 my $chr_bed = $bed_name[-3].$bed_name[-2];
    
17.         my %hash_bed;
    
18.                 open (IN1,"$bedB");
    
19.         %hash_bed = map {chomp;my @data=split;$data[1],1;} <IN1>;#如果后面有碱基信息可以存成value.
    
20.         close IN1;
    
21.         foreach my $gene_A (keys %{$hash_A{$chr_bed}}) {
    
22.                 my @snp_num;
    
23.                 my $num;
    
24.                 open (OUT,">$chr_A-$gene_A");
    
25.                 foreach my $pos_A (@{$hash_A{$chr_A}{$gene_A}}) {
    
26.                                 if ($hash_bed{$pos_A}) {
    
27.                                         push @snp_num,$hash_bed{$pos_A};
    
28.                                 }
    
29.                 }
    
30.                 $num = @snp_num;
    
31.                 print OUT "$gene_A\t$num";
    
32.                 close OUT;
    
33.                 }
    
34. }
    

复制代码

不知道你是不是这个意思···