2013-09-21 31 views
1

我想应用一个函数返回一个矩阵大型data.table对象的每一行(原始文件大约30 GB ,我有80 GB RAM),然后找回一个data.table对象。我想要有效地做到这一点。我目前的做法是:应用函数返回一个data.table,或直接将列表转换为data.table

my.function <- function(x){ 
    alnRanges<-cigarToIRanges(x[6]); 
    alnStarts<-start(alnRanges)+as.numeric(x[4])-1; 
    alnEnds<-end(alnRanges)+as.numeric(x[4])-1; 
    y<-x[-4]; 
    ys<-matrix(rep(y,length(alnRanges)),nrow=length(alnRanges),ncol=length(y),byrow=TRUE); 
    ys<-cbind(ys,alnStarts,alnEnds); 
    return(ys);  # ys is a matrix 
} 

    my.dt<-fread(my.file.name); 
    my.list.of.matrices<-apply(my.dt,1,my.function); 
    new.df<-do.call(rbind.data.frame,my.list.of.matrices); 
    colnames(new.df)[1:14]<-colnames(my.dt)[-4]; 
    new.dt<-as.data.table(new.df); 

注1:我指定的my.function只是为了表明它返回一个矩阵,而我的申请,因此线是矩阵的列表。注2:我不知道我正在做的操作有多慢,但似乎我可以减少行数。例如,将数据框转换为大型对象的数据表是否慢?

重复的例子:

注意,阿伦和罗兰让我想到这个问题,所以我仍然在它的工作...可能是我不需要这些行难...

我想要取得一个sam文件,然后创建一个新的坐标文件,每个读取根据其CIGAR字段进行拆分。

My sam file: 
qname rname pos cigar 
2218 chr1 24613476 42M2S 
2067 chr1 87221030 44M 
2129 chr1 79702717 44M 
2165 chr1 43113438 44M 
2086 chr1 52155089 4M921N40M 

code: 

library("data.table"); 
library("GenomicRanges"); 

sam2bed <- function(x){ 
    alnRanges<-cigarToIRanges(x[4]); 
    alnStarts<-start(alnRanges)+as.numeric(x[3])-1; 
    alnEnds<-end(alnRanges)+as.numeric(x[3])-1; 
    #y<-as.data.frame(x[,pos:=NULL]); 
    #ys<-y[rep(seq_len(nrow(y)),length(alnRanges)),]; 
    y<-x[-3]; 
    ys<-matrix(rep(y,length(alnRanges)),nrow=length(alnRanges),ncol=length(y),byrow=TRUE); 
    ys<-cbind(ys,alnStarts,alnEnds); 
    return(ys); 
} 


sam.chr.dt<-fread(sam.parent.chr.file); 
setnames(sam.chr.dt,old=c("V1","V2","V3","V4"),new=c("qname","rname","pos","cigar")); 
bed.chr.lom<-apply(sam.chr.dt,1,sam2bed); 
> bed.chr.lom 
[[1]] 
          alnStarts alnEnds 
[1,] "2218" "chr1" "42M2S" "24613476" "24613517" 

[[2]] 
         alnStarts alnEnds 
[1,] "2067" "chr1" "44M" "87221030" "87221073" 

[[3]] 
         alnStarts alnEnds 
[1,] "2129" "chr1" "44M" "79702717" "79702760" 

[[4]] 
         alnStarts alnEnds 
[1,] "2165" "chr1" "44M" "43113438" "43113481" 

[[5]] 
           alnStarts alnEnds 
[1,] "2086" "chr1" "4M921N40M" "52155089" "52155092" 
[2,] "2086" "chr1" "4M921N40M" "52156014" "52156053" 

bed.chr.df<-do.call(rbind.data.frame,bed.chr.lom); 

> bed.chr.df 
    V1 V2  V3 alnStarts alnEnds 
1 2218 chr1  42M2S 24613476 24613517 
2 2067 chr1  44M 87221030 87221073 
3 2129 chr1  44M 79702717 79702760 
4 2165 chr1  44M 43113438 43113481 
5 2086 chr1 4M921N40M 52155089 52155092 
6 2086 chr1 4M921N40M 52156014 52156053 

bed.chr.dt<-as.data.table(bed.chr.df); 

> bed.chr.dt 
    V1 V2  V3 alnStarts alnEnds 
1: 2218 chr1  42M2S 24613476 24613517 
2: 2067 chr1  44M 87221030 87221073 
3: 2129 chr1  44M 79702717 79702760 
4: 2165 chr1  44M 43113438 43113481 
5: 2086 chr1 4M921N40M 52155089 52155092 
6: 2086 chr1 4M921N40M 52156014 52156053 
+0

你不应该需要(并使用)'apply' data.table。目标应该是避免副本。请提供可复制的数据并解释您实际尝试做的事情。 – Roland

+0

@Roland,我需要将一个函数应用于data.table的每一行。如果我不应用apply,那么我不知道我还能做什么,我对数据表不熟悉,但我知道它们比数据框要快得多。试着制作一个例子...... – Dnaiel

+1

正如我所说的,提供一个工作示例(包括数据和必要的软件包)并告诉我们你想实现什么。如果你将它与'apply'混合使用,你将无法获得data.table的优势。 – Roland

回答

3

假设ff是你data.table,这个怎么样?

splits <- cigarToIRangesListByAlignment(ff$cigar, ff$pos, reduce.ranges = TRUE) 
widths <- width(attr(splits, 'partitioning')) 
cbind(data.table(qname=rep.int(ff$qname, widths), 
      rname=rep.int(ff$rname, widths)), as.data.frame(splits)) 

    qname rname space start  end width 
1: 2218 chr1  1 24613476 24613517 42 
2: 2067 chr1  2 87221030 87221073 44 
3: 2129 chr1  3 79702717 79702760 44 
4: 2165 chr1  4 43113438 43113481 44 
5: 2086 chr1  5 52155089 52155092  4 
6: 2086 chr1  5 52156014 52156053 40 
+0

非常感谢!很棒的建议。我也在研究rsamtools和iranges,因为也许我可以将所有的代码改为这种方法......谢谢! – Dnaiel

相关问题