免费注册 查看新帖 |

Chinaunix

  平台 论坛 博客 文库
最近访问板块 发新帖
查看: 3866 | 回复: 9
打印 上一主题 下一主题

求助:对数组以5000为步长采样,统计问题! [复制链接]

论坛徽章:
0
跳转到指定楼层
1 [收藏(0)] [报告]
发表于 2011-10-28 11:19 |只看该作者 |倒序浏览
本帖最后由 linghu002127 于 2011-10-28 14:30 编辑

数据:手上有类似数据几百万行,简单贴几行如下(各行之间tab键分开):
H2AK5ac        1000
H2AK5ac        5008
H2AK5ac        10001
H2BK5me        5008
H2BK5me        10004
H4K12ac        10007
H4K12ac        15005
H4R3me2        5010
H4R3me2        10017
此数据为2列,第一列有4个文件名依次为:H2AK5ac,H2BK5me,H4K12ac,H4R3me2。第二列为第一列对应的数值(相同文件名对应的数字都是由小到大),如H2AK5ac对应三个数依次为:1000,5008,10001。H2BK5me对应两个数依次为:5008,10004。(假设$min,$max为整个第二列数据的最小最大值,本例中$min=1000,$max=15005)
在相同文件名下,我想每一次用尺度5000去卡这些数,如果这些数在我5000这个窗内,则结果输出1,如果不在就输出0。相当于数据不动,只是窗在动。当窗滑动到$max范围内时,就不滑动了。比如对H2AK5ac这个文件名来说有三数1000,5008,10001 (如下图红点),滑动5000窗,依次应该输出的数据为:1  1  1  0 ;对H2BK5me有两个数:5008,10004(图中未标出),依次应该输出的数据为:0  1  1  0 ;对于不同的文件名,每次的处理都一样,窗都是从[0~5000]开始,到$max范围时结束(本例中滑到[15000-20000]窗或者最多[20000-25000]窗应该停止)。

即最终得到如下格式的数据(各行之间均tab键隔开):

论坛徽章:
0
2 [报告]
发表于 2011-10-28 11:25 |只看该作者
对最后结果解释如下:
1)首先第一行为表头,由第一列不重复的文件名构成,前面在加上字符“ID(tab键)NAME”,如下:


2)从第2行开始,为真正需要的数据:
前两列分别为ID和NAME:ID值从1开始递增,相当于计数,处理结果有多少行,ID就增加到多少;NAME值从chrY-5000开始,向下依次递增5000。
后四列依次为H2AK5ac,H2BK5me,H4K12ac,H4R3me对应数据处理的结果。如第三列为H2AK5ac的处理结果,第4列为H2BK5me处理结果等等。也就是说,这样能保证第一行数据都在[0~5000]窗口,第二行数据在[5000~10000]窗口,依次类推。

以上数据我折腾快5天了,生物背景数据结构不太强,求各位帮忙看看~~

论坛徽章:
0
3 [报告]
发表于 2011-10-28 11:41 |只看该作者
本帖最后由 iamlimeng 于 2011-10-28 11:48 编辑

这应该是一个简单的问题,容易搞定。LZ应该上传一小段数据。

但我仍然没看懂LZ的要求!

论坛徽章:
0
4 [报告]
发表于 2011-10-28 13:14 |只看该作者
本帖最后由 linghu002127 于 2011-10-28 14:33 编辑

回复 3# iamlimeng


关于数据:我给的这个示例数据(chrY.txt)把数改小了,是方便看出来。示例结果为chrY_result.txt。

原始数据(chrY,用写字板,notepad之类的能打开)很大,而且初始值也很大。原始数据第一列有39个不同的变量名(分别从CD4-H2AK5ac到H4R3me2,表示有39种变化),对于每一个变量名对应的后边的数(第二列,表示变化出现的位置),都是从小到大排列的。

我都一并上传。


原始数据是关于人类染色体Y的数据。Y染色体(记为chrY)长度为500000000(可以想象为一把尺子从1开始编号一直到500000000)。现在在这条染色体上发生了39种变化(也就是第二列记录的数字,如57255836,代表在这把尺子的这个位置发生了变化)。每一种变化互不影响,都在这条染色体上发生。

我主要的意思是:(什么表头呀,ID,NAME那些先不管),对于每一种变化,比如第一种变化(用变量名CD4-H2AK5ac记),在原始数据中它后边的那些所有数字代表这种变化在相应的点上发生了(想象这把尺子上相应位置上已有标记)。现在我以5000长度的窗口去找到这些标记,窗口从0~5000开始(其实就是设置一个变量让它每次循环加5000),如果发现某一个变化刻度在我的这个窗口内,就输出1,窗口内没有变化就输出0,窗口每滑动一次都有输出(0或1),直到滑动到最大一个变化刻度为止(没有必要滑动到染色体的最大长度500000000,这样的话会有很多0的,这也是为什么要先找到整个第二列最大值的原因)。

data.zip

7.11 KB, 下载次数: 37

论坛徽章:
0
5 [报告]
发表于 2011-10-28 13:54 |只看该作者
这是我不成文的代码~
感到逻辑有点乱~·
统计和比较那一块不会写,输出格式也不会调成按列输出。

请各位指点指点~~
  1. #!/usr/bin/perl -w


  2. use strict;

  3. my @temp;
  4. my %hash;

  5. while(<DATA>)
  6. {
  7.         push @{$hash{(split)[0]}},(split)[1];
  8.         push @temp, (split)[1];
  9. }
  10. my $minvalue=&min(@temp);
  11. my $maxvalue=&max(@temp);

  12. open(FH,">>chrY_result.txt");

  13. #<====此段输出表头信息===>#
  14. my @item;
  15. foreach my $cur_key (sort keys %hash) {
  16.                 push @item,$cur_key."\t";
  17. }
  18. @item=reverse @item;
  19. print FH "ID\tNAME\t@item","\n";

  20. #<====此段处理数据===>#
  21. foreach my $mykey1 (keys %hash){
  22.     for (my $i=5000;$i<=$maxvalue;$i+=5000){
  23.                                
  24.                 foreach my $element (@{$hash{$mykey1}}){
  25.                         #这里的循环和输出格式不会写~
  26.                         #也许是我的程序本就出了问题,求帮忙!!
  27.                 }
  28.         }
  29. }

  30. close FH;

  31. sub max{
  32.         my($max1)=shift @_;
  33.         foreach(@_){
  34.         if($_>$max1) {$max1=$_;}
  35.         }
  36.         $max1;
  37. }

  38. sub min{
  39.         my($min1)=shift @_;
  40.         foreach(@_){
  41.         if($_<$min1) {$min1=$_;}
  42.         }
  43.         $min1;
  44. }

  45. __DATA__
  46. H2AK5ac        1000
  47. H2AK5ac        5008
  48. H2AK5ac        10001
  49. H2BK5me        5008
  50. H2BK5me        10004
  51. H4K12ac        10007
  52. H4K12ac        15005
  53. H4R3me2        5010
  54. H4R3me2        10017
复制代码

论坛徽章:
0
6 [报告]
发表于 2011-10-28 14:52 |只看该作者
--看懂了,但输出 10010101 (随便写的,不用理会)--这不会是你最终的要求吧...

把你最终的要求讲清楚,以免同样的东西,大家替你重复地写....

论坛徽章:
0
7 [报告]
发表于 2011-10-28 15:45 |只看该作者
本帖最后由 linghu002127 于 2011-10-28 15:53 编辑

回复 6# Kitaisky


    是的,主要部分就是一些0,1之类的代码。只是另外还有表头呀,ID之类,这主要是一些格式问题~
我在附件data中给出了简单示例chrY.txt的输出结果(文件chrY_result.txt)
贴出来是这样(tab键隔开的,看起来对的不太整齐,不过用notepade就对齐了):
  1. ID        NAME        H2AK5ac        H2BK5me        H4K12ac        H4R3me2
  2. 1        chrY-5000        1        0        0        0
  3. 2        chrY-10000        1        1        0        1
  4. 3        chrY-15000        1        1        1        1
  5. 4        chrY-20000        0        0        1        0
复制代码
不过对于实际数据(chrY).
第一列可以看出有39个键值(从H2AK5ac到H4R3me2,表示有39种变化),对于每一个键值(如H2AK5ac)对应的后边的第二列数(从10629344 到57403474),表示变化出现的位置。所有键值对应的这些数都是从小到大排列好的。

希望设置一个循环:初始值=$min,步长=5000,终止值=$max    (假设$min,$max为整个第二列所有数据的最小最大值)
对39个键值下边的数据都进行同样的处理。
如果发现某一个变化刻度在我的这个窗口内,就输出1,窗口内没有变化就输出0,窗口每滑动一次都有输出(0或1),直到滑动到最大一个变化刻度为止。输出按列输出,比如上边H2AK5ac的输出为它对应的那一列1 1 1 0。 而第二个键值H2BK5me的输出为0 1 1 0

论坛徽章:
0
8 [报告]
发表于 2011-10-28 16:01 |只看该作者
本帖最后由 linghu002127 于 2011-10-28 17:01 编辑

上述循环:初始值=$min,步长=5000,终止值=$max   
初始值不是从0开始,而从$min开始,是因为如果从0开始,实际数据输出结果的前若干行可能全部是0,到时也得去掉。当然如果从0开始,到时去掉也行,多增加一步吧~·

因为对39个键值,我都是用循环:初始值=$min,步长=5000,终止值=$max  去卡它,这样就能保证每一列输出的数据都是一样长度的,而不管原始数据中每一个键值后边对应数据的不同。



各位帮帮忙~~
用其他语言,比如c++等也行

论坛徽章:
2
射手座
日期:2014-10-10 15:59:4715-16赛季CBA联赛之上海
日期:2016-03-03 10:27:14
9 [报告]
发表于 2011-10-29 11:59 |只看该作者
本帖最后由 yinyuemi 于 2011-10-29 12:59 编辑

回复 1# linghu002127


    大概意思,自己调整吧:
  1. #!/bin/perl

  2. use strict;
  3. my ($max_l,$min_l);
  4. my %hash;

  5. my @arr_title=("H2AK5ac","H2BK5me","H4K12ac","H4R3me2");

  6. open (file, "ChrY.txt");

  7. while(<file>){
  8.         chomp;
  9.         my @arr=split /\t/;
  10.         $hash{@arr[0]}.=@arr[1]." ";
  11.         $max_l=$max_l>$arr[1]?$max_l:$arr[1];
  12.         $min_l=$.==1||$min_l>$arr[1]?$arr[1]:$min_l;
  13. }

  14. printf "Name\t";
  15. print join "\t", @arr_title;
  16. print "\n";
  17. my $i;
  18. for ($i=0;$i<=$max_l;$i+=5000){
  19.         my @tmp;
  20.         foreach my $j (@arr_title){       
  21.                 my $n=0;
  22.                 map {if($_>=$i&&$_<$i+5000) {++$n} }(split " ",$hash{$j});
  23.                 push @tmp,$n;
  24.         }
  25.         print "chrY-",$i+5000,"\t";
  26.         print join "\t", @tmp;
  27.         @tmp=();
  28.         print "\n";
  29. }

  30. close(file);

复制代码

论坛徽章:
0
10 [报告]
发表于 2011-10-30 14:21 |只看该作者
回复 9# yinyuemi


    谢谢~~
您需要登录后才可以回帖 登录 | 注册

本版积分规则 发表回复

  

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

清除 Cookies - ChinaUnix - Archiver - WAP - TOP