免费注册 查看新帖 |

Chinaunix

  平台 论坛 博客 文库
最近访问板块 发新帖
查看: 19298 | 回复: 3

perl 随机抽取文件中的一部分行 [复制链接]

论坛徽章:
0
发表于 2018-12-22 17:25 |显示全部楼层
本帖最后由 niuguohao 于 2018-12-22 17:28 编辑

有一个基因组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……
……
第一列是染色体名称 第二列是在该上染色体上的位置,后面是说明。
求助各位大神,应该如何操作呢?

论坛徽章:
12
子鼠
日期:2014-10-11 16:46:482016科比退役纪念章
日期:2018-03-16 10:24:0515-16赛季CBA联赛之山东
日期:2017-11-10 14:32:142016科比退役纪念章
日期:2017-09-02 15:42:4715-16赛季CBA联赛之佛山
日期:2017-08-28 17:11:5515-16赛季CBA联赛之浙江
日期:2017-08-24 16:55:1715-16赛季CBA联赛之青岛
日期:2017-08-17 19:55:2415-16赛季CBA联赛之天津
日期:2017-06-29 10:34:4315-16赛季CBA联赛之四川
日期:2017-05-16 16:38:55黑曼巴
日期:2016-07-19 15:03:112015亚冠之萨济拖拉机
日期:2015-05-22 11:38:5315-16赛季CBA联赛之北京
日期:2019-08-13 17:30:53
发表于 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. }

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

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

  18. __DATA__
  19. scaffold1       127     .       C       T       45.23   PASS AC=2 ;AF=0.011;=-...
  20. scaffold1       167     .       G       A       696.19  PASS AC=4 ;AF=0.02 2;0....
  21. scaffold1       260     .       A       G       1278.59 PASS AC=23 ;AF=0.1 11;5
  22. scaffold1       462     .       T       C       4409.63 PASS AC=49 ;AF=0.2 33;...
  23. scaffold1       636     .       A       G       7061.50 PASS AC=45 ;AF=0.2 12;=...
  24. scaffold1       862     .       C       T       414.11  PASS AC=7 ;AF=0.03 3;;...
  25. scaffold2       999     .       G       A       1135.55 PASS AC=4 ;AF=0.0 1949
  26. scaffold2       128     .       C       T       45.23   PASS AC=2 ;AF=0.011 ;AN==-...
  27. scaffold2       167     .       G       A       696.19  PASS  AC=4 ;AF=0. 022; AN=oeff=0.
  28. scaffold2       699     .       A       G       7061.50 PASS AC=45;  AF=0.2 12;AN;DP=...
  29. scaffold2       888     .       C       T       414.11  PASS AC=7 ;AF= 0.03 3;AN=;...
  30. scaffold2       997     .       G       A       1135.55 PASS AC=4 ;AF=0.01 9;AN...
复制代码

论坛徽章:
0
发表于 2018-12-23 19:04 |显示全部楼层
回复 2# 523066680

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

论坛徽章:
0
发表于 2018-12-23 20:14 |显示全部楼层
本帖最后由 hztj2005 于 2018-12-28 09:39 编辑
  1. <div class="blockcode"><blockquote>
复制代码

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

  3. my @datasource=(
  4. ["scaffold1","127",],
  5. ["scaffold1","167",],
  6. ["scaffold1","260",],
  7. ["scaffold1","462",],
  8. ["scaffold1","636",],
  9. ["scaffold1","862",],
  10. ["scaffold1","999",],#原数据此行scaffold2似有误
  11. ["scaffold2","128",],
  12. ["scaffold2","167",],
  13. ["scaffold2","699",],
  14. ["scaffold2","888",],
  15. ["scaffold2","997",],
  16. );

  17. my $sumline = @datasource;#行数
  18. print $sumline."\n";;

  19. my $endstr =$datasource[$#datasource][0];
  20. if($endstr =~ /(\d+)/)
  21. {
  22.         $scasum = $1;#染色体数量
  23.         #print $scasum."\n";
  24. }

  25. srand;#要先宣告srand函数,才能产生随机数的效果
  26. my @onesca = ();
  27. my $scastart =0;
  28. my $scaend =0;
  29. for(my $k=1;$k<$scasum+1;$k++)
  30. {       
  31.         my $maxvalue = -1;#该染色体最大位点
  32.         for(my $j=$scastart;$j<$sumline;$j++)
  33.         {          
  34.           my $endstr2 =$datasource[$j][0];
  35.           #print "-- $endstr2\n"; #
  36.                 if($endstr2 =~ /(\d+)/)
  37.                 {
  38.                         $scaid = $1;#染色体序号
  39.                         if($scaid == $k)
  40.                         {                       
  41.                          push @onesca,$datasource[$j][1];       
  42.                          $maxvalue =$datasource[$j][1];
  43.                          $scaend = $j+1;#记录处理下个染色体的开始行
  44.                         }else
  45.                         {
  46.                           last;                               
  47.                         }               
  48.                 }               
  49.         }       
  50.        
  51.         print "\n当前染色体最大位点:$maxvalue ||当前染色体开始行:$scastart  下个染色体的开始行:$scaend\n"; #       
  52.         if(@onesca)
  53.         {
  54.                 for(my $baseid = 0; $baseid<$maxvalue+1; $baseid=$baseid+500)
  55.           {       
  56.                           pickline($baseid);
  57.                 }               
  58.         }       
  59.         $scastart = $scaend;         
  60.         print "\n =======第 $k 个染色体处理结束========\n"; #
  61. }

  62. sub pickline
  63. {
  64.             my $baseid=shift;            
  65.             print ("当前染色体数组:@onesca\n");            
  66.             my @part500 = grep $_<$baseid+500, @onesca;
  67.             print ("500区间数组:@part500\n");          
  68.           
  69.                         my $myneed = int(rand(500)) + $baseid; #$myneed是一个0和500之间的整数
  70.                         my @diff= map{abs($_ - $myneed)} @part500;
  71.                         my $near = &min(@diff);  #差值最小,最靠近随机数
  72.                         print "随机数->最小差值: $myneed -> $near\n";
  73.                        
  74.                         for(my $m=$scastart;$m<$scaend;$m++)
  75.                         {
  76.                                 my $posvalue = $datasource[$m][1];
  77.                                 if(abs($posvalue-$myneed) == $near)
  78.                                 {                                       
  79.                                         print "需要的行:".$datasource[$m][0].",$posvalue\n\n";                                                                               
  80.                                 }                                                       
  81.                         }       
  82.                        
  83.                         #print "删除染色体数组中已处理的部分;       
  84.                         while(@onesca and ($onesca[0] < ($baseid+500)))
  85.                         {
  86.                          shift @onesca; #删除
  87.                         }
  88.                        

  89. }

  90. exit;
复制代码


论坛徽章:
0
发表于 2018-12-25 21:55 |显示全部楼层
回复 4# hztj2005

谢谢您!我现在还没有完全看明白,稍后继续向您请教!非常感谢~!
您需要登录后才可以回帖 登录 | 注册

本版积分规则 发表回复

SACC2019中国系统架构师大会

【数字转型 架构演进】SACC2019中国系统架构师大会
2019年10月31日~11月2日第11届中国系统架构师大会(SACC2019)将在北京隆重召开。四大主线并行的演讲模式,1个主会场、20个技术专场、超千人参与的会议规模,100+来自互联网、金融、制造业、电商等领域的嘉宾阵容,将为广大参会者提供一场最具价值的技术交流盛会。




----------------------------------------

大会官网>>
  

北京盛拓优讯信息技术有限公司. 版权所有 16024965号-6 北京市公安局海淀分局网监中心备案编号:11010802020122
中国互联网协会会员  联系我们:huangweiwei@it168.com
感谢所有关心和支持过ChinaUnix的朋友们 转载本站内容请注明原作者名及出处

清除 Cookies - ChinaUnix - Archiver - WAP - TOP