Chinaunix

标题: 菜鸟求助 [打印本页]
作者: jacqueslm2001    时间: 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
作者: sunzhiguolu    时间: 2016-07-28 17:15
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;

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

  8. foreach (sort {@{$hData{$a}} <=> @{$hData{$b}}} 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 "@aData $sum\n";
  19.         next;
  20.     }
  21.     print "\n\n" if (!$first++);
  22.     print join ("\n", @aData[0 .. $#aData-1], ""), "$aData[-1] ", scalar (@aData), " $sum ", ($max - $min), "\n";
  23. }

  25. __DATA__
  26. chrX        hg19_refGene        exon        13336768        13338518        0        -        .        gene_id "NM_001135995"; transcript_id "NM_001135995";
  27. chrX        hg19_refGene        exon        13353360        13353431        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
  28. chrX        hg19_refGene        exon        13353601        13353664        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
  29. chrX        hg19_refGene        exon        13353778        13353816        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
  30. chrX        hg19_refGene        exon        13354498        13354579        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
  31. chrX        hg19_refGene        exon        13356014        13356067        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
  32. chrX        hg19_refGene        exon        13358045        13358218        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
  33. chrX        hg19_refGene        exon        13359850        13359944        0        +        .        gene_id "NR_045260"; transcript_id "NR_045260";
  34. chrX        hg19_refGene        exon        13608411        13608465        0        +        .        gene_id "NR_106734"; transcript_id "NR_106734";
  35. chrX        hg19_refGene        exon        13707240        13707407        0        +        .        gene_id "NM_001195328"; transcript_id "NM_001195328";
  36. chrX        hg19_refGene        exon        13726840        13727944        0        +        .        gene_id "NM_001195328"; transcript_id "NM_001195328";
  37. chrX        hg19_refGene        exon        13707240        13707407        0        +        .        gene_id "NM_004251"; transcript_id "NM_004251";
  38. chrX        hg19_refGene        exon        13721933        13722021        0        +        .        gene_id "NM_004251"; transcript_id "NM_004251";
  39. chrX        hg19_refGene        exon        13726840        13727944        0        +        .        gene_id "NM_004251"; transcript_id "NM_004251";
复制代码
作者: jacqueslm2001    时间: 2016-07-28 18:36
非常感谢，试了一下，可以用。这几天快逼疯了，不是学计算机专业的，写几行极简单的代码还可以，这个完全干不了（目前还看不懂高手写的东东，当然我会努力学习的）
原始data.txt有49万多行，因为写不来代码，被逼用excel硬搞，算到30多万行，excel打开就崩溃
目前算了半个小时了，还没出结果；看来是个大工程呀！
万分感谢1楼
作者: jacqueslm2001    时间: 2016-07-28 18:45
回复 2# sunzhiguolu

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

能不能让电脑读取data.txt，然后直接print出结果1.txt和2.txt。我是很白很白的菜鸟。
   
作者: sunzhiguolu    时间: 2016-07-28 18:49
回复 3# jacqueslm2001
问下, 实际的文件大小?

   
作者: sunzhiguolu    时间: 2016-07-28 18:52
另外, 按照实际数据文件的格式贴出一部分示例文本. 出现问题的原因可能与文本格式有关.
作者: jacqueslm2001    时间: 2016-07-28 18:58
是一个有493587行，*9列的文本；大小48521KB。
数据还有一个特点，就是第9列“gene_id***”排列是规则的，不会第1-10行出现了，在第10001-10010行又会冒出来，这个应该能缩写运算量
作者: sunzhiguolu    时间: 2016-07-28 19:12
本帖最后由 sunzhiguolu 于 2016-07-28 19:14 编辑
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;

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

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

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

  27. close ($fhF1);
  28. 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
作者: jacqueslm2001    时间: 2016-07-28 22:08
非常非常感谢，解决了我的一个大麻烦。
我也查出来，为啥开始第一次会出现结果错误。数据里面存在一类第9列一致，但是第1列不一致的极个别的例子；因为有49万行，开始不知道有这样的特例在里面。所以，这次的结果还是有错，但是应该这样的特例不多，至少excel应该能找出来；此外，我也试着看看能不能用perl解决。

还是特别感谢大神！
作者: jacqueslm2001    时间: 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
作者: sunzhiguolu    时间: 2016-07-29 11:44
即，第1列有的内容是"字符串_字符串”；第9列是最后一个引号内“字符串_dup*”；这里的例子是单一，很多时候也是重复的，即“chr*_hap1”第二行仍是”chr*_hap1"。第9列也是一样，可能好几个“*_dup1”。数据仍是规律排列。

你将特殊的情况一同贴出一小部分, 然后将希望的到的结果贴出来瞧瞧. 看的更直观些.
作者: jacqueslm2001    时间: 2016-07-29 12:25
我刚刚贴的附件里面test.rar里面的txt文件，所有类型都包含了。
8-15行，是第9列“*_dup1”的例子，有、无“_dup1”各重复4次；21-28行是“*_dup(1..7)”单一不重复的例子。
最后14行，是第1列变化类型，有、无“_*_*”各重复7次；16-20行是不重复的例子。
因此，test.txt已经包含了所有的情况。

非常感谢您的及时回复。
作者: sunzhiguolu    时间: 2016-07-29 13:19
本帖最后由 sunzhiguolu 于 2016-07-29 13:20 编辑
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;

  5. my ($regex, %hData) = qr/\b(\d+)\h+(\d+)\b/;
  6. push (@{$hData{qq(@{[split]}[-1])}}, $_) while (<>);

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

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

  27. close ($fhF1);
  28. close ($fhF2);
复制代码
perl abc.pl file
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

cat 1.txt
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
chr4    hg19_refGene    exon    9345874 9347466 0       +       .       gene_id "NM_001242331"; transcript_id "NM_001242331_dup7";  1592
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
chr4    hg19_refGene    exon    9350619 9352211 0       +       .       gene_id "NM_001242331"; transcript_id "NM_001242331";  1592
chr4    hg19_refGene    exon    9326891 9328483 0       +       .       gene_id "NM_001242331"; transcript_id "NM_001242331_dup3";  1592
chr4    hg19_refGene    exon    9364855 9366447 0       +       .       gene_id "NM_001242331"; transcript_id "NM_001242331_dup2";  1592
chr4    hg19_refGene    exon    9331637 9333229 0       +       .       gene_id "NM_001242331"; transcript_id "NM_001242331_dup4";  1592
chr4    hg19_refGene    exon    9355364 9356956 0       +       .       gene_id "NM_001242331"; transcript_id "NM_001242331_dup1";  1592
chr4    hg19_refGene    exon    9336384 9337976 0       +       .       gene_id "NM_001242331"; transcript_id "NM_001242331_dup5";  1592
chr4    hg19_refGene    exon    9341129 9342721 0       +       .       gene_id "NM_001242331"; transcript_id "NM_001242331_dup6";  1592

cat 2.txt
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
chrX    hg19_refGene    exon    119206241       119206810       0       -       .       gene_id "NM_032498"; transcript_id "NM_032498";
chrX    hg19_refGene    exon    119209865       119209910       0       -       .       gene_id "NM_032498"; transcript_id "NM_032498";
chrX    hg19_refGene    exon    119210842       119211253       0       -       .       gene_id "NM_032498"; transcript_id "NM_032498";
chrX    hg19_refGene    exon    119211439       119211707       0       -       .       gene_id "NM_032498"; transcript_id "NM_032498";  4 1296 5466
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    119292467       119292735       0       +       .       gene_id "NM_032498"; transcript_id "NM_032498_dup1";
chrX    hg19_refGene    exon    119292921       119293332       0       +       .       gene_id "NM_032498"; transcript_id "NM_032498_dup1";
chrX    hg19_refGene    exon    119294264       119294309       0       +       .       gene_id "NM_032498"; transcript_id "NM_032498_dup1";
chrX    hg19_refGene    exon    119297364       119297933       0       +       .       gene_id "NM_032498"; transcript_id "NM_032498_dup1";  4 1296 5466
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
chr6    hg19_refGene    exon    31633657        31633985        0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6    hg19_refGene    exon    31634598        31634680        0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6    hg19_refGene    exon    31635645        31635747        0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6    hg19_refGene    exon    31636316        31636431        0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6    hg19_refGene    exon    31636874        31636949        0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6    hg19_refGene    exon    31637105        31637285        0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6    hg19_refGene    exon    31637613        31637847        0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6_cox_hap2   hg19_refGene    exon    3143277 3143605 0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6_cox_hap2   hg19_refGene    exon    3144218 3144300 0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6_cox_hap2   hg19_refGene    exon    3145265 3145367 0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6_cox_hap2   hg19_refGene    exon    3145936 3146051 0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6_cox_hap2   hg19_refGene    exon    3146494 3146569 0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6_cox_hap2   hg19_refGene    exon    3146725 3146905 0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6_cox_hap2   hg19_refGene    exon    3147233 3147467 0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";  14 2245 4190
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";  5 319 63
作者: jacqueslm2001    时间: 2016-07-29 13:25
sunzhiguolu 发表于 2016-07-29 13:19
perl abc.pl file
---------------------------------------------------------------------------------- ...

还是不对，第1列不一样，应该也区分开。也就是数据要用第1列先判断，再用第9列判断，用双重标准回归最后的两个文件
还需要再麻烦您
作者: sunzhiguolu    时间: 2016-07-29 13:29
你看下结果是不是你想要的那样, 我只是按照你的要求 将最后一对 双引号内的字符串作为 ID 进行分类统计的.
作者: sunzhiguolu    时间: 2016-07-29 13:32
回复 14# jacqueslm2001
你所说的双重标准不是很理解, 能否将这个示例文件预期结果贴出来看下 双重标准到底是个啥样?

   
作者: jacqueslm2001    时间: 2016-07-29 13:55
sunzhiguolu 发表于 2016-07-29 13:32
回复 14# jacqueslm2001
你所说的双重标准不是很理解, 能否将这个示例文件预期结果贴出来看下 双重标准[/ ...

比如：
chr6    hg19_refGene    exon    31633657        31633985        0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6    hg19_refGene    exon    31634598        31634680        0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6    hg19_refGene    exon    31635645        31635747        0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6    hg19_refGene    exon    31636316        31636431        0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6    hg19_refGene    exon    31636874        31636949        0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6    hg19_refGene    exon    31637105        31637285        0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6    hg19_refGene    exon    31637613        31637847        0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385"; 7 1122 4190
chr6_cox_hap2   hg19_refGene    exon    3143277 3143605 0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6_cox_hap2   hg19_refGene    exon    3144218 3144300 0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6_cox_hap2   hg19_refGene    exon    3145265 3145367 0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6_cox_hap2   hg19_refGene    exon    3145936 3146051 0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6_cox_hap2   hg19_refGene    exon    3146494 3146569 0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6_cox_hap2   hg19_refGene    exon    3146725 3146905 0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6_cox_hap2   hg19_refGene    exon    3147233 3147467 0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";  7 1122 4190
还有最后5行，且应该回归到1.txt：
chr6    hg19_refGene    exon    31804853        31804916        0       +       .       gene_id "NR_002742"; transcript_id "NR_002742"; 1 63 63
chr6_apd_hap1   hg19_refGene    exon    3119608 3119671 0       +       .       gene_id "NR_002742"; transcript_id "NR_002742"; 1 63 63
chr6_cox_hap2   hg19_refGene    exon    3314427 3314490 0       +       .       gene_id "NR_002742"; transcript_id "NR_002742"; 1 63 63
chr6_dbb_hap3   hg19_refGene    exon    3090439 3090502 0       +       .       gene_id "NR_002742"; transcript_id "NR_002742"; 1 63 63
chr6_qbl_hap6   hg19_refGene    exon    3098509 3098572 0       +       .       gene_id "NR_002742"; transcript_id "NR_002742";  1 63 63

最初以为只是第9列会不同，结果运算完回到数据本身发现第1列也不同，因此先要比对第1列是不是一样，再比对第9列，问题变成这么复杂以后，我感觉自己能力完全不具备修改您的代码，实现最终的正确结果的输出
真的麻烦您了
作者: sunzhiguolu    时间: 2016-07-29 14:05
本帖最后由 sunzhiguolu 于 2016-07-29 14:17 编辑
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;

  5. my ($regex, %hData) = qr/\b(\d+)\h+(\d+)\b/;
  6. push (@{$hData{qq(@{[split]}[0,-1])}}, $_) while (<>);

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

  11. foreach (keys %hData){
  12.     chomp (my @aData = @{$hData{$_}});
  13.     my ($sum, $min, $max) = (0, 9e9, 0);
  14.     foreach my $offset (0 .. $#aData){
  15.         my ($n1, $n2) = $aData[$offset] =~ /$regex/;
  16.         my @aSort = sort {$a <=> $b} ($n1, $n2);
  17.         $min = $aSort[0] if ($aSort[0] < $min);
  18.         $max = $aSort[-1] if ($aSort[-1] > $max);
  19.         $sum += $offset ? ($2 - $1 + 1) : ($2 - $1);
  20.     }
  21.     if (@aData == 1){
  22.         print $fhF1 "@aData $sum\n";
  23.         next;
  24.     }
  25.     print $fhF2 join ("\n", @aData[0 .. $#aData-1], ""), "$aData[-1] ", scalar (@aData), " $sum ", ($max - $min), "\n";
  26. }

  28. close ($fhF1);
  29. close ($fhF2);
复制代码
perl abc.pl file

cat 1.txt
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
chr4    hg19_refGene    exon    9331637 9333229 0       +       .       gene_id "NM_001242331"; transcript_id "NM_001242331_dup4";  1592
chrX        hg19_refGene        exon        13608411        13608465        0        +        .        gene_id "NR_106734"; transcript_id "NR_106734";  54
chr6_cox_hap2   hg19_refGene    exon    3314427 3314490 0       +       .       gene_id "NR_002742"; transcript_id "NR_002742";  63
chr4    hg19_refGene    exon    9326891 9328483 0       +       .       gene_id "NM_001242331"; transcript_id "NM_001242331_dup3";  1592
chr4    hg19_refGene    exon    9336384 9337976 0       +       .       gene_id "NM_001242331"; transcript_id "NM_001242331_dup5";  1592
chr6_dbb_hap3   hg19_refGene    exon    3090439 3090502 0       +       .       gene_id "NR_002742"; transcript_id "NR_002742";  63
chr4    hg19_refGene    exon    9364855 9366447 0       +       .       gene_id "NM_001242331"; transcript_id "NM_001242331_dup2";  1592
chr4    hg19_refGene    exon    9355364 9356956 0       +       .       gene_id "NM_001242331"; transcript_id "NM_001242331_dup1";  1592
chr4    hg19_refGene    exon    9350619 9352211 0       +       .       gene_id "NM_001242331"; transcript_id "NM_001242331";  1592
chr6_apd_hap1   hg19_refGene    exon    3119608 3119671 0       +       .       gene_id "NR_002742"; transcript_id "NR_002742";  63
chr4    hg19_refGene    exon    9345874 9347466 0       +       .       gene_id "NM_001242331"; transcript_id "NM_001242331_dup7";  1592
chr6_qbl_hap6   hg19_refGene    exon    3098509 3098572 0       +       .       gene_id "NR_002742"; transcript_id "NR_002742";  63
chr4    hg19_refGene    exon    9341129 9342721 0       +       .       gene_id "NM_001242331"; transcript_id "NM_001242331_dup6";  1592
chrX        hg19_refGene        exon        13336768        13338518        0        -        .        gene_id "NM_001135995"; transcript_id "NM_001135995";  1750
chr6    hg19_refGene    exon    31804853        31804916        0       +       .       gene_id "NR_002742"; transcript_id "NR_002742";  63

cat 2.txt
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
chrX    hg19_refGene    exon    119206241       119206810       0       -       .       gene_id "NM_032498"; transcript_id "NM_032498";
chrX    hg19_refGene    exon    119209865       119209910       0       -       .       gene_id "NM_032498"; transcript_id "NM_032498";
chrX    hg19_refGene    exon    119210842       119211253       0       -       .       gene_id "NM_032498"; transcript_id "NM_032498";
chrX    hg19_refGene    exon    119211439       119211707       0       -       .       gene_id "NM_032498"; transcript_id "NM_032498";  4 1296 5466
chr6_cox_hap2   hg19_refGene    exon    3143277 3143605 0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6_cox_hap2   hg19_refGene    exon    3144218 3144300 0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6_cox_hap2   hg19_refGene    exon    3145265 3145367 0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6_cox_hap2   hg19_refGene    exon    3145936 3146051 0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6_cox_hap2   hg19_refGene    exon    3146494 3146569 0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6_cox_hap2   hg19_refGene    exon    3146725 3146905 0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6_cox_hap2   hg19_refGene    exon    3147233 3147467 0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";  7 1122 4190
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    119292467       119292735       0       +       .       gene_id "NM_032498"; transcript_id "NM_032498_dup1";
chrX    hg19_refGene    exon    119292921       119293332       0       +       .       gene_id "NM_032498"; transcript_id "NM_032498_dup1";
chrX    hg19_refGene    exon    119294264       119294309       0       +       .       gene_id "NM_032498"; transcript_id "NM_032498_dup1";
chrX    hg19_refGene    exon    119297364       119297933       0       +       .       gene_id "NM_032498"; transcript_id "NM_032498_dup1";  4 1296 5466
chr6    hg19_refGene    exon    31633657        31633985        0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6    hg19_refGene    exon    31634598        31634680        0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6    hg19_refGene    exon    31635645        31635747        0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6    hg19_refGene    exon    31636316        31636431        0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6    hg19_refGene    exon    31636874        31636949        0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6    hg19_refGene    exon    31637105        31637285        0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";
chr6    hg19_refGene    exon    31637613        31637847        0       +       .       gene_id "NM_001282385"; transcript_id "NM_001282385";  7 1122 4190
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
作者: sunzhiguolu    时间: 2016-07-29 14:18
你再测试下结果是否正确,,,
作者: jacqueslm2001    时间: 2016-07-29 14:54
回复 19# sunzhiguolu

粗略扫过，excel再筛选检查，应该是没错了，感谢之至。
作者: jacqueslm2001    时间: 2016-08-01 03:12
回复 19# sunzhiguolu
继续碰到问题，整个N*9列的数据太过复杂。我用您写的代码，分析数据后发现大概有400多个例子（大概8000多行）仍然属于复杂情况，于是我写了很粗糙的代码剥离这部分数据。如下，见笑：
open IN,"< Ref.txt";
@Ref = <IN>;
close IN;

open IN,"< data.txt";
@Data = <IN>;
close IN;

open (my $fhF1, '>', './Seperative.txt');
foreach $Data(@Data){
                chomp $Data;
                @temp2 = split(/\t/,$Data);
                foreach $Ref(@Ref){
                                chomp $Ref;
                                @temp1 = split(/\t/,$Ref);
                                if($temp2[0] eq $temp1[0] && $temp2[-1] eq $temp1[-1]){
                                        print $fhF1 "@temp2\n";
                                        next;
                                }
                }
}

close ($fhF1);
就是，我对data.txt，如下：
chr1        hg19_refGene        exon        154474695        154475505        0        +        .        gene_id "NM_001098475"; transcript_id "NM_001098475";
chr1        hg19_refGene        start_codon        154479380        154479381        0        +        .        gene_id "NM_001098475"; transcript_id "NM_001098475";
chr1        hg19_refGene        CDS        154479380        154479381        0        +        0        gene_id "NM_001098475"; transcript_id "NM_001098475";
chr1        hg19_refGene        exon        154479353        154479381        0        +        .        gene_id "NM_001098475"; transcript_id "NM_001098475";
chr1        hg19_refGene        stop_codon        154522914        154522916        0        -        .        gene_id "NM_017582"; transcript_id "NM_017582";
chr1        hg19_refGene        CDS        154522917        154522945        0        -        2        gene_id "NM_017582"; transcript_id "NM_017582";
chr1        hg19_refGene        exon        154521051        154522945        0        -        .        gene_id "NM_017582"; transcript_id "NM_017582";
chr1        hg19_refGene        start_codon        161087967        161087968        0        +        .        gene_id "NM_001185092"; transcript_id "NM_001185092";
chr1        hg19_refGene        CDS        161087967        161087968        0        +        0        gene_id "NM_001185092"; transcript_id "NM_001185092";
chr1        hg19_refGene        exon        161087862        161087968        0        +        .        gene_id "NM_001185092"; transcript_id "NM_001185092";
设计了一个，ref.txt，如下：
chr1        hg19_refGene        exon        154474695        154475505        0        +        .        gene_id "NM_001098475"; transcript_id "NM_001098475";
chr1        hg19_refGene        exon        161087862        161087968        0        +        .        gene_id "NM_001185092"; transcript_id "NM_001185092";
通过，这条“$temp2[0] eq $temp1[0] && $temp2[-1] eq $temp1[-1]”表示只有第1列，第9列都相同的数据，分割出来（Seperative.txt）。原本我的目的是直接生成两个txt文件（都相同的生产一个文件；剩下的生产另一个文件），可是只有这个都形同的结果输出可以，是正确的7行；但是else的结果不对，不是3行数据（没贴在这儿）。
此外，即便正确的，输出的数组原来的“\t”格式不知道为什么只被识别成了空格（9列数据变成合在一起的1列数据输出了），修改了好多遍，解决不了。其它辅助方法，通过UE将空格转“\t”，由于[-1]gene**后面好几个空格，结果第9列被拆成了好几列。
1.我想学习下，如果是您会如何设计这样的代码？
2.显然，我的代码拆分两个txt就会出错；而且格式被修改了。因此，只能再次求助了
作者: sunzhiguolu    时间: 2016-08-01 09:01
本帖最后由 sunzhiguolu 于 2016-08-01 09:05 编辑

你再测试下,
open IN,"< Ref.txt";
@Ref = <IN>;
close IN;

open IN,"< data.txt";
@Data = <IN>;
close IN;

open (my $fhF1, '>', './Seperative.txt');
foreach $Data(@Data){
    chomp $Data;
    @temp2 = split(' ',$Data);
    $match = 0;
    foreach $Ref(@Ref){
        chomp $Ref;
        @temp1 = split(' ',$Ref);
        $match = 1 if($temp2[0] eq $temp1[0] && $temp2[-1] eq $temp1[-1]);
    }
    if ($match == 1){
        print $fhF1 join ("\t", @temp2), "\n";
        next;
    }
    print join ("\t", @temp2), "\n";
}
作者: jacqueslm2001    时间: 2016-08-01 12:34
稍微修了一点点，可以用了；开始我也用了一个check=0，后来想想好像这个判断可有可无，最后就删掉了，结果。。。哎！还是凭着20年前basic的底子，真的干不动
谢谢



欢迎光临 Chinaunix (http://bbs.chinaunix.net/) Powered by Discuz! X3.2