免费注册 查看新帖 |

Chinaunix

  平台 论坛 博客 文库
123下一页
最近访问板块 发新帖
查看: 43683 | 回复: 22

菜鸟求助 [复制链接]

论坛徽章:
0
发表于 2016-07-28 15:03 |显示全部楼层
本人初学编程,无奈遇到棘手问题,遍寻网络,只好在此求教。望高人略施援手,不盛感激之至。
部分txt示例文档如下:
chrX        hg19_refGene        exon        13336768        13338518        0        -        .        gene_id "NM_001135995"; transcript_id "NM_001135995";
chrX        hg19_refGene        exon        13353360        13353431        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
chrX        hg19_refGene        exon        13353601        13353664        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
chrX        hg19_refGene        exon        13353778        13353816        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
chrX        hg19_refGene        exon        13354498        13354579        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
chrX        hg19_refGene        exon        13356014        13356067        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
chrX        hg19_refGene        exon        13358045        13358218        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
chrX        hg19_refGene        exon        13359850        13359944        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
chrX        hg19_refGene        exon        13608411        13608465        0        +        .        gene_id "NR_106734"; transcript_id "NR_106734";
chrX        hg19_refGene        exon        13707240        13707407        0        +        .        gene_id "NM_001195328"; transcript_id "NM_001195328";
chrX        hg19_refGene        exon        13726840        13727944        0        +        .        gene_id "NM_001195328"; transcript_id "NM_001195328";
chrX        hg19_refGene        exon        13707240        13707407        0        +        .        gene_id "NM_004251"; transcript_id "NM_004251";
chrX        hg19_refGene        exon        13721933        13722021        0        +        .        gene_id "NM_004251"; transcript_id "NM_004251";
chrX        hg19_refGene        exon        13726840        13727944        0        +        .        gene_id "NM_004251"; transcript_id "NM_004251";
求助以第9列“gene_id*”的内容生成两个文件,
文件一,第9列内容仅出现一次的集合,并且输出第5列减去第4列数字的数值;
如第一行就是不重复的:
输出文件为:
chrX        hg19_refGene        exon        13336768        13338518        0        -        .        gene_id "NM_001135995"; transcript_id "NM_001135995";  1750
其中,1750 = 13338518-13336768;
第二个文件,同样依据第九列内容把出现重复的归到此类,并且在第10列输出重复次数;第11列输出每一次重复时,第5列减去第4列数值,相加总和再加1;第12列为重复的最后一行最大数值减去第一行重复最小数值后的数值。
如gene_id "NM_001195328",重复了2次,
第10列的数值是:1272 = 13707407-13707240+13727944-13726840+1
第11列的数值是:20704 = 13727944-13707240
因此,这个示例txt最终返回以下两个文件,即:
1.txt
chrX        hg19_refGene        exon        13336768        13338518        0        -        .        gene_id "NM_001135995"; transcript_id "NM_001135995";  1750
chrX        hg19_refGene        exon        13608411        13608465        0        +        .        gene_id "NR_106734"; transcript_id "NR_106734";  54

2.txt
chrX        hg19_refGene        exon        13353360        13353431        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
chrX        hg19_refGene        exon        13353601        13353664        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
chrX        hg19_refGene        exon        13353778        13353816        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
chrX        hg19_refGene        exon        13354498        13354579        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
chrX        hg19_refGene        exon        13356014        13356067        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
chrX        hg19_refGene        exon        13358045        13358218        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
chrX        hg19_refGene        exon        13359850        13359944        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";  7   579  6584
chrX        hg19_refGene        exon        13707240        13707407        0        +        .        gene_id "NM_001195328"; transcript_id "NM_001195328";
chrX        hg19_refGene        exon        13726840        13727944        0        +        .        gene_id "NM_001195328"; transcript_id "NM_001195328";  2  1272  20704
chrX        hg19_refGene        exon        13707240        13707407        0        +        .        gene_id "NM_004251"; transcript_id "NM_004251";
chrX        hg19_refGene        exon        13721933        13722021        0        +        .        gene_id "NM_004251"; transcript_id "NM_004251";
chrX        hg19_refGene        exon        13726840        13727944        0        +        .        gene_id "NM_004251"; transcript_id "NM_004251";   3   1361   20704

论坛徽章:
307
程序设计版块每周发帖之星
日期:2016-04-08 00:41:33操作系统版块每日发帖之星
日期:2015-09-02 06:20:00每日论坛发贴之星
日期:2015-09-02 06:20:00程序设计版块每日发帖之星
日期:2015-09-04 06:20:00每日论坛发贴之星
日期:2015-09-04 06:20:00每周论坛发贴之星
日期:2015-09-06 22:22:00程序设计版块每日发帖之星
日期:2015-09-09 06:20:00程序设计版块每日发帖之星
日期:2015-09-19 06:20:00程序设计版块每日发帖之星
日期:2015-09-20 06:20:00每日论坛发贴之星
日期:2015-09-20 06:20:00程序设计版块每日发帖之星
日期:2015-09-22 06:20:00程序设计版块每日发帖之星
日期:2015-09-24 06:20:00
发表于 2016-07-28 17:15 |显示全部楼层
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;

  4. my ($regex, %hData, $first) = qr/\b(\d+)\h+(\d+)\b/;
  5. /("[^"]+";)/ and push (@{$hData{$1}}, $_) while (<DATA>);

  6. foreach (sort {@{$hData{$a}} <=> @{$hData{$b}}} keys %hData){
  7.     chomp (my @aData = @{$hData{$_}});
  8.     my ($sum, $min, $max) = (0) x 3;
  9.     foreach my $offset (0 .. $#aData){
  10.         my ($n1, $n2) = $aData[$offset] =~ /$regex/;
  11.         $min = $n1 if (!$offset);
  12.         $max = $n2 if ($offset == $#aData);
  13.         $sum += $offset ? ($2 - $1 + 1) : ($2 - $1);
  14.     }
  15.     if (@aData == 1){
  16.         print "@aData $sum\n";
  17.         next;
  18.     }
  19.     print "\n\n" if (!$first++);
  20.     print join ("\n", @aData[0 .. $#aData-1], ""), "$aData[-1] ", scalar (@aData), " $sum ", ($max - $min), "\n";
  21. }

  22. __DATA__
  23. chrX        hg19_refGene        exon        13336768        13338518        0        -        .        gene_id "NM_001135995"; transcript_id "NM_001135995";
  24. chrX        hg19_refGene        exon        13353360        13353431        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
  25. chrX        hg19_refGene        exon        13353601        13353664        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
  26. chrX        hg19_refGene        exon        13353778        13353816        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
  27. chrX        hg19_refGene        exon        13354498        13354579        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
  28. chrX        hg19_refGene        exon        13356014        13356067        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
  29. chrX        hg19_refGene        exon        13358045        13358218        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
  30. chrX        hg19_refGene        exon        13359850        13359944        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
  31. chrX        hg19_refGene        exon        13608411        13608465        0        +        .        gene_id "NR_106734"; transcript_id "NR_106734";
  32. chrX        hg19_refGene        exon        13707240        13707407        0        +        .        gene_id "NM_001195328"; transcript_id "NM_001195328";
  33. chrX        hg19_refGene        exon        13726840        13727944        0        +        .        gene_id "NM_001195328"; transcript_id "NM_001195328";
  34. chrX        hg19_refGene        exon        13707240        13707407        0        +        .        gene_id "NM_004251"; transcript_id "NM_004251";
  35. chrX        hg19_refGene        exon        13721933        13722021        0        +        .        gene_id "NM_004251"; transcript_id "NM_004251";
  36. chrX        hg19_refGene        exon        13726840        13727944        0        +        .        gene_id "NM_004251"; transcript_id "NM_004251";
复制代码

论坛徽章:
0
发表于 2016-07-28 18:36 |显示全部楼层
非常感谢,试了一下,可以用。这几天快逼疯了,不是学计算机专业的,写几行极简单的代码还可以,这个完全干不了(目前还看不懂高手写的东东,当然我会努力学习的)
原始data.txt有49万多行,因为写不来代码,被逼用excel硬搞,算到30多万行,excel打开就崩溃
目前算了半个小时了,还没出结果;看来是个大工程呀!
万分感谢1楼

论坛徽章:
0
发表于 2016-07-28 18:45 |显示全部楼层
回复 2# sunzhiguolu

大神,再次请教,运行完结果出错耶!我只能看到最后一行,居然“加和”出现负值,看来程序还是有问题

能不能让电脑读取data.txt,然后直接print出结果1.txt和2.txt。我是很白很白的菜鸟。
   

论坛徽章:
307
程序设计版块每周发帖之星
日期:2016-04-08 00:41:33操作系统版块每日发帖之星
日期:2015-09-02 06:20:00每日论坛发贴之星
日期:2015-09-02 06:20:00程序设计版块每日发帖之星
日期:2015-09-04 06:20:00每日论坛发贴之星
日期:2015-09-04 06:20:00每周论坛发贴之星
日期:2015-09-06 22:22:00程序设计版块每日发帖之星
日期:2015-09-09 06:20:00程序设计版块每日发帖之星
日期:2015-09-19 06:20:00程序设计版块每日发帖之星
日期:2015-09-20 06:20:00每日论坛发贴之星
日期:2015-09-20 06:20:00程序设计版块每日发帖之星
日期:2015-09-22 06:20:00程序设计版块每日发帖之星
日期:2015-09-24 06:20:00
发表于 2016-07-28 18:49 |显示全部楼层
回复 3# jacqueslm2001
问下, 实际的文件大小?

   

论坛徽章:
307
程序设计版块每周发帖之星
日期:2016-04-08 00:41:33操作系统版块每日发帖之星
日期:2015-09-02 06:20:00每日论坛发贴之星
日期:2015-09-02 06:20:00程序设计版块每日发帖之星
日期:2015-09-04 06:20:00每日论坛发贴之星
日期:2015-09-04 06:20:00每周论坛发贴之星
日期:2015-09-06 22:22:00程序设计版块每日发帖之星
日期:2015-09-09 06:20:00程序设计版块每日发帖之星
日期:2015-09-19 06:20:00程序设计版块每日发帖之星
日期:2015-09-20 06:20:00每日论坛发贴之星
日期:2015-09-20 06:20:00程序设计版块每日发帖之星
日期:2015-09-22 06:20:00程序设计版块每日发帖之星
日期:2015-09-24 06:20:00
发表于 2016-07-28 18:52 |显示全部楼层
另外, 按照实际数据文件的格式贴出一部分示例文本. 出现问题的原因可能与文本格式有关.

论坛徽章:
0
发表于 2016-07-28 18:58 |显示全部楼层
是一个有493587行,*9列的文本;大小48521KB。
数据还有一个特点,就是第9列“gene_id***”排列是规则的,不会第1-10行出现了,在第10001-10010行又会冒出来,这个应该能缩写运算量

论坛徽章:
307
程序设计版块每周发帖之星
日期:2016-04-08 00:41:33操作系统版块每日发帖之星
日期:2015-09-02 06:20:00每日论坛发贴之星
日期:2015-09-02 06:20:00程序设计版块每日发帖之星
日期:2015-09-04 06:20:00每日论坛发贴之星
日期:2015-09-04 06:20:00每周论坛发贴之星
日期:2015-09-06 22:22:00程序设计版块每日发帖之星
日期:2015-09-09 06:20:00程序设计版块每日发帖之星
日期:2015-09-19 06:20:00程序设计版块每日发帖之星
日期:2015-09-20 06:20:00每日论坛发贴之星
日期:2015-09-20 06:20:00程序设计版块每日发帖之星
日期:2015-09-22 06:20:00程序设计版块每日发帖之星
日期:2015-09-24 06:20:00
发表于 2016-07-28 19:12 |显示全部楼层
本帖最后由 sunzhiguolu 于 2016-07-28 19:14 编辑
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;

  4. my ($regex, %hData) = qr/\b(\d+)\h+(\d+)\b/;
  5. /("[^"]+";)/ and push (@{$hData{$1}}, $_) while (<>);

  6. open (my $fhF1, '>', './1.txt');
  7. open (my $fhF2, '>', './2.txt');

  8. foreach (keys %hData){
  9.     chomp (my @aData = @{$hData{$_}});
  10.     my ($sum, $min, $max) = (0) x 3;
  11.     foreach my $offset (0 .. $#aData){
  12.         my ($n1, $n2) = $aData[$offset] =~ /$regex/;
  13.         $min = $n1 if (!$offset);
  14.         $max = $n2 if ($offset == $#aData);
  15.         $sum += $offset ? ($2 - $1 + 1) : ($2 - $1);
  16.     }
  17.     if (@aData == 1){
  18.         print $fhF1 "@aData $sum\n";
  19.         next;
  20.     }
  21.     print $fhF2 join ("\n", @aData[0 .. $#aData-1], ""), "$aData[-1] ", scalar (@aData), " $sum ", ($max - $min), "\n";
  22. }

  23. close ($fhF1);
  24. close ($fhF2);
复制代码
perl ab.pl file

cat 1.txt
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
chrX        hg19_refGene        exon        13336768        13338518        0        -        .        gene_id "NM_001135995"; transcript_id "NM_001135995"; 1750
chrX        hg19_refGene        exon        13608411        13608465        0        +        .        gene_id "NR_106734"; transcript_id "NR_106734"; 54

cat 2.txt
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
chrX        hg19_refGene        exon        13707240        13707407        0        +        .        gene_id "NM_001195328"; transcript_id "NM_001195328";
chrX        hg19_refGene        exon        13726840        13727944        0        +        .        gene_id "NM_001195328"; transcript_id "NM_001195328"; 2 1272 20704
chrX        hg19_refGene        exon        13707240        13707407        0        +        .        gene_id "NM_004251"; transcript_id "NM_004251";
chrX        hg19_refGene        exon        13721933        13722021        0        +        .        gene_id "NM_004251"; transcript_id "NM_004251";
chrX        hg19_refGene        exon        13726840        13727944        0        +        .        gene_id "NM_004251"; transcript_id "NM_004251"; 3 1361 20704
chrX        hg19_refGene        exon        13353360        13353431        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
chrX        hg19_refGene        exon        13353601        13353664        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
chrX        hg19_refGene        exon        13353778        13353816        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
chrX        hg19_refGene        exon        13354498        13354579        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
chrX        hg19_refGene        exon        13356014        13356067        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
chrX        hg19_refGene        exon        13358045        13358218        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
chrX        hg19_refGene        exon        13359850        13359944        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260"; 7 579 6584

论坛徽章:
0
发表于 2016-07-28 22:08 |显示全部楼层
非常非常感谢,解决了我的一个大麻烦。
我也查出来,为啥开始第一次会出现结果错误。数据里面存在一类第9列一致,但是第1列不一致的极个别的例子;因为有49万行,开始不知道有这样的特例在里面。所以,这次的结果还是有错,但是应该这样的特例不多,至少excel应该能找出来;此外,我也试着看看能不能用perl解决。

还是特别感谢大神!

论坛徽章:
0
发表于 2016-07-29 11:29 |显示全部楼层
sunzhiguolu 发表于 2016-07-28 17:15

昨晚试了半天,加了些小零件,把输出格式部分调整到最理想了。但是核心问题,没有解决;
我发现49万数据里面,第1列和第9列都可以出现变化。
诸如:
chr6        hg19_refGene        exon        31804853        31804916        0        +        .        gene_id "NR_002742"; transcript_id "NR_002742";
chr6_apd_hap1        hg19_refGene        exon        3119608        3119671        0        +        .        gene_id "NR_002742"; transcript_id "NR_002742";
chr6_cox_hap2        hg19_refGene        exon        3314427        3314490        0        +        .        gene_id "NR_002742"; transcript_id "NR_002742";
chr6_dbb_hap3        hg19_refGene        exon        3090439        3090502        0        +        .        gene_id "NR_002742"; transcript_id "NR_002742";
chr6_qbl_hap6        hg19_refGene        exon        3098509        3098572        0        +        .        gene_id "NR_002742"; transcript_id "NR_002742";
chr4        hg19_refGene        exon        9350619        9352211        0        +        .        gene_id "NM_001242331"; transcript_id "NM_001242331";
chr4        hg19_refGene        exon        9355364        9356956        0        +        .        gene_id "NM_001242331"; transcript_id "NM_001242331_dup1";
chr4        hg19_refGene        exon        9364855        9366447        0        +        .        gene_id "NM_001242331"; transcript_id "NM_001242331_dup2";
chr4        hg19_refGene        exon        9326891        9328483        0        +        .        gene_id "NM_001242331"; transcript_id "NM_001242331_dup3";
chr4        hg19_refGene        exon        9331637        9333229        0        +        .        gene_id "NM_001242331"; transcript_id "NM_001242331_dup4";
chr4        hg19_refGene        exon        9336384        9337976        0        +        .        gene_id "NM_001242331"; transcript_id "NM_001242331_dup5";
chr4        hg19_refGene        exon        9341129        9342721        0        +        .        gene_id "NM_001242331"; transcript_id "NM_001242331_dup6";
chr4        hg19_refGene        exon        9345874        9347466        0        +        .        gene_id "NM_001242331"; transcript_id "NM_001242331_dup7";
即,第1列有的内容是"字符串_字符串”;第9列是最后一个引号内“字符串_dup*”;这里的例子是单一,很多时候也是重复的,即“chr*_hap1”第二行仍是”chr*_hap1"。第9列也是一样,可能好几个“*_dup1”。数据仍是规律排列。

只能再麻烦大神帮忙修一下代码,不同列的多重比较,对我实在是太难;我贴一个把主贴和新问题集合在一起的test.txt数据文件,希望您不吝赐教。

test.rar

703 Bytes, 下载次数: 2

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

本版积分规则 发表回复

  

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

清除 Cookies - ChinaUnix - Archiver - WAP - TOP