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

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中。
你再好好检查下吧。

picbhan

回复 45# yanglc2013


    你觉得这样子输出有具体意义吗？inc中任何一个region都会被输出到non_include中去，只要你的refGene足够多。对于后面三种有交集的情况，如果要输出到多个文件中，你把下面的代码替换一下就行。

1. my %result; # one of ref_include_inc, inc_include_ref, overlap or undef
   
2.             while (@ref_data) {
   
3.                 my @rd = @{+shift @ref_data};
   
4.                 # igore genes on the left side of inc region
   
5.                 next if $rd[3] < $id[2];
   
6.                 push @store, [@rd];
   
7.                 # stop when this gene is on the right side of inc region
   
8.                 last if $rd[2] > $id[3];
   
9.                 # this gene must overlap with this inc region
   
10.                 if ($id[2] >= $rd[2] && $id[3] <= $rd[3]) {
    
11.                     $result{ref_include_inc} = 1;
    
12.                 }
    
13.                 elsif ($rd[2] >= $id[2] && $rd[3] <= $id[3]) {
    
14.                     $result{inc_include_ref} = 1;
    
15.                 }
    
16.                 else {
    
17.                     $result{overlap} = 1;
    
18.                 }
    
19.             }
    
20.             # output result
    
21.             if (%result) {
    
22.                 if ($result{ref_include_inc}) {
    
23.                     print {$ref_include_inc} join("\t", @id), "\n";
    
24.                 }
    
25.                 if ($result{inc_include_ref}) {
    
26.                     print {$inc_include_ref} join("\t", @id), "\n";
    
27.                 }
    
28.                 if ($result{overlap}) {
    
29.                     print {$overlap} join("\t", @id), "\n";
    
30.                 }
    
31.             }
    
32.             else { # if $result is undef, means no overlapped genes
    
33.                 print {$non_include} join("\t", @id), "\n";
    
34.             }
    
35.             
    
36.             # put this genes back into @ref_data, use to check next inc region
    
37.             unshift @ref_data, @store;
    
38.             @store = ();
    
39.         }

复制代码

picbhan

回复 54# yanglc2013


    你在下面的代码中加入输出的信息（包括@rd和@id），并且把原来的三个赋值语句改为一个表示已经输出过的flag，这样在后面判断是否是non-overlap的时候就只需要看这个flag是不是真值了。

1. 
   
2.                 if ($id[2] >= $rd[2] && $id[3] <= $rd[3]) {
   
3.                     $result{ref_include_inc} = 1;
   
4.                 }
   
5.                 elsif ($rd[2] >= $id[2] && $rd[3] <= $id[3]) {
   
6.                     $result{inc_include_ref} = 1;
   
7.                 }
   
8.                 else {
   
9.                     $result{overlap} = 1;
   
10.                 }
    
11.                 if (%result) {
    
12.                 if ($result{ref_include_inc}) {
    
13.                     print {$ref_include_inc} join("\t", @id), "\n";
    
14.                 }
    
15.                 if ($result{inc_include_ref}) {
    
16.                     print {$inc_include_ref} join("\t", @id), "\n";
    
17.                 }
    
18.                 if ($result{overlap}) {
    
19.                     print {$overlap} join("\t", @id), "\n";
    
20.                 }
    
21.             }
    
22.             else { # if $result is undef, means no overlapped genes
    
23.                 print {$non_include} join("\t", @id), "\n";
    
24.             }

复制代码

picbhan

yanglc2013 发表于 2013-04-08 19:17
回复 2# picbhan

谢谢啊，你的程序非常好用。还有个难题想请教你，我也写了几个，一塌糊涂啊。

首先，我没有看到你说的你也写了几个，我只看到你的问题。
其次，其实你的问题很简单，如果你真认真学了Perl并且自己动手写了，肯定是可以完成的，这根本不是什么难题，仅仅只是简单的格式转换而已。
最后，没有人可以一直帮你写代码，而且这种问题都等着别人给你写代码的话只会降低你自己的效率。好好学perl基础编程吧，足以在生物里用了。

1. #!perl
   
2. use List::MoreUtils 'pairwise';
   
3. 
   
4. $fin = 'refGene.txt';
   
5. open my $in, '<', $fin or die "Can't open $fin for reading. $!\n";
   
6. $fout = 'result.txt';
   
7. open my $to, '>', $fout or die "Can't create $fout for writing. $!\n";
   
8. 
   
9. while (<$in>) {
   
10.     chomp;
    
11.     my @t = split;
    
12.     @start = split /,/, $t[1];
    
13.     @end = split /,/, $t[2];
    
14.     pairwise { print {$to} "$t[0]\t$a\t$b\n" } @start, @end;
    
15. }
    
16. close $in;
    
17. close $to;

复制代码