【求助改程序高手帮我改下这个程序】，只需输出结果中加入一些数据

yanglc2013

因为他输出的格式是（只有lnc中的数据）：
lnc-C1orf227-1:2        chr1        213028819        213031452        -

但我想在同行加上对应的refGene中的数据（附件中，下面的程序也是基于附件中的两个文件的比较。）
lnc-C1orf227-1:2        chr1        213028819        213031452        -      NM_014053        chr1        213031596        213072705        +        2210  （overlap关系）

原帖地址：http://bbs.chinaunix.net/thread-4073924-1-1.html（一下程序已经解决原帖的问题，现在想输出结果把refGene中对应的数据输出）

1. #!/usr/bin/perl -w
   
2. use strict;
   
3. use Data::Dumper;
   
4. my($ref,$lnc,$contain,$overlap,$gapped,$reverse);
   
5. open($ref,$ARGV[0]) or die $!;
   
6. open($lnc,$ARGV[1]) or die $!;
   
7. open($contain,'>','contain.txt') or die $!;
   
8. open($overlap,'>','overlap.txt') or die $!;
   
9. open($gapped ,'>','gapped.txt')  or die $!;
   
10. open($reverse,'>','reverse.txt') or die $!;
    
11. 
    
12. my %hRefDB;
    
13. while (<$ref>){
    
14.    my ($chr, $beg, $end) = (split)[1,2,3]; #Tab拆分，并取第二，第三，第四个元素。
    
15.       push @{ $hRefDB{$chr}}, [$beg,$end]; #存入哈希数组。
    
16. }
    
17. 
    
18. while (<$lnc>) {
    
19.    my ($chr,$beg,$end) = (split)[1,2,3];    #Tab拆分，并取第二，第三，第四个元素。
    
20.    if (!exists $hRefDB{$chr}) {            
    
21.        print "not exists in refDB\n";
    
22.        next;
    
23.    }
    
24.    my $lnc = [$beg,$end];
    
25.    my ($is_contain,$is_overlap,$is_reverse) = (0,0,0);
    
26.    foreach my $ref (@{$hRefDB{$chr}}){       #遍历哈希数组的每个元素。
    
27.        last if $is_contain and $is_overlap and $is_reverse;
    
28.        unless ($is_contain){                 #除非$is_contain的值为0，否则，将is_contain($lnc,$ref)的返回值赋值给 $is_contain。
    
29.             $is_contain = is_contain($lnc,$ref)
    
30.             and next;
    
31.        }
    
32.        unless ($is_reverse){
    
33.             $is_reverse = is_reverse($lnc,$ref)
    
34.             and next;
    
35.        }
    
36.        unless ($is_overlap){
    
37.             $is_overlap = is_overlap($lnc,$ref)
    
38.             and next;
    
39.        }
    
40.    }
    
41.    my $is_gapped = $is_contain || $is_reverse || $is_overlap ? 0:1;
    
42.    print $contain $_ if $is_contain;
    
43.    print $overlap $_ if $is_overlap;
    
44.    print $gapped  $_ if $is_gapped;
    
45.    print $reverse $_ if $is_reverse;
    
46. }
    
47. sub is_ln {
    
48.     my($num,$ref) = @_;
    
49.     return $num >= $ref->[0] && $num <=$ref->[1] ? 1:0;
    
50. }
    
51. sub is_contain {
    
52.     my ($lnc,$ref) = @_;
    
53.     return (is_ln($lnc->[0],$ref) && is_ln($lnc->[1],$ref)) ? 1:0;
    
54. }
    
55. sub is_reverse {
    
56.     my ($lnc,$ref) = @_;
    
57.     return is_contain($ref,$lnc);
    
58. }
    
59. sub is_overlap {
    
60.     my ($lnc,$ref) = @_;
    
61.     return (is_ln($lnc->[0],$ref) && ! is_ln($lnc->[1],$ref)) ||
    
62.            (is_ln($lnc->[1],$ref) && ! is_ln($lnc->[0],$ref)) ? 1:0;
    
63. }
    
64. sub is_gapped {
    
65.     my ($lnc,$ref) = @_;
    
66.     return (! is_ln($lnc->[0],$ref) && ! is_ln($lnc->[1],$ref)) ? 1:0;
    
67. }
    
68. 
    

复制代码

rubyish

overlay, contain, reverse, gapped
lnc
|| ref
|| ref
...
gapped => BIG FILE

1. #!/usr/bin/perl
   
2. use 5.016;
   
3. my @in    = qw/refGene.txt lncipedia_location.txt/;
   
4. my @out   = qw/gapped contain overlap  reverse/;
   
5. my %index = qw/gap 0 con 1 ove 2  rev 3/;
   
6. my ( $r, $l ) = map { open my ($f), $_; $f } @in;
   
7. my @fh = map { open my ($f), '>', $_; $f } @out;
   
8. my ( %r, %l );
   
9. map { @_ = split /\t/, $_, 5; push @{ $r{ $_[1] } }, [ @_[ 2, 3 ], $_ ] } <$r>;
   
10. map { @_ = split; push @{ $l{ $_[1] } }, [ @_[ 2, 3 ], $_ ] } <$l>;
    
11. 
    
12. while ( my ( $k, $v ) = each %l ) {
    
13.     next unless $r{$k};
    
14.     my @L = sort { $a->[0] <=> $b->[0] } @$v;
    
15.     my @R = sort { $b->[1] <=> $a->[1] } @{ $r{$k} };
    
16.     my %RO = map { $_ => $R[$_][2] } 0 .. $#R;
    
17.     for my $l (@L) {
    
18.         my ( $A, $B, $L, %Q ) = @$l;
    
19.         my %R = %RO;
    
20.         for my $ri ( 0 .. $#R ) {
    
21.             my ( $C, $D, $R ) = @{ $R[$ri] };
    
22.             next if $B < $C;
    
23.             last if $A > $D;
    
24.             my $cor =
    
25.                 ( $C <= $A and $B <= $D ) ? 'con'
    
26.               : ( $A < $C  and $D < $B )  ? 'rev'
    
27.               :                             'ove';
    
28.             push @{ $Q{$cor} }, $R;
    
29.             delete $R{$ri};
    
30.         }
    
31.         $Q{gap} = [ values %R ] if %R;
    
32.         local $, = "|| ";
    
33.         print { $fh[ $index{$_} ] } $L, @{ $Q{$_} } for keys %Q;
    
34.     }
    
35. }

复制代码

overlay, contain, reverse
lnc
|| ref
|| ref
...
gapped
lnc

1. #!/usr/bin/perl
   
2. use 5.016;
   
3. my @in    = qw/refGene.txt lncipedia_location.txt/;
   
4. my @out   = qw/gapped contain overlap  reverse/;
   
5. my %index = qw/gap 0 con 1 ove 2  rev 3/;
   
6. my ( $r, $l ) = map { open my ($f), $_; $f } @in;
   
7. my @fh = map { open my ($f), '>', $_; $f } @out;
   
8. my ( %r, %l );
   
9. map { @_ = split /\t/, $_, 5; push @{ $r{ $_[1] } }, [ @_[ 2, 3 ], $_ ] } <$r>;
   
10. map { @_ = split; push @{ $l{ $_[1] } }, [ @_[ 2, 3 ], $_ ] } <$l>;
    
11. 
    
12. while ( my ( $k, $v ) = each %l ) {
    
13.     next unless $r{$k};
    
14.     my @L = sort { $a->[0] <=> $b->[0] } @$v;
    
15.     my @R = sort { $b->[1] <=> $a->[1] } @{ $r{$k} };
    
16.     for my $l (@L) {
    
17.         my ( $A, $B, $L, %Q ) = @$l;
    
18.         for my $r ( @R ) {
    
19.             my ( $C, $D, $R ) = @$r;
    
20.             next if $B < $C;
    
21.             last if $A > $D;
    
22.             my $cor =
    
23.                 ( $C <= $A and $B <= $D ) ? 'con'
    
24.               : ( $A < $C  and $D < $B )  ? 'rev'
    
25.               :                             'ove';
    
26.             push @{ $Q{$cor} }, $R;
    
27.         }
    
28.         my $all = 0;
    
29.         $all += @$_ for values %Q;
    
30.         $Q{gap} = [] if @R > $all;
    
31.         local $, = "|| ";
    
32.         print { $fh[ $index{$_} ] } $L, @{ $Q{$_} } for keys %Q;
    
33.     }
    
34. }

复制代码

yanglc2013

回复 13# rubyish
老大还有gapped的呢？


   

rubyish

yanglc2013 发表于 2013-04-04 11:41
回复 11# rubyish
行啊。谢谢老大啊

BU 行啊。
有 ERROR!!
输出有少了!

rubyish

看懂 le

1. #!/usr/bin/perl
   
2. use 5.016;
   
3. my @in    = qw/refGene.txt lncipedia_location.txt/;
   
4. my @out   = qw/gapped contain overlap  reverse/;
   
5. my %index = qw/gap 0 con 1 ove 2  rev 3/;
   
6. my ( $r, $l ) = map { open my ($f), $_; $f } @in;
   
7. my @fh = map { open my ($f), '>', $_; $f } @out;
   
8. my ( %r, %l );
   
9. map { @_ = split /\t/, $_, 5; push @{ $r{ $_[1] } }, [ @_[ 2, 3 ], $_ ] } <$r>;
   
10. map { @_ = split; push @{ $l{ $_[1] } }, [ @_[ 2, 3 ], $_ ] } <$l>;
    
11. 
    
12. while ( my ( $k, $v ) = each %l ) {
    
13.     next unless $r{$k};
    
14.     my @L = sort { $a->[0] <=> $b->[0] } @$v;
    
15.     my @R = sort { $b->[1] <=> $a->[1] } @{ $r{$k} };
    
16.     for my $l (@L) {
    
17.         my ( $A, $B, $L, %Q ) = @$l;
    
18.         for my $r ( @R ) {
    
19.             my ( $C, $D, $R ) = @$r;
    
20.             next if $B < $C;
    
21.             last if $A > $D;
    
22.             my $cor =
    
23.                 ( $C <= $A and $B <= $D ) ? 'con'
    
24.               : ( $A < $C  and $D < $B )  ? 'rev'
    
25.               :                             'ove';
    
26.             $cor eq 'ove' ? push @{ $Q{$cor} }, $R : $Q{$cor}++;
    
27.         }
    
28.         my $all = 0;
    
29.         $all += ref $_ eq 'ARRAY' ? @$_ : $_ for values %Q;
    
30.         $Q{gap}++ if @R > $all;
    
31.         for ( keys %Q ) {
    
32.             print { $fh[ $index{$_} ] } $L and next if $_ ne 'ove';
    
33.             my @p = ( $L, @{ $Q{$_} } );
    
34.             s/\r?\n$// for @p;
    
35.             say { $fh[ $index{$_} ] } "@p";
    
36.         }
    
37.     }
    
38. }

复制代码

yanglc2013

回复 11# rubyish
行啊。谢谢老大啊

rubyish

3Q~
是否适合:

1. 
   
2. ERROR!!
   
3. delete
   
4. 
   

复制代码

yanglc2013

回复 8# rubyish

就是lnc中的其中一个染色体，比如chr1，但是ref中却存在很多chr1，他们的区间都不一样。


   

yanglc2013

回复 8# rubyish

因为在refGene中有比较多重复的染色体。比如有两个chr1，那么lnc中的每个chr1就必须和refGene中的两个chr1都比较。出现的情况就可能是gapped，或者和另一个出现overlap。