亚洲av成人无遮挡网站在线观看,少妇性bbb搡bbb爽爽爽,亚洲av日韩精品久久久久久,兔费看少妇性l交大片免费,无码少妇一区二区三区

Chinaunix

標(biāo)題: perl 隨機抽取文件中的一部分行 [打印本頁]
作者: niuguohao    時間: 2018-12-22 17:25
標(biāo)題: perl 隨機抽取文件中的一部分行
本帖最后由 niuguohao 于 2018-12-22 17:28 編輯

有一個基因組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……
……
第一列是染色體名稱 第二列是在該上染色體上的位置，后面是說明。
求助各位大神，應(yīng)該如何操作呢？
作者: 523066680    時間: 2018-12-23 16:14
暫時先這樣。

  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. }

  15. my $part1 = $hash{'scaffold1'}{'part1'};
  16. my $part2 = $hash{'scaffold1'}{'part2'};

  18. printf "%s\n", $part1->[rand($#{$part1}+1)];
  19. printf "%s\n", $part2->[rand($#{$part2}+1)];

  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...
復(fù)制代碼
作者: niuguohao    時間: 2018-12-23 19:04
回復(fù) 2# 523066680

多謝大神指點，但是我這個文件的第二列數(shù)字還有很多大于1000的，最大的可能有6，7位數(shù)，有什么辦法改改么？
作者: hztj2005    時間: 2018-12-23 20:14
本帖最后由 hztj2005 于 2018-12-28 09:39 編輯
  1. <div class="blockcode"><blockquote>
復(fù)制代碼

  1. use List::Util qw/min/;
  2. use utf8;

  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. );

  19. my $sumline = @datasource;#行數(shù)
  20. print $sumline."\n";;

  22. my $endstr =$datasource[$#datasource][0];
  23. if($endstr =~ /(\d+)/)
  24. {
  25.         $scasum = $1;#染色體數(shù)量
  26.         #print $scasum."\n";
  27. }

  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當(dāng)前染色體最大位點：$maxvalue ||當(dāng)前染色體開始行：$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. }

  67. sub pickline
  68. {
  69.             my $baseid=shift;            
  70.             print ("當(dāng)前染色體數(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.                        

  95. }

  97. exit;
復(fù)制代碼


作者: niuguohao    時間: 2018-12-25 21:55
回復(fù) 4# hztj2005

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




歡迎光臨 Chinaunix (http://www.72891.cn/) Powered by Discuz! X3.2