Mercurial > repos > crs4 > exomedepth

comparison exomedepth-a9701c421408/exomedepth.R @ 0:6ea78cede7aa draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author crs4
date Tue, 05 Jun 2018 18:53:13 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:6ea78cede7aa
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 }