2012-09-20 66 views
0

我想將我的一個Perl腳本轉換爲R腳本。我在R A數據幀,它看起來像(忽略的列名) -如果列相同,合併R dataframe中的連續行

CHR  START   END  TYPE 
chr1 945493   945593 normal 
chr1 945593  947374 normal 
chr1 947374  947474 normal 
chr1 947474  947574 gain 
chr1 947574  947674 gain 
chr1 947674  960364 gain 
chr1 960364  960464 normal 
chr22 17290491 17290591 normal 
chr22 17290591 17290691 normal 
chr22 17290691 17290791 gain 
chr22 17290791 17292513 gain 
chr22 17292513 17292613 gain 
chr22 17292613 17292713 gain 
chr22 17292713 17293046 gain 
chr22 17293346 17298475 gain 
chr22 17298475 17298575 gain 
chr22 17298575 17298675 normal 
chr22 17298675 17303632 normal 
chr22 17303632 17303732 loss 
chr22 17303732 17303832 normal 
chrX 154162621 154181221 normal 
chrX 154181221 154181321 normal 
chrX 154181321 154181421 loss 
chrX 154181421 154181521 loss 
chrX 154181521 154181621 loss 
chrX 154181621 154181721 loss 
chrX 154181721 154216867 loss 
chrX 154216867 154216967 normal 
chrX 154216967 154217067 normal 
chrX 154217067 154217167 normal

如果至少連續5行具有在「CHR」列和「類型」列相同的值,然後將所有這些行中一個行,以便START列應該具有第一行的值,END列具有最後一行的值,最後只返回具有「增益」或「損失」類型的行。因此,所需的輸出是:

chr22 17290691  17298575  gain 
chrX 154181321  154216867  loss

我在做什麼現在的問題是:

  1. 保存與 「write.table」 數據幀。
  2. 使用該perl腳本:

    open $first, "<",$ARGV[0] or die "Unable to open input file: $!"; 
        my $count=1; 
        $_ = <$first>; 
        chomp; 
        my ($p_key, $p_col1, $p_col2,$p_cnv) = split; 
    
        while(<$first>) { 
         chomp; 
         my ($key, $col1, $col2,$cnv) = split; 
         if ($key eq $p_key and $cnv eq $p_cnv) { 
         $p_col2 = $col2; 
         $count++; 
    
         } elsif ($count > 4){ 
    
    
         print $p_key,"\t", $p_col1,"\t", $p_col2,"\t", $p_cnv,"\n" if($p_cnv eq "gain" or $p_cnv eq "loss"); 
         ($p_key, $p_col1, $p_col2, $p_cnv) = ($key, $col1, $col2, $cnv); 
         $count=1; 
         } 
    
         else { 
    
        ($p_key, $p_col1, $p_col2, $p_cnv) = ($key, $col1, $col2, $cnv); 
         $count=1; 
         } 
    } 
    

我認爲這是一個額外的步驟先保存數據框,然後使用Perl腳本。任何人都可以請建議一個更簡單的方法來做到這一點在R - 任何包或任何其他技巧?

+0

作爲第一步,爲什麼不簡單地將Perl腳本翻譯成R(看起來相當簡單直接的代碼),並查看是否遇到任何問題? – joran

+0

這就是我想要做的。在Perl中看起來很容易,但我無法弄清楚我該如何在R中執行它。 – Vikas

+0

在你的perl腳本中,你使用'&'(按位AND)而不是'&&'(C形邏輯AND)。我會推薦使用'或'和'和',因爲它們更透明並且優先級更低。 – TLP

回答

1

我擔心如果在一個染色體中存在TYPE的插入替代值,您應該想要中斷序列(即將它們視爲不同)。你沒有特別說明,但我認爲生物學會保證有額外的要求。因此需要創建另一個變量。如果沒有相反的建議,我們將假定數據框名爲cdat。這看起來在TYPE的連續運行中,應用測試,並在開始處綁定CHR和START,並在最後一個元素的END和TYPE處綁定。

cdat$conseq <-cumsum(c(1, cdat$TYPE[-1] != cdat$TYPE[-length(cdat$TYPE)])) 
do.call(rbind, 
    by(cdat, list(cdat$CHR, cdat$conseq), 
     function(df) 
      if(NROW(df) >=5 & df$TYPE[1] %in% c("gain", "loss")) { 
       cbind(df[1, c("CHR", "START")] , df[NROW(df), c("END", "TYPE")]) 
       } else{NULL})) 
    CHR  START  END TYPE 
10 chr22 17290691 17298575 gain 
23 chrX 154181321 154216867 loss 

的conseq矢量由下式值與事先值和cumsum建立() - 荷蘭國際集團沿其整個長度的新值的外觀。由於這些變量是一個更短的元素。在開始時將1添加爲佔位符,以使其與數據框對齊。

+0

完美。非常感謝。你能告訴我一點關於第一行嗎?我認爲(正如你所說)你正在創造新的變量。 – Vikas

+0

它總結了TYPE值的變化。與已經移動一個位置的向量子集的比較或計算是R中相當常見的數據操作策略。參見上文。 –

+0

AAh!我想我明白了。使用第一行代碼,如果連續,則創建此變量以將相同的數字分配給相同的TYPE。 – Vikas

2

沿着這些線?

library(plyr) 
ddply(x[x$TYPE %in% c("gain", "loss"), ], 
     .(CHR, TYPE), 
     function(z){if(nrow(z) < 5) NULL else z[range(seq_len(nrow(z))), ]} 
    ) 

    CHR  START  END TYPE 
1 chr22 17290691 17290791 gain 
2 chr22 17298475 17298575 gain 
3 chrX 154181321 154181421 loss 
4 chrX 154181721 154216867 loss 
+0

+1爲此,在行號上使用'range'非常光滑。除了看起來他只想要一行。我還想知道,被其他值中斷的「收益」或「損失」是否應該被視爲同一部分的一部分。目前你的方法會將所有這些結合在一起。 –

+0

@DWin謝謝。而且,是的,說實話,我沒有試圖完全理解OP的要求,也沒有想出完美的結果。但我明白你的意思了。我可能會在幾分鐘內編輯這個答案。 – Andrie

0

如果您熟悉SQL,想要做了很多數據幀的處理的,另一種選擇是sqldf庫,它可以讓你R上dataframes執行SQL查詢。它使得這樣的操作非常簡單。

還有R-Perl interface,它可以讓你保留現有的Perl代碼,然後用結果做R事情。