|
5
|
1 args<-commandArgs(TRUE)
|
|
|
2
|
|
|
3 input=args[1]
|
|
|
4 tmp_dir=args[2]
|
|
|
5 nbcall=as.numeric(args[3])
|
|
|
6 cellularity=as.numeric(args[4])
|
|
|
7 output=args[5]
|
|
|
8 method=args[6]
|
|
|
9 userId=args[7]
|
|
|
10 signalType=args[8]
|
|
|
11
|
|
|
12 library(MPAgenomics)
|
|
|
13 workdir=file.path(tmp_dir, "mpagenomics",userId)
|
|
|
14 setwd(workdir)
|
|
|
15
|
|
|
16 CN=read.table(input,header=TRUE)
|
|
|
17 uniqueChr=unique(CN$chromosome)
|
|
|
18 drops=c("chromosome","position","probeName")
|
|
|
19 CNsignal=CN[,!(names(CN)%in% drops),drop=FALSE]
|
|
|
20
|
|
|
21 samples=names(CNsignal)
|
|
|
22
|
|
|
23 if (signalType=="CN")
|
|
|
24 {
|
|
|
25
|
|
|
26 result=data.frame(sampleNames=character(0),chrom=character(0),chromStart=numeric(0),chromEnd=numeric(0),probes=numeric(0),means=numeric(0),calls=character(0),stringsAsFactors=FALSE)
|
|
|
27
|
|
|
28 for (chr in uniqueChr)
|
|
|
29 {
|
|
|
30 currentSubset=subset(CN, chromosome==chr)
|
|
|
31 currentPositions=currentSubset["position"]
|
|
|
32 for (sample in samples)
|
|
|
33 {
|
|
|
34 currentSignal=currentSubset[sample]
|
|
|
35 if (length(which(!is.na(unlist(currentSignal))))>1)
|
|
|
36 {
|
|
|
37 currentSeg=segmentation(signal=unlist(currentSignal),position=unlist(currentPositions),method=method)
|
|
|
38 callobj= callingObject(copynumber=currentSeg$signal, segmented=currentSeg$segmented,chromosome=rep(chr,length(currentSeg$signal)), position=currentSeg$startPos,sampleNames=sample)
|
|
|
39 currentCall=callingProcess(callobj,nclass=nbcall,cellularity=cellularity,verbose=TRUE)
|
|
|
40 currentResult=currentCall$segment
|
|
|
41 currentResult["sampleNames"]=c(rep(sample,length(currentCall$segment$chrom)))
|
|
|
42 result=rbind(result,currentResult)
|
|
|
43 }
|
|
|
44 }
|
|
|
45 }
|
|
|
46 finalResult=data.frame(sampleNames=result["sampleNames"],chrom=result["chrom"],chromStart=result["chromStart"],chromEnd=result["chromEnd"],probes=result["probes"],means=result["means"],calls=result["calls"],stringsAsFactors=FALSE)
|
|
|
47 sink(output)
|
|
|
48 print(format(finalResult))
|
|
|
49 sink()
|
|
|
50 #write.table(finalResult,output,row.names = FALSE, quote=FALSE, sep = "\t")
|
|
|
51 } else {
|
|
|
52 result=data.frame(sampleNames=character(0),chrom=character(0),start=numeric(0),end=numeric(0),points=numeric(0),means=numeric(0),stringsAsFactors=FALSE)
|
|
|
53
|
|
|
54 for (chr in uniqueChr)
|
|
|
55 {
|
|
|
56 cat(paste0("chromosome ",chr,"\n"))
|
|
|
57 currentSubset=subset(CN, chromosome==chr)
|
|
|
58 currentPositions=currentSubset["position"]
|
|
|
59 for (sample in samples)
|
|
|
60 {
|
|
|
61 cat(paste0(" sample ",sample,"..."))
|
|
|
62 currentSignal=currentSubset[sample]
|
|
|
63 if (length(which(!is.na(unlist(currentSignal))))>1)
|
|
|
64 {
|
|
|
65 currentSeg=segmentation(signal=unlist(currentSignal),position=unlist(currentPositions),method=method)
|
|
|
66 currentResult=currentSeg$segment
|
|
|
67 currentResult["chrom"]=c(rep(chr,length(currentSeg$segment$means)))
|
|
|
68 currentResult["sampleNames"]=c(rep(sample,length(currentSeg$segment$means)))
|
|
|
69 result=rbind(result,currentResult)
|
|
|
70
|
|
|
71 }
|
|
|
72 cat(paste0("OK\n"))
|
|
|
73 }
|
|
|
74 }
|
|
|
75 finalResult=data.frame(sampleNames=result["sampleNames"],chrom=result["chrom"],chromStart=result["start"],chromEnd=result["end"],probes=result["points"],means=result["means"],stringsAsFactors=FALSE)
|
|
|
76 sink(output)
|
|
|
77 print(format(finalResult))
|
|
|
78 sink()
|
|
|
79 #write.table(finalResult,output,row.names = FALSE, quote=FALSE, sep = "\t")
|
|
|
80 }
|
|
|
81
|
