求助！列内部区间的个数统计

summeraing

例如表 t 的数据如下：

ChrID   d1   d2   d3   d4   d5  d6
Chr1:1     A    B     A    C     N   A
Chr1:2     -    C     N    B     A   N
Chr1:5     B    -     -     -     N   A
Chr1:8     C    -     B    A     B   A
Chr1:10     A    -     -    N    A   C
Chr1:15     -    B     A    C     N   A
Chr1:16     A    C     N    B     A   N
Chr1:19     -    N     N    B     N   A
Chr1:22     C    -     B    -     A   B
Chr1:23     -    N     B    N    C   C
Chr1:25     A    -     A    -     N   A
Chr1:26     A    C     N    -     A   -
Chr1:27     B    N     C    -     -   A
Chr1:28     A    A     B    A     -   B
Chr1:29     A    C     B    N    C   C
Chr2:1     A    B     A    C     N   A
Chr2:2     -    C     N    B     A   N
Chr2:3     B    -     -     -     N   A
Chr2:5     C    -     B    A     B   A
Chr2:7     A    -     -    N    A   C
Chr2:10     -    B     A    C     N   A

表中 只有 A  B  C  N -   5种值
需要查出 表t中   分别求出每列的 A B C N在每个Chr的每10区间的个数？Chr的总数和所有总数

分别找 A B C N

例如结果如下：
  找A的区间个数和总数
A.txt
site        d1   d2   d3    d4   d5    d6
Chr1:1-10       2     0     1      1     2     3
Chr1:10-20     1     0     1      0     1     2
Chr1:20-30     4     1     1      1     2     2
Chr2:1-10     2     0     2      1     2     4


Chr1total     7     1     3      2     5     7
Chr2total     2     0     2      1     2     4
totalA          9     1     5      3     7     11


B.txt


C.txt


N.txt

jason680

回复 1# summeraing

1. 
   
2. $ perl gen_ABCN.pl
   
3. Usage  : gen_ABCN.pl Table_file
   
4. Example: gen_ABCN.pl t.txt
   
5. 
   
6. === table file ===
   
7. ChrID     d1   d2    d3   d4    d5  d6
   
8. Chr1:1     A    B     A    C     N   A
   
9. Chr1:2     -    C     N    B     A   N
   
10. ...
    
11. 
    
12. 
    
13. $ perl gen_ABCN.pl t.txt
    
14. Output file: A.txt
    
15. Output file: B.txt
    
16. Output file: C.txt
    
17. Output file: N.txt

复制代码

1. $ grep '.*' {A,B,C,N}.txt
   
2. A.txt:site             d1    d2    d3    d4    d5    d6
   
3. A.txt:Chr1:01-10        2     0     1     1     2     2
   
4. A.txt:Chr1:11-20        1     0     1     0     1     2
   
5. A.txt:Chr1:21-30        4     1     1     1     2     2
   
6. A.txt:Chr2:01-10        2     0     2     1     2     4
   
7. A.txt:
   
8. A.txt:Chr1total         7     1     3     2     5     6
   
9. A.txt:Chr2total         2     0     2     1     2     4
   
10. A.txt:totalA            9     1     5     3     7    10
    
11. B.txt:site             d1    d2    d3    d4    d5    d6
    
12. B.txt:Chr1:01-10        0     1     1     1     1     0
    
13. B.txt:Chr1:11-20        0     1     0     2     0     0
    
14. B.txt:Chr1:21-30        1     0     4     0     0     2
    
15. B.txt:Chr2:01-10        1     2     1     1     1     0
    
16. B.txt:
    
17. B.txt:Chr1total         1     2     5     3     1     2
    
18. B.txt:Chr2total         1     2     1     1     1     0
    
19. B.txt:totalB            2     4     6     4     2     2
    
20. C.txt:site             d1    d2    d3    d4    d5    d6
    
21. C.txt:Chr1:01-10        1     1     0     1     0     1
    
22. C.txt:Chr1:11-20        0     1     0     1     0     0
    
23. C.txt:Chr1:21-30        1     2     1     0     2     2
    
24. C.txt:Chr2:01-10        1     1     0     2     0     1
    
25. C.txt:
    
26. C.txt:Chr1total         2     4     1     2     2     3
    
27. C.txt:Chr2total         1     1     0     2     0     1
    
28. C.txt:totalC            3     5     1     4     2     4
    
29. N.txt:site             d1    d2    d3    d4    d5    d6
    
30. N.txt:Chr1:01-10        0     0     1     1     1     1
    
31. N.txt:Chr1:11-20        0     1     2     0     2     1
    
32. N.txt:Chr1:21-30        0     2     1     2     1     0
    
33. N.txt:Chr2:01-10        0     0     1     1     3     1
    
34. N.txt:
    
35. N.txt:Chr1total         0     3     4     3     4     2
    
36. N.txt:Chr2total         0     0     1     1     3     1
    
37. N.txt:totalN            0     3     5     4     7     3

复制代码

$ cat gen_ABCN.pl

1. use strict;
   
2. use warnings;
   
3. 
   
4. my $sFile_display = 0;
   
5. my $sItem_check = 1;
   
6. my @aItem = qw(A B C N -);
   
7. my %hItem = map {$_ => 1} @aItem;
   
8. 
   
9. my @aHead;
   
10. my %hChr;
    
11. my $sHead;
    
12. 
    
13. 
    
14. sub usage{
    
15.   print <<__USAGE__;
    
16. Usage  : $0 Table_file
    
17. Example: $0 t.txt
    
18. 
    
19. === table file ===
    
20. ChrID     d1   d2    d3   d4    d5  d6
    
21. Chr1:1     A    B     A    C     N   A
    
22. Chr1:2     -    C     N    B     A   N
    
23. ...
    
24. 
    
25. __USAGE__
    
26.   exit 1;
    
27. }
    
28. 
    
29. usage() if(@ARGV == 0);
    
30. 
    
31. my $sCnt_head = 0;
    
32. my $sCnt = 0;
    
33. my $sLen = 0;
    
34. while(<>){
    
35.   ++$sCnt;
    
36.   chomp;
    
37.   my $sLine = $_;
    
38.   s/^\s+|\s$//g;
    
39.   next if(m/^(#|$)/);
    
40. 
    
41.   # ChrID   d1   d2   d3   d4   d5  d6
    
42.   # Chr1:1   A    B    A    C    N   A
    
43.   my($sChr_num, @aData) = split;
    
44.   if(m/^ChrID\s/){
    
45.     @aHead  = @aData;
    
46.     $sLen  = scalar @aHead;
    
47.     $sHead = $sLine;
    
48.     $sCnt_head = $sCnt;
    
49.     next;
    
50.   }
    
51.   die "*** error: Cannot get 'ChrID' header\n" if($sLen == 0);
    
52.   if(scalar @aData != $sLen){
    
53.     print "*** error *** Line[$sCnt]: get wrong items number\n";
    
54.     printf("  Line[%03d]: $sHead\n",$sCnt_head);
    
55.     printf("  Line[%03d]: $sLine\n",$sCnt);
    
56.     # exit 1;
    
57.   }
    
58. 
    
59.   # Chr1:1   A    B    A    C    N   A
    
60.   #   Chr1 :  1
    
61.   my($sChr, $sNum) = split(/:/, $sChr_num);
    
62.   # @aData =  A    B    A    C    N   A
    
63.   foreach(0 .. $#aData){
    
64.     my $sItem = $aData[$_];
    
65.     if($sItem_check != 0 && (! exists $hItem{$sItem})){
    
66.       print "*** error *** Line[$sCnt]: get wrong item '$sItem'\n";
    
67.       print "  check items: @aItem\n";
    
68.       printf("  Line[%03d]: $sLine\n", $sCnt);
    
69.       # exit 1;
    
70.     }
    
71.     # for initial value (0)
    
72.     if(! exists $hChr{$sItem}{"Total"}{$sChr}){
    
73.       @{$hChr{$sItem}{"Total"}{$sChr}} = (0) x $sLen;
    
74.     }
    
75.     if(! exists $hChr{$sItem}{"Number"}{$sChr}{int(($sNum-1)/10)}){
    
76.       @{$hChr{$sItem}{"Number"}{$sChr}{int(($sNum-1)/10)}} = (0) x $sLen;
    
77.     }
    
78.     if(! exists $hChr{$sItem}{"All"}){
    
79.       @{$hChr{$sItem}{"All"}} = (0) x $sLen;
    
80.     }
    
81.     # for Chr1:1-10, Chr2:11-20, ...
    
82.     ++$hChr{$sItem}{"Number"}{$sChr}{int(($sNum-1)/10)}[$_];
    
83.     # for Chr1total, Chr2total
    
84.     ++$hChr{$sItem}{"Total"}{$sChr}[$_];
    
85.     # for totalA, totalB, ...
    
86.     ++$hChr{$sItem}{"All"}[$_];
    
87.   }  
    
88. }
    
89. 
    
90. my $sFmt = "%-12s " . (" %5s" x @aHead) . "\n";
    
91. foreach my $sItem(sort keys %hChr){
    
92.   next if($sItem eq "-");
    
93.   # for head, site  d1   d2   d3    d4   d5    d6
    
94.   my $sOut = sprintf($sFmt, "site", @aHead);
    
95.   foreach my $sChr(sort keys %{$hChr{$sItem}{"Number"}}){
    
96.     my $rhChr = $hChr{$sItem}{"Number"}{$sChr};
    
97.     # for Chr1:01-10 2 0 1 1 2 3
    
98.     # for Chr1:11-20 1 0 1 0 1 2
    
99.     foreach(sort{ $a<=>$b } keys %{$rhChr}){
    
100.       #    Chr1:01-10  <==     Chr1: 01 - 10
     
101.       my $sChr_num = sprintf("$sChr:%02d-%02d",$_*10+1,$_*10+10);
     
102.       $sOut .= sprintf($sFmt, $sChr_num, @{$rhChr->{$_}});
     
103.     }
     
104.   }
     
105.   $sOut .= "\n";
     
106.   # for Chr1total 7 1 3 2 5 7
     
107.   # for Chr2total ...
     
108.   foreach my $sChr(sort keys %{$hChr{$sItem}{"Total"}}){
     
109.     $sOut .= sprintf($sFmt, "${sChr}total", @{$hChr{$sItem}{"Total"}{$sChr}});
     
110.   }
     
111.   # for totalA ..., totalB ..., ...
     
112.   $sOut .= sprintf($sFmt, "total$sItem", @{$hChr{$sItem}{"All"}});
     
113. 
     
114.   # output file
     
115.   my $sFout = "$sItem.txt";
     
116.   open(my $FHout, ">", $sFout) or die "cannot open $sFout file\n";
     
117.   print {$FHout} $sOut;
     
118.   close $FHout;
     
119.   print "Output file: $sFout\n";
     
120.   print $sOut if($sFile_display);
     
121.   
     
122. }

复制代码

summeraing

回复 2# jason680

太感谢大神，完全满足我的格式，完美！~！

献上膝盖和打赏！~~~！

rubyish

perl abc.pl data.txt

1. #!/usr/bin/perl -w
   
2. use 5.010;
   
3. 
   
4. my $data = $ARGV[0] || 'data.txt';
   
5. open my $DATA, $data or die;
   
6. 
   
7. my ( $site, $tot, $X ) = qw[site total -];
   
8. my ( undef, @D ) = split /\s+/, <$DATA>;
   
9. my @O = (0) x @D;
   
10. my ( %CHR, %TOT, %TAL );
    
11. 
    
12. while (<$DATA>) {
    
13.     my ( $chri, @d ) = split;
    
14.     my ( $chr, $i ) = split ':', $chri;
    
15.     my $I = int( ( $i - 1 ) / 10 );
    
16.    
    
17.     for my $i ( 0 .. $#d ) {
    
18.         next if $d[$i] eq $X;
    
19.         
    
20.         $CHR{ $d[$i] }{$chr}[$I] //= [@O];
    
21.         $TOT{ $d[$i] }{$chr}     //= [@O];
    
22.         $TAL{ $d[$i] }           //= [@O];
    
23. 
    
24.         $CHR{ $d[$i] }{$chr}[$I][$i]++;
    
25.         $TOT{ $d[$i] }{$chr}[$i]++;
    
26.         $TAL{ $d[$i] }[$i]++;
    
27. 
    
28.     }
    
29. }
    
30. 
    
31. my ( undef, $K ) = %TOT;
    
32. my @ORD =
    
33.   map  { $_->[1] }
    
34.   sort { $a->[0] <=> $b->[0] }
    
35.   map  { [ /(\d+)/, $_ ] } keys %$K;
    
36. 
    
37. $site = join "\t", $site, @D;
    
38. 
    
39. while ( my ( $C, $V ) = each %CHR ) {
    
40.     open my $file, '>', $C . '.txt' or die;
    
41.     say STDERR "write\t$C.txt";
    
42.     say $file $site;
    
43.    
    
44.     for my $c (@ORD) {
    
45.         my @range = ( 1, '-', 10 );
    
46.         
    
47.         for my $d ( @{ $V->{$c} } ) {
    
48.             if ($d) {
    
49.                 my $head = $c . ':' . join( '', @range );
    
50.                 say $file join "\t", $head, @$d;
    
51.             }
    
52.             @range[ 0, -1 ] = ( $range[-1], $range[-1] + 10 );
    
53.         }
    
54.     }
    
55.     say $file "\n";
    
56. 
    
57.     for my $chr (@ORD) {
    
58.         say $file join "\t", $chr . $tot, @{ $TOT{$C}{$chr} };
    
59.     }
    
60.    
    
61.     say $file join "\t", $tot . $C, @{ $TAL{$C} };
    
62.     close $file;
    
63. }
    
64. 
    
65. say STDERR 'OK';
    
66. 
    
67. __DATA__
    
68. $_
    

复制代码