Nothing
### This function detects chromosomal breakpoints along genome
### Copyright (C) 2003 Institut Curie
### Author(s): Philippe Hupé (Institut Curie) 2003
### Contact: glad@curie.fr
chrBreakpoints <- function(...)
{
UseMethod("chrBreakpoints")
}
chrBreakpoints.profileCGH <- function(profileCGH, smoothfunc="lawsglad", OnlyOptimCall = FALSE, base=FALSE, sigma=NULL,
model="Gaussian", bandwidth=10, round=1.5, verbose=FALSE,
breaksFdrQ = 0.0001, haarStartLevel = 1, haarEndLevel = 5, weights.name = NULL, ...)
{
if (verbose)
{
print("chrBreakpoints: starting function")
call <- match.call()
print(paste("Call function:", call))
}
if (smoothfunc != "laws" &&
smoothfunc != "aws" &&
smoothfunc != "lawsglad" &&
smoothfunc != "haarseg" &&
smoothfunc != "haarsegbychr")
{
stop("Choose either aws, laws or haarseg for smoothfunc")
}
if (base == TRUE)
{
if(smoothfunc != "lawsglad" && smoothfunc != "haarseg")
{
stop("Choose either aws, or laws when base=TRUE")
}
}
if (is.null(sigma))
{
resetsigma <- TRUE
}
else
{
resetsigma <- FALSE
}
### le round est fait a de nombreux endroits différents: cela peut
### surement être changé
#### local functions
IQRdiff <- function(y) IQR(diff(y)) / 1.908
roundglad <- function(x, digits=2)
{
## en C return floor(n * pow(10., d) + .5) / pow(10., d);
dec <- (digits-trunc(digits))
if (dec==0)
dec <- 1
r <- 10^(-trunc(digits))*dec
x <- r*round(x/r)
return(x)
}
## roundtest <- function(x, d=2)
## {
## return(floor(x * 10^d + .5) / 10^d)
## }
if(smoothfunc != "haarseg")
{
indice <- 1:length(profileCGH$profileValues[[1]])
ChrIndice <- split(indice,profileCGH$profileValues[["Chromosome"]])
ChrName <- names(ChrIndice)
if (is.numeric(profileCGH$profileValues[["Chromosome"]][1]))
{
labelChr <- as.numeric(ChrName)
}
else
{
labelChr <- as.character(ChrName)
}
NbChr <- length(labelChr)
lg <- rep(0,NbChr)
PosOrderRange <- data.frame(Chromosome = lg, MinPosOrder = lg, MaxPosOrder = lg)
## initialization of region number to 0
nbregion <- 0
## initialization of level number to 0
nblevel <- 0
profileCGH$BkpDetected <- data.frame(Chromosome = as.integer(labelChr), BkpDetected = 0)
IQRvalue <- IQRChr <- NULL
for (i in 1:NbChr)
{
if (verbose) print(paste("chrBreakpoints: starting chromosome", labelChr[i]))
## location of data related to each chromosome
indexChr <- ChrIndice[[i]]
subsetdata <- profileCGH$profileValues[indexChr,]
## information sur les bornes d'un chromosome
PosOrderRange$Chromosome[i] <- labelChr[i]
rangePos <- range(subsetdata[["PosOrder"]])
subsetdata[["MinPosOrder"]] <- PosOrderRange$MinPosOrder[i] <- rangePos[1]
subsetdata[["MaxPosOrder"]] <- PosOrderRange$MaxPosOrder[i] <- rangePos[2]
if (length(indexChr)>1)
{
if (resetsigma)
{
sigma <- IQRdiff(subsetdata[["LogRatio"]])^2
IQRvalue[i] <- sigma^(0.5)
IQRChr[i] <- labelChr[i]
}
else
{
IQRvalue[i] <- sigma^0.5
IQRChr[i] <- labelChr[i]
}
if(sigma == 0)
{
print("Warnings: sigma equal 0")
print("Number of probes:", length(indexChr))
print("sigma was automatically set to 1")
sigma <- 1
}
if (base == TRUE)
{
x <- subsetdata[["PosBase"]]
datarange <- range(x)
hmax <- diff(datarange)*bandwidth
hinit <- median(diff(x))
## smoothing
if (smoothfunc == "laws")
{
dim(x) <- c(1,length(x)) #à supprimer dans la nouvelle version du package AWS
awsres <- aws::laws(y=subsetdata[["LogRatio"]], x=x,
hinit=hinit, hmax=hmax, shape=sigma, NN=FALSE,
symmetric=TRUE, model=model, ...)$theta
if (is.null(awsres) == FALSE)
{
subsetdata[["Smoothing"]] <- roundglad(awsres, round)
}
else
{
subsetdata[["Smoothing"]] <- 99999
}
}
if (smoothfunc == "aws")
{
awsres <- aws::aws(y=subsetdata[["LogRatio"]], x=x, hinit=hinit, hmax=hmax, sigma2=sigma, NN=FALSE, symmetric=TRUE, ...)$theta
if (is.null(awsres) == FALSE)
{
subsetdata[["Smoothing"]] <- roundglad(awsres, round)
}
else
{
subsetdata[["Smoothing"]] <- 99999
}
}
}
else
{
hinit <- 1
hmax <- length(subsetdata[["PosOrder"]])*bandwidth #si on laisse hmax à la valeur de length(data$PosOrder[indexChr]) ce n'est pas assez et les créneaux ne sont pas bien fittés
if (smoothfunc == "lawsglad")
{
awsres <- lawsglad(y=subsetdata[["LogRatio"]],
hinit=hinit, hmax=hmax, shape=sigma, model=model, ...)
if(is.null(awsres) == FALSE)
{
subsetdata[["Smoothing"]] <- roundglad(awsres, round)
}
else
{
subsetdata[["Smoothing"]] <- 99999
}
}
if (smoothfunc == "laws")
{
awsres <- aws::laws(y=subsetdata[["LogRatio"]], hinit=hinit,
hmax=hmax, shape=sigma, symmetric=TRUE, model=model,
...)$theta
if(is.null(awsres) == FALSE)
{
subsetdata[["Smoothing"]] <- roundglad(awsres, round)
}
else
{
subsetdata[["Smoothing"]] <- 99999
}
}
if (smoothfunc == "aws")
{
awsres <- aws::aws(y=subsetdata[["LogRatio"]], hinit=hinit, hmax=hmax, sigma2=sigma, symmetric=TRUE, ...)$theta
if(is.null(awsres)==FALSE)
{
subsetdata[["Smoothing"]] <- roundglad(awsres, round)
}
else
{
subsetdata[["Smoothing"]] <- 99999
}
}
if (smoothfunc == "haarsegbychr")
{
awsres <- HaarSegGLADCPP(subsetdata[["LogRatio"]], breaksFdrQ=breaksFdrQ, haarStartLevel=haarStartLevel , haarEndLevel=haarEndLevel)
if (is.null(awsres) == FALSE)
{
subsetdata[["Smoothing"]] <- awsres
## subsetdata$Smoothing <- roundglad(awsres, round)
}
else
{
subsetdata[["Smoothing"]] <- 99999
}
}
}
nbregion <- nbregion + 1
l <- length(subsetdata[["Smoothing"]])
putLevel <- .C("putLevel_awsBkp",
## variables pour putLevel
Smoothing = as.double(subsetdata[["Smoothing"]]), ## valeur de sortie
as.double(subsetdata[["LogRatio"]]),
Level = integer(l), ## valeur de sortie
nblevel = as.integer(nblevel), ## valeur de sortie
as.integer(l),
## variables pour awsBkp
OutliersAws = integer(l), ## valeur de sortie
nbregion = as.integer(nbregion), ## valeur de sortie
regionChr = integer(l), ## valeur de sortie
Breakpoints = integer(l), ## valeur de sortie
BkpDetected = integer(1), ## valeur de sortie
PACKAGE="GLAD")
subsetdata[c("Smoothing", "Level", "Region", "OutliersAws", "Breakpoints")] <- putLevel[c("Smoothing", "Level", "regionChr", "OutliersAws", "Breakpoints")]
nblevel <- putLevel$nblevel
profileCGH$BkpDetected$BkpDetected[i] <- putLevel$BkpDetected
nbregion <- putLevel$nbregion
}
else
{
subsetdata[["Region"]] <- -1
subsetdata[["Level"]] <- -1
subsetdata[["Breakpoints"]] <- 0
IQRvalue[i] <- 0
IQRChr[i] <- labelChr[i]
}
profileCGH$profileValues[indexChr,] <- subsetdata
if (verbose) print(paste("chrBreakpoints: ending chromosome", labelChr[i]))
} ## fin de la boucle par chromosome
profileCGH$Sigma <- data.frame(Chromosome=IQRChr,Value=IQRvalue)
profileCGH$PosOrderRange <- PosOrderRange
}
else
{
if(!is.null(weights.name))
{
W <- profileCGH$profileValues[[weights.name]]
}
else
{
W <- NULL
}
NbChr <- length(unique(profileCGH$profileValues[["Chromosome"]]))
l <- length(profileCGH$profileValues[["LogRatio"]])
res <- .C("chrBreakpoints_haarseg",
as.double(profileCGH$profileValues[["LogRatio"]]),
as.integer(profileCGH$profileValues[["Chromosome"]]),
Smoothing = as.double(profileCGH$profileValues[["Smoothing"]]), ## valeur de sortie
Level = integer(l), ## valeur de sortie
OutliersAws = integer(l), ## valeur de sortie
regionChr = integer(l), ## valeur de sortie
Breakpoints = integer(l), ## valeur de sortie
sizeChr = integer(NbChr), ## taille de chaque chromosome
startChr = integer(NbChr), ## position pour le début des valeurs de chaque chromosome
IQRChr = integer(NbChr), ## numéro du chromosome pour le calcul de l'IQR
IQRValue = double(NbChr), ## valeur de l'IQR
BkpDetected = integer(NbChr), ## valeur de sortie
## paramètres pour Haarseg
as.double(breaksFdrQ),
as.integer(haarStartLevel),
as.integer(haarEndLevel),
as.integer(NbChr), ## nombre de chromosomes
as.integer(l),
as.double(W),
as.integer(OnlyOptimCall),
PACKAGE = "GLAD")
## ##############################
## récupération des informations
## ##############################
## Range de position pour chaque chromosome
MinPosOrder <- res[["startChr"]] + 1
MaxPosOrder <- res[["sizeChr"]] + MinPosOrder - 1
profileCGH$PosOrderRange <- data.frame(Chromosome = res[["IQRChr"]],
MinPosOrder = MinPosOrder,
MaxPosOrder = MaxPosOrder)
## valeur de l'IQR
profileCGH$Sigma <- data.frame(Chromosome = res[["IQRChr"]],
Value = res[["IQRValue"]])
## valeur BkpDetected
profileCGH$BkpDetected <- data.frame(Chromosome = res[["IQRChr"]],
BkpDetected = res[["BkpDetected"]])
## valeurs de segmentation
profileCGH$profileValues[c("Smoothing", "Level", "Region", "OutliersAws", "Breakpoints")] <- res[c("Smoothing", "Level", "regionChr", "OutliersAws", "Breakpoints")]
}
if (verbose) print("chrBreakpoints: ending function")
return(profileCGH)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.