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