|
2
|
1 # Load ExomeDepth library (without warnings)
|
|
|
2 suppressMessages(library(ExomeDepth))
|
|
|
3
|
|
|
4 # Import parameters from xml wrapper (args_file)
|
|
|
5 args <- commandArgs(trailingOnly=TRUE)
|
|
4
|
6 eval(parse(text=args[[1]]))
|
|
|
7 param <- read.table(mypars,sep="=", as.is=TRUE)
|
|
2
|
8
|
|
|
9 # Set common parameters
|
|
|
10 target <- param[match("target",param[,1]),2]
|
|
|
11 trans_prob <- as.numeric(param[match("trans_prob",param[,1]),2])
|
|
|
12 output <- param[match("output",param[,1]),2]
|
|
|
13 test_vs_ref <- as.logical(param[match("test_vs_ref",param[,1]),2])
|
|
|
14
|
|
|
15 # Create symbolic links for multiple bam and bai
|
|
|
16 bam <- param[param[,1]=="bam",2]
|
|
|
17 bam_bai <- param[param[,1]=="bam_bai",2]
|
|
|
18 bam_label <- param[param[,1]=="bam_label",2]
|
|
|
19 bam_label <- gsub(" ", "_", bam_label)
|
|
|
20
|
|
|
21 for(i in 1:length(bam)){
|
|
|
22 stopifnot(file.symlink(bam[i], paste(bam_label[i], "bam", sep=".")))
|
|
|
23 stopifnot(file.symlink(bam_bai[i], paste(bam_label[i], "bam.bai", sep=".")))
|
|
|
24 }
|
|
|
25
|
|
|
26 # Generate read count data
|
|
|
27 BAMFiles <- paste(bam_label, "bam", sep=".")
|
|
|
28 sink("/dev/null")
|
|
|
29 ExomeCount <- suppressMessages(getBamCounts(bed.file=target, bam.files = BAMFiles))
|
|
|
30 sink()
|
|
|
31
|
|
|
32 # Convert counts in a data frame
|
|
|
33 ExomeCount.dafr <- as(ExomeCount[, colnames(ExomeCount)], 'data.frame')
|
|
|
34
|
|
|
35 # Prepare the main matrix of read count data
|
|
|
36 ExomeCount.mat <- as.matrix(ExomeCount.dafr[, grep(names(ExomeCount.dafr), pattern='.bam')])
|
|
|
37
|
|
|
38 # Remove .bam from sample name
|
|
|
39 colnames(ExomeCount.mat) <- gsub(".bam", "", colnames(ExomeCount.mat))
|
|
|
40
|
|
|
41 # Set nsamples == 1 if mode is test vs reference, assuming test is sample 1
|
|
|
42 nsamples <- ifelse(test_vs_ref, 1, ncol(ExomeCount.mat))
|
|
|
43
|
|
|
44 # Loop over samples
|
|
|
45 for (i in 1:nsamples){
|
|
|
46
|
|
|
47 # Create the aggregate reference set for this sample
|
|
|
48 my.choice <- suppressWarnings(suppressMessages(
|
|
|
49 select.reference.set(test.counts = ExomeCount.mat[,i],
|
|
|
50 reference.counts = subset(ExomeCount.mat, select=-i),
|
|
|
51 bin.length = (ExomeCount.dafr$end - ExomeCount.dafr$start)/1000,
|
|
|
52 n.bins.reduced = 10000)))
|
|
|
53
|
|
|
54 my.reference.selected <- apply(X = ExomeCount.mat[, my.choice$reference.choice, drop=FALSE],
|
|
|
55 MAR = 1,
|
|
|
56 FUN = sum)
|
|
|
57
|
|
|
58 # Now create the ExomeDepth object for the CNVs call
|
|
|
59 all.exons<-suppressWarnings(suppressMessages(new('ExomeDepth',
|
|
|
60 test = ExomeCount.mat[,i],
|
|
|
61 reference = my.reference.selected,
|
|
|
62 formula = 'cbind(test,reference)~1')))
|
|
|
63
|
|
|
64
|
|
|
65 # Now call the CNVs
|
|
|
66 result <- try(all.exons<-suppressMessages(CallCNVs(x=all.exons,
|
|
|
67 transition.probability = trans_prob,
|
|
|
68 chromosome = ExomeCount.dafr$space,
|
|
|
69 start = ExomeCount.dafr$start,
|
|
|
70 end = ExomeCount.dafr$end,
|
|
|
71 name = ExomeCount.dafr$names)), silent=T)
|
|
|
72
|
|
|
73 # Next if CNVs are not detected
|
|
|
74 if (class(result)=="try-error"){
|
|
|
75 next
|
|
|
76 }
|
|
|
77
|
|
|
78 # Compute correlation between ref and test
|
|
|
79 my.cor <- cor(all.exons@reference, all.exons@test)
|
|
|
80 n.call <- nrow(all.exons@CNV.calls)
|
|
|
81
|
|
|
82 # Write results
|
|
|
83 my.results <- cbind(all.exons@CNV.calls[,c(7,5,6,3)],
|
|
|
84 sample=colnames(ExomeCount.mat)[i],
|
|
|
85 corr=my.cor,
|
|
|
86 all.exons@CNV.calls[,c(4,9,12)])
|
|
|
87
|
|
|
88 # Re-order by chr and position
|
|
|
89 chrOrder<-c(paste("chr",1:22,sep=""),"chrX","chrY","chrM")
|
|
|
90 my.results[,1] <- factor(my.results[,1], levels=chrOrder)
|
|
|
91 my.results <- my.results[order(my.results[,1], my.results[,2], my.results[,3]),]
|
|
|
92
|
|
|
93 write.table(sep='\t', quote=FALSE, file = output,
|
|
|
94 x = my.results,
|
|
|
95 row.names = FALSE, col.names = FALSE, dec=".", append=TRUE)
|
|
|
96 }
|
