perl 随机抽取文件中的一部分行

niuguohao

有一个基因组vcf文件，里面位点非常多（每个位点一行），现在想从这些位点中按照染色体上的位置，每500bp范围内随机抽取一个位点出来（比如下面在scaffold1中属于1-500位置的是前四行，然后从中抽一行，属于501-1000的是5，6，7三行，再从中随机抽出一行，以此类推），想了半天不知道该如何实现。
vcf文件格式如下(tab分隔)：

scaffold1       127     .       C       T       45.23   PASS AC=2 ;AF=0.011;=-……（后面还有很多故省略）
scaffold1       167     .       G       A       696.19  PASS AC=4 ;AF=0.02 2;0.……
scaffold1       260     .       A       G       1278.59 PASS AC=23 ;AF=0.1 11;5
scaffold1       462     .       T       C       4409.63 PASS AC=49 ;AF=0.2 33;……
scaffold1       636     .       A       G       7061.50 PASS AC=45 ;AF=0.2 12;=……
scaffold1       862     .       C       T       414.11  PASS AC=7 ;AF=0.03 3;;……
scaffold2       999     .       G       A       1135.55 PASS AC=4 ;AF=0.0 1949
scaffold2       128     .       C       T       45.23   PASS AC=2 ;AF=0.011 ;AN==-……
scaffold2       167     .       G       A       696.19  PASS  AC=4 ;AF=0. 022; AN=oeff=0.
scaffold2       699     .       A       G       7061.50 PASS AC=45;  AF=0.2 12;AN;DP=……
scaffold2       888     .       C       T       414.11  PASS AC=7 ;AF= 0.03 3;AN=;……
scaffold2       997     .       G       A       1135.55 PASS AC=4 ;AF=0.01 9;AN……
……
第一列是染色体名称 第二列是在该上染色体上的位置，后面是说明。
求助各位大神，应该如何操作呢？

523066680

暂时先这样。

1. my %hash;
   
2. for my $line ( <DATA> )
   
3. {
   
4.     next unless ( $line=~/(\w+)\s+(\w+)/ );
   
5.     #next unless ( $1 eq "scaffold1" );
   
6.     chomp $line;
   
7.     if ( $2 > 1 and $2 <= 500 ) {
   
8.         push @{$hash{ $1 }{'part1'}}, $line;
   
9.     }
   
10.     elsif ( $2 > 500 and $2 <= 1000 ) {
    
11.         push @{$hash{ $1 }{'part2'}}, $line;
    
12.     }
    
13. }
    
14. 
    
15. my $part1 = $hash{'scaffold1'}{'part1'};
    
16. my $part2 = $hash{'scaffold1'}{'part2'};
    
17. 
    
18. printf "%s\n", $part1->[rand($#{$part1}+1)];
    
19. printf "%s\n", $part2->[rand($#{$part2}+1)];
    
20. 
    
21. __DATA__
    
22. scaffold1       127     .       C       T       45.23   PASS AC=2 ;AF=0.011;=-...
    
23. scaffold1       167     .       G       A       696.19  PASS AC=4 ;AF=0.02 2;0....
    
24. scaffold1       260     .       A       G       1278.59 PASS AC=23 ;AF=0.1 11;5
    
25. scaffold1       462     .       T       C       4409.63 PASS AC=49 ;AF=0.2 33;...
    
26. scaffold1       636     .       A       G       7061.50 PASS AC=45 ;AF=0.2 12;=...
    
27. scaffold1       862     .       C       T       414.11  PASS AC=7 ;AF=0.03 3;;...
    
28. scaffold2       999     .       G       A       1135.55 PASS AC=4 ;AF=0.0 1949
    
29. scaffold2       128     .       C       T       45.23   PASS AC=2 ;AF=0.011 ;AN==-...
    
30. scaffold2       167     .       G       A       696.19  PASS  AC=4 ;AF=0. 022; AN=oeff=0.
    
31. scaffold2       699     .       A       G       7061.50 PASS AC=45;  AF=0.2 12;AN;DP=...
    
32. scaffold2       888     .       C       T       414.11  PASS AC=7 ;AF= 0.03 3;AN=;...
    
33. scaffold2       997     .       G       A       1135.55 PASS AC=4 ;AF=0.01 9;AN...

复制代码

niuguohao

回复 2# 523066680

多谢大神指点，但是我这个文件的第二列数字还有很多大于1000的，最大的可能有6，7位数，有什么办法改改么？

hztj2005

1. <div class="blockcode"><blockquote>
   

复制代码

1. use List::Util qw/min/;
   
2. use utf8;
   
3. 
   
4. my @datasource=(
   
5. ["scaffold1","127",],
   
6. ["scaffold1","167",],
   
7. ["scaffold1","260",],
   
8. ["scaffold1","462",],
   
9. ["scaffold1","636",],
   
10. ["scaffold1","862",],
    
11. ["scaffold1","999",],#原数据此行scaffold2似有误
    
12. ["scaffold2","128",],
    
13. ["scaffold2","167",],
    
14. ["scaffold2","699",],
    
15. ["scaffold2","888",],
    
16. ["scaffold2","997",],
    
17. );
    
18. 
    
19. my $sumline = @datasource;#行数
    
20. print $sumline."\n";;
    
21. 
    
22. my $endstr =$datasource[$#datasource][0];
    
23. if($endstr =~ /(\d+)/)
    
24. {
    
25.         $scasum = $1;#染色体数量
    
26.         #print $scasum."\n";
    
27. }
    
28. 
    
29. srand;#要先宣告srand函数，才能产生随机数的效果
    
30. my @onesca = ();
    
31. my $scastart =0;
    
32. my $scaend =0;
    
33. for(my $k=1;$k<$scasum+1;$k++)
    
34. {       
    
35.         my $maxvalue = -1;#该染色体最大位点
    
36.         for(my $j=$scastart;$j<$sumline;$j++)
    
37.         {          
    
38.           my $endstr2 =$datasource[$j][0];
    
39.           #print "-- $endstr2\n"; #
    
40.                 if($endstr2 =~ /(\d+)/)
    
41.                 {
    
42.                         $scaid = $1;#染色体序号
    
43.                         if($scaid == $k)
    
44.                         {                       
    
45.                          push @onesca,$datasource[$j][1];       
    
46.                          $maxvalue =$datasource[$j][1];
    
47.                          $scaend = $j+1;#记录处理下个染色体的开始行
    
48.                         }else
    
49.                         {
    
50.                           last;                               
    
51.                         }               
    
52.                 }               
    
53.         }       
    
54.        
    
55.         print "\n当前染色体最大位点：$maxvalue ||当前染色体开始行：$scastart  下个染色体的开始行：$scaend\n"; #       
    
56.         if(@onesca)
    
57.         {
    
58.                 for(my $baseid = 0; $baseid<$maxvalue+1; $baseid=$baseid+500)
    
59.           {       
    
60.                           pickline($baseid);
    
61.                 }               
    
62.         }       
    
63.         $scastart = $scaend;         
    
64.         print "\n =======第 $k 个染色体处理结束========\n"; #
    
65. }
    
66. 
    
67. sub pickline
    
68. {
    
69.             my $baseid=shift;            
    
70.             print ("当前染色体数组：@onesca\n");            
    
71.             my @part500 = grep $_<$baseid+500, @onesca;
    
72.             print ("500区间数组：@part500\n");          
    
73.           
    
74.                         my $myneed = int(rand(500)) + $baseid; #$myneed是一个0和500之间的整数
    
75.                         my @diff= map{abs($_ - $myneed)} @part500;
    
76.                         my $near = &min(@diff);  #差值最小，最靠近随机数
    
77.                         print "随机数->最小差值: $myneed -> $near\n";
    
78.                        
    
79.                         for(my $m=$scastart;$m<$scaend;$m++)
    
80.                         {
    
81.                                 my $posvalue = $datasource[$m][1];
    
82.                                 if(abs($posvalue-$myneed) == $near)
    
83.                                 {                                       
    
84.                                         print "需要的行：".$datasource[$m][0].",$posvalue\n\n";                                                                               
    
85.                                 }                                                       
    
86.                         }       
    
87.                        
    
88.                         #print "删除染色体数组中已处理的部分;       
    
89.                         while(@onesca and ($onesca[0] < ($baseid+500)))
    
90.                         {
    
91.                          shift @onesca; #删除
    
92.                         }
    
93.                        
    
94. 
    
95. }
    
96. 
    
97. exit;

复制代码

niuguohao

回复 4# hztj2005

谢谢您！我现在还没有完全看明白，稍后继续向您请教！非常感谢~！