20积分悬赏，求个perl程序，新手，写了好几个都无果。

yanglc2013

回复 19# rubyish
能加上注释吗

yanglc2013

回复 32# picbhan
不考虑正负啊。我只是想学学，应该怎么改程序。加上注释，我就可以看明白了。有些地方我不太懂。


   

picbhan

1. #!/usr/bin/perl
   
2. use strict;
   
3. use warnings;
   
4. 
   
5. my (%ref, %data);
   
6. 
   
7. open my $in , '<', 'refGene.txt'
   
8.     or die "Can't open 'refGene.txt' for reading. $!\n";
   
9. while (<$in>) {
   
10.     chomp;
    
11.     my @t = split;
    
12.     push @{$ref{$t[1]}}, [@t];
    
13. }
    
14. close $in;
    
15. 
    
16. open $in , '<', 'lncipedia_location.txt'
    
17.     or die "Can't open 'lncipedia_location.txt' for reading. $!\n";
    
18. while (<$in>) {
    
19.     chomp;
    
20.     my @t = split;
    
21.     push @{$data{$t[1]}}, [@t];
    
22. }
    
23. close $in;
    
24. 
    
25. open my $ref_include_inc , '>', 'ref_include_inc.txt'
    
26.     or die "Can't create 'ref_include_inc.txt' for writing. $!\n";
    
27. open my $inc_include_ref , '>', 'inc_include_ref.txt'
    
28.     or die "Can't create 'inc_include_ref.txt' for writing. $!\n";
    
29. open my $non_include , '>', 'non_include.txt'
    
30.     or die "Can't create 'non_include.txt' for writing. $!\n";
    
31. open my $overlap , '>', 'overlap.txt'
    
32.     or die "Can't create 'overlap.txt' for writing. $!\n";
    
33. open my $to , '>', 'ref_sorted.txt' or die "Can't create  for writing. $!\n";
    
34. 
    
35. 
    
36. while (my ($chr, $v) = each %ref) {
    
37.     # sort region using start sites, if equal, using end sites
    
38.     my @ref_data = sort {$a->[2] <=> $b->[2] or $a->[3] <=> $b->[3]} @$v;
    
39.     # output sorted refGenes
    
40.     print {$to} join("\t", @$_) . "\n" for @ref_data;
    
41.     # if no such chromosome in inc, ignore and try next chr
    
42.     next unless exists $data{$chr};
    
43.     # sort region using start sites, if equal, using end sites
    
44.     my @inc_data = sort {$a->[2] <=> $b->[2] or $a->[3] <=> $b->[3]} @{$data{$chr}};
    
45.     my @store; # use to store some ref regions, since inc regions may overlapped
    
46.     for (@inc_data) {
    
47.         my @id = @$_;
    
48.         if (@ref_data) {
    
49.             my $result; # one of ref_include_inc, inc_include_ref, overlap or undef
    
50.             while (@ref_data) {
    
51.                 my @rd = @{+shift @ref_data};
    
52.                 # igore genes on the left side of inc region
    
53.                 next if $rd[3] < $id[2];
    
54.                 # record each overlapped genes
    
55.                 push @store, [@rd];
    
56.                 # stop when this gene is on the right side of inc region
    
57.                 last if $rd[2] > $id[3];
    
58.                 # this gene must overlap with this inc region
    
59.                 if ($id[2] >= $rd[2] && $id[3] <= $rd[3]) {
    
60.                     $result = 'ref_include_inc';
    
61.                 }
    
62.                 elsif ($rd[2] >= $id[2] && $rd[3] <= $id[3]) {
    
63.                     $result = 'inc_include_ref';
    
64.                 }
    
65.                 else {
    
66.                     $result = 'overlap';
    
67.                 }
    
68.             }
    
69.             # output result
    
70.             if ($result) {
    
71.                 if ($result eq 'ref_include_inc') {
    
72.                     print {$ref_include_inc} join("\t", @id), "\n";
    
73.                 }
    
74.                 elsif ($result eq 'inc_include_ref') {
    
75.                     print {$inc_include_ref} join("\t", @id), "\n";
    
76.                 }
    
77.                 else {
    
78.                     print {$overlap} join("\t", @id), "\n";
    
79.                 }
    
80.             }
    
81.             else { # if $result is undef, means no overlapped genes
    
82.                 print {$non_include} join("\t", @id), "\n";
    
83.             }
    
84.             
    
85.             # put this genes back into @ref_data, use to check next inc region
    
86.             unshift @ref_data, @store;
    
87.             @store = ();
    
88.         }
    
89.         else {
    
90.             # no regions in ref left, print all regions left in inc to non_include
    
91.             print {$non_include} join("\t", @id), "\n";
    
92.             next;
    
93.         }
    
94.     }
    
95.     delete $data{$chr};
    
96. }
    
97. 
    
98. while (my ($chr, $v) = each %data) {
    
99.     my @inc_data = @$v;
    
100.     for (@inc_data) {
     
101.         print {$non_include} join("\t", @$_), "\n";
     
102.     }
     
103. }
     
104. 
     
105. close $to;
     
106. close $ref_include_inc;
     
107. close $inc_include_ref;
     
108. close $non_include;
     
109. close $overlap;

复制代码

yanglc2013

回复 35# picbhan

这个程序思路确实很好，希望能加上注释。这样比较容易看。

yanglc2013

回复 35# picbhan

抱歉啊，刚复制到editplus才看到有注释，非常感谢啊。


   

picbhan

回复 36# yanglc2013


    还可以更改一下，在判断是否属于包含还是overlap的时候加上last，这样可以在同一个inc的region与refgene存在多种情况时优先选择包含的情况，也可以减少运行时间。

1. if ($id[2] >= $rd[2] && $id[3] <= $rd[3]) {
   
2.                     $result = 'ref_include_inc';
   
3.                     last;
   
4.                 }
   
5.                 elsif ($rd[2] >= $id[2] && $rd[3] <= $id[3]) {
   
6.                     $result = 'inc_include_ref';
   
7.                     last;
   
8.                 }

复制代码

yanglc2013

回复 38# picbhan

大哥，输出结果中，有两种情况（non_include和overlap）少了数据，我没找到原因。请教你下。

    test for help.rar (3.17 KB, 下载次数: 3)

2013-03-29 14:48 上传

点击文件名下载附件

dahe_1984

picbhan 发表于 2013-03-28 09:28
回复 8# dahe_1984

我自己写的是把ref中ch*作为Key值，然后foreach($ref[2]..$ref[3]),$hash{$[1]}++;  

Inc中判断exists($hash{}), 然后判断in，out，还是overlap， 生成的时候就跑死了，跟您的程度无关。

picbhan

回复 40# dahe_1984


    我只能说你太暴力，比inc中每一个region对ref中整个chr*的比较还要暴力。呵呵

picbhan

回复 39# yanglc2013


    不是因为少了，而是inc中每个region我都确保只输出到一个文件中，没有输出到多个文件中，如果某个region与ref中的gene存在完全包含关系，就不会输出到overlap中。比如你说overlap中少了 lnc-TAF12-1:4        chr2        9        15        - ， 那是因为ref中 NM_001097613        chr2        11 13        -        1697完全落在了inc的区间内，所以那个region就只输出到了inc_include_ref这个文件中。你如果要想一个region输出到多个文件中（对每种情况都输出）的话，可以稍微改一下代码就行，把标量$result改为%result,然后用不同情况下的字符串作为key，最近依据key值输出到不同文件中就行。但是这种情况下不能在找到包含关系之后就立即last（就是我最后帮你加的那两个last语句）。

另外你说non_include 少了 lnc-HES5-1:2        chr3        14 100        - ，这是因为你说不需要考虑正负链，那么这个region与ref中 NM_016505        chr3        50 128        +        827 是overlap的，所以输出在overlap中。
你再好好检查下吧。