Mercurial > repos > crs4 > exomedepth
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 |
| rev | line source |
|---|---|
| 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) | |
|
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 | 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 } |