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

yanglc2013

其中ref图中的数据为库，在同一条染色体上（两个图片中蓝色区域必须先相同），如果lnc图中红色矩形内的区域包含在ref图中红色区域内输出文件1（contain），

两个红色矩形框内的数据有交集输出文件2（overlap），无交集输出文件3（gapped）.如果lnc区间比ref大，就是ref图中红色区域包含在lnc图中红色区域输出文件4（reserve）.

希望输出的格式：（比如下面是overlap，不仅要输出lnc的数据，还有在同一行加上ref对应的数据。）

lnc-C1orf227-1:2        chr1        213028819        213031452        -    NM_014053        chr1        213031596        213072705        +        2210

后面的高手只是输出了lnc的数据。

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;

复制代码

yanglc2013

回复 2# picbhan

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

难题.rar (1.89 MB, 下载次数: 9)

2013-04-08 19:16 上传

点击文件名下载附件

   

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.             }

复制代码

yanglc2013

回复 51# picbhan

大哥，还有个问题请教你。
如果我想在输出的结果中的结果中，同时加上对应的ref中的数据，比如。
lnc-C1orf227-1:2        chr1        213028819        213031452        -    NM_014053        chr1        213031596        213072705        +        2210

应该怎么修改print。

现在输出的结果形式如下：
lnc-C1orf227-1:2        chr1        213028819        213031452        -


   

dahe_1984

rubyish 发表于 2013-03-31 07:41

你写的code怎么没有了？ 我还想学习下呢

yanglc2013

回复 51# picbhan

不会出现这样情况的。因为ref中每个染色体，比如很多chr1都是来源于有一个DNA，这样发生的剪切，如果存在contain，那么另一种情况只能是gapped。


   

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.         }

复制代码

rubyish

{:3_188:}

1. 
   
2. Delete
   
3. 有错误
   

复制代码

rubyish

{:3_188:}

1. 
   
2. Delete
   
3. 有错误
   

复制代码

yanglc2013

回复 47# rubyish
不是啊，输出的是lnc中的数据。（比如lnc-HES5-1:2        chr3        14 100        -）中的 [14,100]这个区间和ref中的（NM_016505        chr3        50 128        +        827
）中的[50,128]这个区间比较，它们有一个相交的区间[50,100]，所以在overlap文件中输出。如果和（NR_026818        chr3        110 150        -        585）比较，即区间[110,150]，没有公共部分，就输出到文件gapped。如果和（NR_026835        chr3        12 150        -        586）比较，lnc包含在ref区间中，输出到contain中，如果反过来，ref包含在lnc中，就输出到reserve文件中，输出的文件都是（lnc-HES5-1:2        chr3        14 100        -），如果四种情况都存在，那么，输出的四个文件中，都应该有（lnc-HES5-1:2        chr3        14 100        -）