perl 隨機抽取文件中的一部分行

niuguohao

有一個基因組vcf文件，里面位點非常多（每個位點一行），現(xiàn)在想從這些位點中按照染色體上的位置，每500bp范圍內(nèi)隨機抽取一個位點出來（比如下面在scaffold1中屬于1-500位置的是前四行，然后從中抽一行，屬于501-1000的是5，6，7三行，再從中隨機抽出一行，以此類推），想了半天不知道該如何實現(xiàn)。
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

多謝大神指點，但是我這個文件的第二列數(shù)字還有很多大于1000的，最大的可能有6，7位數(shù)，有什么辦法改改么？

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",],#原數(shù)據(jù)此行scaffold2似有誤
    
12. ["scaffold2","128",],
    
13. ["scaffold2","167",],
    
14. ["scaffold2","699",],
    
15. ["scaffold2","888",],
    
16. ["scaffold2","997",],
    
17. );
    
18. 
    
19. my $sumline = @datasource;#行數(shù)
    
20. print $sumline."\n";;
    
21. 
    
22. my $endstr =$datasource[$#datasource][0];
    
23. if($endstr =~ /(\d+)/)
    
24. {
    
25.         $scasum = $1;#染色體數(shù)量
    
26.         #print $scasum."\n";
    
27. }
    
28. 
    
29. srand;#要先宣告srand函數(shù)，才能產(chǎn)生隨機數(shù)的效果
    
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 個染色體處理結(jié)束========\n"; #
    
65. }
    
66. 
    
67. sub pickline
    
68. {
    
69.             my $baseid=shift;            
    
70.             print ("當前染色體數(shù)組：@onesca\n");            
    
71.             my @part500 = grep $_<$baseid+500, @onesca;
    
72.             print ("500區(qū)間數(shù)組：@part500\n");          
    
73.           
    
74.                         my $myneed = int(rand(500)) + $baseid; #$myneed是一個0和500之間的整數(shù)
    
75.                         my @diff= map{abs($_ - $myneed)} @part500;
    
76.                         my $near = &min(@diff);  #差值最小，最靠近隨機數(shù)
    
77.                         print "隨機數(shù)->最小差值: $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 "刪除染色體數(shù)組中已處理的部分;       
    
89.                         while(@onesca and ($onesca[0] < ($baseid+500)))
    
90.                         {
    
91.                          shift @onesca; #刪除
    
92.                         }
    
93.                        
    
94. 
    
95. }
    
96. 
    
97. exit;

復制代碼

niuguohao

回復 4# hztj2005

謝謝您！我現(xiàn)在還沒有完全看明白，稍后繼續(xù)向您請教！非常感謝~！