|
8
|
1 # Load ExomeDepth library (without warnings)
|
|
|
2 suppressMessages(library(ExomeDepth))
|
|
|
3
|
|
|
4 # Import parameters
|
|
|
5 args <- commandArgs(trailingOnly=TRUE)
|
|
|
6 eval(parse(text=args[[1]]))
|
|
|
7 param <- read.table(mypars,sep="=", as.is=TRUE)
|
|
|
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], MAR = 1, FUN = sum)
|
|
|
55
|
|
|
56 # Now create the ExomeDepth object for the CNVs call
|
|
|
57 all.exons<-suppressWarnings(suppressMessages(new('ExomeDepth',
|
|
|
58 test = ExomeCount.mat[,i],
|
|
|
59 reference = my.reference.selected,
|
|
|
60 formula = 'cbind(test,reference)~1')))
|
|
|
61
|
|
|
62
|
|
|
63 # Now call the CNVs
|
|
|
64 result <- try(all.exons<-suppressMessages(CallCNVs(x=all.exons,
|
|
|
65 transition.probability = trans_prob,
|
|
|
66 chromosome = ExomeCount.dafr$space,
|
|
|
67 start = ExomeCount.dafr$start,
|
|
|
68 end = ExomeCount.dafr$end,
|
|
|
69 name = ExomeCount.dafr$names)), silent=T)
|
|
|
70
|
|
|
71 # Next if CNVs are not detected
|
|
|
72 if (class(result)=="try-error"){
|
|
|
73 next
|
|
|
74 }
|
|
|
75
|
|
|
76 # Compute correlation between ref and test
|
|
|
77 my.cor <- cor(all.exons@reference, all.exons@test)
|
|
|
78 n.call <- nrow(all.exons@CNV.calls)
|
|
|
79
|
|
|
80 # Write results
|
|
|
81 my.results <- cbind(all.exons@CNV.calls[,c(7,5,6,3)],
|
|
|
82 sample=colnames(ExomeCount.mat)[i],
|
|
|
83 corr=my.cor,
|
|
|
84 all.exons@CNV.calls[,c(4,9,12)])
|
|
|
85
|
|
|
86 # Re-order by chr and position
|
|
|
87 chrOrder<-c(paste("chr",1:22,sep=""),"chrX","chrY","chrM")
|
|
|
88 my.results[,1] <- factor(my.results[,1], levels=chrOrder)
|
|
|
89 my.results <- my.results[order(my.results[,1], my.results[,2], my.results[,3]),]
|
|
|
90
|
|
|
91 write.table(sep='\t', quote=FALSE, file = output,
|
|
|
92 x = my.results,
|
|
|
93 row.names = FALSE, col.names = FALSE, dec=".", append=TRUE)
|
|
|
94 }
|
