20積分懸賞，求個perl程序，新手，寫了好幾個都無果。

yanglc2013

其中ref圖中的數(shù)據(jù)為庫，在同一條染色體上（兩個圖片中藍色區(qū)域必須先相同），如果lnc圖中紅色矩形內(nèi)的區(qū)域包含在ref圖中紅色區(qū)域內(nèi)輸出文件1（contain），

兩個紅色矩形框內(nèi)的數(shù)據(jù)有交集輸出文件2（overlap），無交集輸出文件3（gapped）.如果lnc區(qū)間比ref大，就是ref圖中紅色區(qū)域包含在lnc圖中紅色區(qū)域輸出文件4（reserve）.

希望輸出的格式：（比如下面是overlap，不僅要輸出lnc的數(shù)據(jù)，還有在同一行加上ref對應的數(shù)據(jù)。）

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

后面的高手只是輸出了lnc的數(shù)據(jù)。

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 中的數(shù)據(jù)，計算在哪個范圍，是否有交集，是在 ref 那個對應字段全部范圍內(nèi)計算嗎？

mcshell

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

yanglc2013

回復 2# Perlvim
一個染色體可以分好幾段，雖然你說的有重復，但是其起始位置和終止位置不一樣，所以不要合并啊。這是求gene的location信息。判斷l(xiāng)nc中的染色體（chr）其起始位置和終止位置是不是包含在ref中的起始位置和終止位置中。就是你指的第三列（起始位置）和第四列（終止位置）。

Perlvim

回復 5# yanglc2013


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

yanglc2013

回復 6# Perlvim
對啊。就是比對染色體，如果lnc中的那一段染色體在ref中就輸出，不在輸出待另一個文件，有一部分在又輸出一個文件。

dahe_1984

數(shù)據(jù)好大，機器都快跑死了，你這個夠牛的，不好做啊。

blackfur

回復 1# yanglc2013


    神馬叫‘落在’，是兩個數(shù)值相同的意思嗎，生物不大懂，建議樓主先掃盲

dahe_1984

blackfur 發(fā)表于 2013-03-28 00:53
回復 1# yanglc2013

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