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

rubyish

ref

NM_000721        chr1        100        110        +        246
NM_020394        chr1        110        300        -        2470

lnc

lnc-FAM20C-5:2        chr1        100        150        +

輸出文件1
OR
輸出文件2

yanglc2013

回復(fù) 43# rubyish

lnc和ref比較，ref作為數(shù)據(jù)庫。輸出四個文件，輸出的內(nèi)容是lnc中的序列：

文件一：contain：就是ref圖中的紅色區(qū)間包含lnc圖中紅色區(qū)間。
文件二：overlap：就是lnc圖中的紅色區(qū)間和ref圖中紅色區(qū)間相交。
文件三：gapped：就是lnc圖中的紅色區(qū)間和ref圖中紅色區(qū)間包含，也不相交，是獨立分開的。
文件四：reverse：就是lnc圖中的紅色區(qū)間包含ref圖中的紅色區(qū)間，和文件一剛好相反。

   

yanglc2013

回復(fù) 42# picbhan

大哥，其實是這樣的，lnc中每個區(qū)間都要單獨和ref中的區(qū)間比較，每一個染色體（比如lnc-HES5-1:2        chr3        14 100        -），它和NM_016505        chr3        50 128        +        827是overlap的，但是它和NR_026818        chr3        110 150        -        585卻是non_include，所以輸出的non_include文件中應(yīng)該也包含lnc-HES5-1:2        chr3        14 100        - 啊。（就是lnc中同一個區(qū)間，他可能有很多種情況，都要分別輸出。這也沒有考慮正負(fù)鏈啊）

yanglc2013

回復(fù) 42# picbhan
我意思是對每種情況都輸出啊。就是一個region，它可能是contain，可能是gapped，可能是overlap也可能lnc包含ref，如果這四種情況都存在，那么輸出的四個文件中，都要包含這個region。大哥，真抱歉啊，我沒表達(dá)清楚。你那個程序，我還沒研究透，不是很會改你的程序。還是會出錯。


   

rubyish

比如
lnc-HES5-1:2        chr3        14 100        -

NM_016505        chr3        50 128        +        827
NR_026818        chr3        110 150        -        585
NR_026835        chr3        12 150        -        586
NR_026822        chr3        16 50        -        534
NR_0268200       chr3        6 12        -        531

輸出文件?

yanglc2013

回復(fù) 47# rubyish
不是啊，輸出的是lnc中的數(shù)據(jù)。（比如lnc-HES5-1:2        chr3        14 100        -）中的 [14,100]這個區(qū)間和ref中的（NM_016505        chr3        50 128        +        827
）中的[50,128]這個區(qū)間比較，它們有一個相交的區(qū)間[50,100]，所以在overlap文件中輸出。如果和（NR_026818        chr3        110 150        -        585）比較，即區(qū)間[110,150]，沒有公共部分，就輸出到文件gapped。如果和（NR_026835        chr3        12 150        -        586）比較，lnc包含在ref區(qū)間中，輸出到contain中，如果反過來，ref包含在lnc中，就輸出到reserve文件中，輸出的文件都是（lnc-HES5-1:2        chr3        14 100        -），如果四種情況都存在，那么，輸出的四個文件中，都應(yīng)該有（lnc-HES5-1:2        chr3        14 100        -）
   

rubyish

{:3_188:}

1. 
   
2. Delete
   
3. 有錯誤
   

復(fù)制代碼

rubyish

{:3_188:}

1. 
   
2. Delete
   
3. 有錯誤
   

復(fù)制代碼

picbhan

回復(fù) 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.         }

復(fù)制代碼

yanglc2013

回復(fù) 51# picbhan

不會出現(xiàn)這樣情況的。因為ref中每個染色體，比如很多chr1都是來源于有一個DNA，這樣發(fā)生的剪切，如果存在contain，那么另一種情況只能是gapped。