vegaMC_R <- function(mtx, output_file_name="output",
beta=0.5, min_region_bp_size=1000, correction=FALSE,
loss_threshold=-0.2, gain_threshold=0.2, baf=FALSE,
loh_threshold=0.75, loh_frequency=0.8, bs=1000,
pval_threshold=0.05, html=TRUE, getGenes=TRUE,
mart_database="ensembl",
ensembl_dataset="hsapiens_gene_ensembl"){
if( output_file_name == "" ||
substr(output_file_name,
nchar(output_file_name), nchar(output_file_name)) == "/" ){
message("ERROR: Invalid Output File Name")
return(FALSE)
}
n_samples=0
n_probes=0
n_chromosomes=0
if(baf){
baf = 1
}else{
baf=0
}
mtx <- as.matrix(mtx[,-1])
res <- .C("run_vegaMC",
M=as.double(mtx),
row=as.integer(nrow(mtx)),
col=as.integer(ncol(mtx)),
out=as.character(output_file_name),
b= as.double(beta),
mrbs = as.integer(min_region_bp_size),
losst = as.double(loss_threshold),
gaint = as.double(gain_threshold),
ba = as.integer(baf),
loht = as.double(loh_threshold),
lohf = as.double(loh_frequency),
bsp = as.integer(bs),
ns = as.integer(n_samples),
np = as.integer(n_probes),
nc = as.integer(n_chromosomes))
n_samples = res$ns
n_probes = res$np
n_chromosomes = res$nc
segmentation <- read.table(output_file_name, sep="\t",
header=TRUE, as.is=TRUE)
segmentation[which(is.na(segmentation))] <- 0
if(correction==TRUE){
segmentation[,6] <- qvalue(as.numeric(segmentation[,6]))
segmentation[,7] <- qvalue(as.numeric(segmentation[,7]))
segmentation[,8] <- qvalue(as.numeric(segmentation[,8]))
}
segmentation[,6] <- round(as.numeric(segmentation[,6]), 5)
segmentation[,7] <- round(as.numeric(segmentation[,7]), 5)
segmentation[,8] <- round(as.numeric(segmentation[,8]), 5)
f_l <-
as.numeric(segmentation[,9]) * abs(
as.numeric(segmentation[,13]))/ as.numeric(segmentation[,12 ])
f_g <-
as.numeric(segmentation[,10]) * abs(
as.numeric(segmentation[,14]))/ as.numeric(segmentation[,12 ])
f_loh<-
as.numeric(segmentation[,11]) * abs(
as.numeric(segmentation[,15]))/ as.numeric(segmentation[,12 ])
segmentation <- cbind(segmentation, f_l, f_g, f_loh)
segmentation[which(is.na(segmentation))] <- 0
segmentation[,9] <- round(as.numeric(segmentation[,9])*100, 1)
segmentation[,10] <- round(as.numeric(segmentation[,10])*100, 1)
segmentation[,11] <- round(as.numeric(segmentation[,11])*100, 1)
segmentation[,9] <- paste(as.numeric(segmentation[,9]), "%", sep="")
segmentation[,10] <- paste(as.numeric(segmentation[,10]), "%", sep="")
segmentation[,11] <- paste(as.numeric(segmentation[,11]), "%", sep="")
colnames(segmentation)[12:18] <- c("Probe Size", "Loss Mean", "Gain Mean",
"LOH Mean", "Focal-score Loss",
"Focal-score Gain", "Focal-score LOH")
#segmentation[,1] <- gsub("X", "23", segmentation[,1])
#maxINT <- 2147483647
#ind_overflow <- which(as.numeric(segmentation[,2])<1)
#segmentation[ind_overflow,2] <- segmentation[ind_overflow,2] + (2 * (maxINT + 1))
#segmentation[ind_overflow,3] <- segmentation[ind_overflow,3] + (2 * (maxINT + 1))
write.table(segmentation, output_file_name, sep="\t", row.names=FALSE,
col.names=TRUE, quote=FALSE, eol="\n")
return(segmentation)
}
#' getBreaksVegaMC Get SCEVAN segmentation of the matrix.
#'
#' @param mtx count matrix
#' @param chr_vect Vector specifying for each gene the chromosome where it is located
#' @param sample sample name (optional)
#' @param beta_vega specifies beta parameter for segmentation, higher beta for more coarse-grained segmentation. (default 0.5)
#'
#' @return breakpoints
#' @export
#'
#' @examples
getBreaksVegaMC <- function(mtx, chr_vect, sample = "", beta_vega = 0.5){
#write.table(mtx, paste("./output/", sample, "_mtx_vega.txt", sep=""), sep="\t", row.names = FALSE, quote = F)
res_vega <- vegaMC_R(mtx = mtx, output_file_name=paste("./output/", sample,"vega_output"), beta = beta_vega);
BR <- unlist(lapply(res_vega$Start, function(x) which(chr_vect == x)[1]))
n <- nrow(mtx)
BR <- sort(unique(c(1, BR, n)))
return(BR)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.