annotate exomedepth.R @ 7:45af4a9748cf draft

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