0

我正在尝试将我的 Perl 脚本之一转换为 R 脚本。我在 R 中有一个数据框，看起来像（忽略列名）-

CHR 起始端类型
chr1 945493 945593 正常
chr1 945593 947374 正常
chr1 947374 947474 正常
chr1 947474 947574 增益
chr1 947574 947674 增益
chr1 947674 960364 增益
chr1 960364 960464 正常
chr22 17290491 17290591 正常
chr22 17290591 17290691 正常
chr22 17290691 17290791 增益
chr22 17290791 17292513 增益
chr22 17292513 17292613 增益
chr22 17292613 17292713 增益
chr22 17292713 17293046 增益
chr22 17293346 17298475 增益
chr22 17298475 17298575 增益
chr22 17298575 17298675 正常
chr22 17298675 17303632 正常
chr22 17303632 17303732 损失
chr22 17303732 17303832 正常
chrX 154162621 154181221 正常
chrX 154181221 154181321 正常
chrX 154181321 154181421 损失
chrX 154181421 154181521 损失
chrX 154181521 154181621 损失
chrX 154181621 154181721 损失
chrX 154181721 154216867 损失
chrX 154216867 154216967 正常
chrX 154216967 154217067 正常
chrX 154217067 154217167 正常

如果至少 5 个连续行在“CHR”列和“TYPE”列中具有相同的值，则将所有这些行合并为一行，以便 START 列应具有第一行的值，END 列应具有最后一行的值，最后只返回具有“增益”或“损失”类型的行。所以想要的输出是：

chr22 17290691 17298575 增益
chrX 154181321 154216867 损失

我现在正在做的是：

  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 中执行此操作的更简单方法 - 任何包或任何其他技巧？

4

3 回答 3

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
于 2012-09-20T13:16:31.567 回答
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 向量是通过将下一个 TYPE 值与其先前值进行比较并 cumsum() 建立新值沿其全长出现的。因为这些变量要短一个元素。1 在开头添加为占位符，以使其与数据框对齐。

于 2012-09-20T14:40:46.210 回答
0

如果您熟悉 SQL 并希望对数据帧进行大量操作，则另一个选择是sqldf库，它允许您对 R 数据帧执行 SQL 查询。它使这样的操作非常容易。

还有R-Perl 接口，它允许你保留现有的 Perl 代码，然后用结果做 R 事情。

于 2012-09-20T13:21:55.420 回答