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

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 $include , '>', 'include.txt'
    
26.     or die "Can't create 'include.txt' for writing. $!\n";
    
27. open my $non_include , '>', 'non_include.txt'
    
28.     or die "Can't create 'non_include.txt' for writing. $!\n";
    
29. open my $overlap , '>', 'overlap.txt'
    
30.     or die "Can't create 'overlap.txt' for writing. $!\n";
    
31. 
    
32. 
    
33. while (my ($chr, $v) = each %ref) {
    
34.     my @ref_data = sort {$a->[2] <=> $b->[2]} @$v;
    
35.     next unless exists $data{$chr};
    
36.     my @inc_data = sort {$a->[2] <=> $b->[2]} @{$data{$chr}};
    
37.     for (@ref_data) {
    
38.         my @rd = @$_;
    
39.         last unless @inc_data;
    
40.         while (1) {
    
41.             last unless @inc_data;
    
42.             my @id = @{+shift @inc_data};
    
43.             if ($id[3] < $rd[2]) {
    
44.                 print {$non_include} join("\t", @id), "\n";
    
45.                 next;
    
46.             }
    
47.             if ($id[2] > $rd[3]) {
    
48.                 unshift @inc_data, [@id];
    
49.                 last;
    
50.             }
    
51.             if ($id[2] >= $rd[2] && $id[3] <= $rd[3]) {
    
52.                 print {$include} join("\t", @id), "\n";
    
53.             }
    
54.             else {
    
55.                 print {$overlap} join("\t", @id), "\n";
    
56.             }
    
57.         }
    
58.     }
    
59.     delete $data{$chr};
    
60. }
    
61. 
    
62. while (my ($chr, $v) = each %data) {
    
63.     my @inc_data = @$v;
    
64.     for (@inc_data) {
    
65.         print {$non_include} join("\t", @$_), "\n";
    
66.     }
    
67. }
    
68. 
    
69. close $include;
    
70. close $non_include;
    
71. close $overlap;

复制代码

Perlvim

是测试DNA相似度吗？

1. refGene 文件中的
NR_026818        chr1        34610        36081        -        585
NR_033908        chr1        700244        714068        -        590
6个字段，是否只使用第2个到第4个，其他没有用？

2. refGene 文件中 第二个字段相同的，是否要合并范围？计算并集？

3. lnc 中的数据，计算在哪个范围，是否有交集，是在 ref 那个对应字段全部范围内计算吗？

mcshell

不懂生物 真苦逼{:3_186:}

yanglc2013

回复 2# Perlvim
一个染色体可以分好几段，虽然你说的有重复，但是其起始位置和终止位置不一样，所以不要合并啊。这是求gene的location信息。判断lnc中的染色体（chr）其起始位置和终止位置是不是包含在ref中的起始位置和终止位置中。就是你指的第三列（起始位置）和第四列（终止位置）。

Perlvim

回复 5# yanglc2013


    意思是一个染色体中一个位置有很多片段，要每个片段都查找，是否有交集吗？

yanglc2013

回复 6# Perlvim
对啊。就是比对染色体，如果lnc中的那一段染色体在ref中就输出，不在输出待另一个文件，有一部分在又输出一个文件。

dahe_1984

数据好大，机器都快跑死了，你这个够牛的，不好做啊。

blackfur

回复 1# yanglc2013


    神马叫‘落在’，是两个数值相同的意思吗，生物不大懂，建议楼主先扫盲

dahe_1984

blackfur 发表于 2013-03-28 00:53
回复 1# yanglc2013

落在就是包含呗。 比如 4..6 落在 1..10中