R/allMD.R
In SoniaTC/NOISeq: Exploratory analysis and differential expression for RNA-seq data
allMD <- function (input, factor, conditions, k = 0.5, replicates, norm = "rpkm",
pnr = 0.2, nss = 5, v = 0.02, lc = 0)
# input: Set of data of type Input
# conditions: Levels of the factor to be compared (when the factor has more than 2 levels)
# k: When counts = 0, 0 will be changed to k. By default, k = 0.5.
# norm: Normalization method. It can be one of "rpkm" (default), "uqua"
# (upper quartile), "tmm" (trimmed mean of M) or "n" (no normalization).
# pnr: Percentage of total reads (seq.depth) for each simulated sample.
# Only needed when noise = "simul". By default, pnr = 1.
# nss: Number of simulated samples (>= 2). By default, nss = 5.
# If nss = 0, real samples are used to compute noise.
# v: Variability in sample total reads used to simulate samples.
# By default, v = 0.02. Sample total reads is computed as a
# random value from a uniform distribution in the interval
# [(pnr-v)*sum(counts), (pnr+v)*sum(counts)]
# lc: Length correction in done by dividing expression by length^lc.
# By default, lc = 1.
{
# n1 <- ncol(as.matrix(datos1))
# n2 <- ncol(as.matrix(datos2))
# g.sinL <- names(which(is.na(long)))
# Check if the factor introduced is already defined
# If the factor introduced is defined and has more than 2 conditions, it will check if the conditions specified are defined too
condition_fac = FALSE
condition_lev = FALSE
datos1 <- datos2 <- matrix()
for (i in colnames(pData(input))) {
if (factor == i) {
condition_fac = TRUE
if (!is.factor(pData(input)[,i])) pData(input)[,i] = as.factor(pData(input)[,i])
if (length(levels(pData(input)[,i])) == 2) {
if (!is.null(assayData(input)$exprs)) {
datos1 <- assayData(input)$exprs[,which(pData(input)[,i] ==levels(pData(input)[,i])[1]), drop = FALSE]
datos2 <- assayData(input)$exprs[,which(pData(input)[,i] ==levels(pData(input)[,i])[2]), drop = FALSE]
} else {
datos1 <- assayData(input)$counts[,which(pData(input)[,i] ==levels(pData(input)[,i])[1]), drop = FALSE]
datos2 <- assayData(input)$counts[,which(pData(input)[,i] ==levels(pData(input)[,i])[2]), drop = FALSE]
}
# Define the comparison string
comparison <- paste(levels(pData(input)[,i])[1], levels(pData(input)[,i])[2], sep=" - ")
condition_lev = TRUE
}
else {
if (is.null(conditions))
stop("Error. You must specify which conditions you wish to compare when the factor has two or more conditions.\n")
if (length(conditions) != 2)
stop("Error. The argument conditions must contain the 2 conditions you wish to compare.")
l <- conditions %in% pData(input)[,i]
# If they are defined, they will be TRUE
if (l[1] == TRUE && l[2] == TRUE) {
if (!is.null(assayData(input)$exprs)) {
datos1 <- assayData(input)$exprs[,which(pData(input)[,i] == conditions[1]), drop = FALSE]
datos2 <- assayData(input)$exprs[,which(pData(input)[,i] == conditions[2]), drop = FALSE]
} else {
datos1 <- assayData(input)$counts[,which(pData(input)[,i] == conditions[1]), drop = FALSE]
datos2 <- assayData(input)$counts[,which(pData(input)[,i] == conditions[2]), drop = FALSE]
}
# Define the comparison string
comparison <- paste(conditions[1],conditions[2], sep=" - ")
condition_lev = TRUE
}
}
}
}
if (condition_fac == FALSE)
stop("The factor you have written does not correspond with any of the ones you have defined.")
if (condition_lev == FALSE)
stop("The conditions you have written don't exist in the factor specified.\n")
# Correction to make it work when there are simulated samples
if (replicates == "no")
replicates = "technical"
# if (description(input)@samples[[1]] == "no")
# description(input)@samples[[1]] = "technical"
n1 <- ncol(as.matrix(datos1))
n2 <- ncol(as.matrix(datos2))
if (norm == "n") { # no normalization
datos1 <- round(datos1, 100)
datos2 <- round(datos2, 100)
}
if (is.null(k)) {
m1 <- min(datos1[noceros(datos1, num = FALSE)], na.rm = TRUE)
m2 <- min(datos2[noceros(datos2, num = FALSE)], na.rm = TRUE)
mm <- min(m1, m2)
k <- mm/2
}
# Total counts for each gene:
suma1 <- rowSums(as.matrix(datos1))
suma2 <- rowSums(as.matrix(datos2))
# All genes
todos <- rownames(as.matrix(datos1))
# Genes with counts in any condition
concounts <- names(which(suma1+suma2 > 0))
long <- 1000
g.sinL <- NULL
if (!is.null(featureData(input)@data$Length)) {
g.sinL <- names(which(is.na(featureData(input)@data$Length)))
if (any(!is.na(featureData(input)@data$Length)) == TRUE)
long <- featureData(input)@data[concounts, "Length"]
}
if (replicates == "technical") { ### technical replicates
suma1 <- suma1[concounts]
suma2 <- suma2[concounts]
#-------------------------------------------------------------------------#
# Normalization of counts for each condition (aggregating replicates)
if (norm == "rpkm") { # RPKM
suma1.norm <- rpkm(suma1, long = long, k = k, lc = lc)
suma2.norm <- rpkm(suma2, long = long, k = k, lc = lc)
}
if (norm == "uqua") {
suma.norm <- uqua(cbind(suma1, suma2), long = long, lc = lc, k = k)
suma1.norm <- as.matrix(suma.norm[ ,1])
suma2.norm <- as.matrix(suma.norm[ ,2])
}
if (norm == "tmm") {
suma.norm <- tmm(as.matrix(cbind(suma1, suma2)), long = long,
lc = lc, k = k)
suma1.norm <- as.matrix(suma.norm[ ,1])
suma2.norm <- as.matrix(suma.norm[ ,2])
}
}
#-------------------------------------------------------------------------#
## Noise distribution
if ((n1+n2)>2) { # with real samples
datitos <- cbind(datos1, datos2)
datitos <- datitos[concounts,]
gens.sin0 <- setdiff(concounts, g.sinL)
if (norm == "n") { # no normalization
datitos.0 <- sinceros(datitos, k = k)
datitos.norm <- datitos.0[gens.sin0, ]
}
if (norm == "rpkm") { # RPKM
datitos.0 <- rpkm(datitos, long = long, k = k, lc = lc)
datitos.norm <- datitos.0[gens.sin0, ]
}
if (norm == "uqua") { # Upper Quartile
datitos.0 <- uqua(datitos, long = long, lc = lc, k = k)
datitos.norm <- datitos.0[gens.sin0, ]
}
if (norm == "tmm") {
datitos.0 <- tmm(datitos, long = long, lc = lc, k = k)
datitos.norm <- datitos.0[gens.sin0, ]
}
datos1.norm <- datitos.norm[ ,1:n1]
datos2.norm <- datitos.norm[ ,(n1+1):(n1+n2)]
if (n1 > 1) {
MD1 <- MD(dat = datos1.norm)
} else { MD1 <- NULL }
if (n2 > 1) {
MD2 <- MD(dat = datos2.norm)
} else { MD2 <- NULL }
} else { # with simulated samples
if (nss == 0) {
nss <- 5
}
datos.sim <- sim.samples(counts1 = sinceros(suma1, k = k),
counts2 = sinceros(suma2, k = k),
pnr = pnr, nss = nss, v = v)
nn <- sapply(datos.sim, ncol)
dat.sim.norm <- vector("list", length = 2)
datitos <- cbind(datos.sim[[1]], datos.sim[[2]])
rownames(datitos) = names(suma1)
sumita <- rowSums(datitos)
g.sin0 <- names(which(sumita > 0))
gens.sin0 <- setdiff(g.sin0, g.sinL)
if (norm == "n") { # no normalization
datitos.0 <- sinceros(datitos, k = k)
datitos.norm <- datitos.0[gens.sin0, ]
}
if (norm == "rpkm") { # RPKM
datitos.0 <- rpkm(datitos, long = long, k = k, lc = lc)
datitos.norm <- datitos.0[gens.sin0, ]
}
if (norm == "uqua") { # Upper Quartile
datitos.0 <- uqua(datitos, long = long, lc = lc, k = k)
datitos.norm <- datitos.0[gens.sin0, ]
}
if (norm == "tmm") {
datitos.0 <- tmm(datitos, long = long, lc = lc, k = k)
datitos.norm <- datitos.0[gens.sin0, ]
}
dat.sim.norm[[1]] <- datitos.norm[ ,1:nn[1]]
dat.sim.norm[[2]] <- datitos.norm[ ,(nn[1]+1):sum(nn)]
MD1 <- MD(dat = dat.sim.norm[[1]])
MD2 <- MD(dat = dat.sim.norm[[2]])
}
Mr <- c(as.numeric(MD1$M), as.numeric(MD2$M))
Dr <- c(as.numeric(MD1$D), as.numeric(MD2$D))
#-------------------------------------------------------------------------#
## M and D for different experimental conditions
if (replicates == "technical" & norm != "n") {
MDs <- MD(dat = cbind(suma1.norm, suma2.norm))
lev1 <- suma1.norm[,1]
lev1 <- lev1[todos]
lev2 <- suma2.norm[,1]
lev2 <- lev2[todos]
} else {
if ((n1+n1) == 2) {
datos1.norm <- sinceros(as.matrix(datos1)[concounts,], k = k)
datos2.norm <- sinceros(as.matrix(datos2)[concounts,], k = k)
}
resum1.norm <- rowMeans(as.matrix(datos1.norm))
resum2.norm <- rowMeans(as.matrix(datos2.norm))
lev1 <- resum1.norm[todos]
lev2 <- resum2.norm[todos]
MDs <- MD(dat = cbind(resum1.norm, resum2.norm))
}
## Completing M and D
names(lev1) <- names(lev2) <- todos
Ms <- as.numeric(MDs$M)
names(Ms) <- rownames(MDs$M)
Ms <- Ms[todos]
names(Ms) <- todos
Ds <- as.numeric(MDs$D)
names(Ds) <- rownames(MDs$D)
Ds <- Ds[todos]
names(Ds) <- todos
## Results
list("k" = k, "comp" = comparison, "Level1" = lev1, "Level2" = lev2, "Ms" = Ms, "Ds" = Ds, "Mn" = Mr, "Dn" = Dr)
}
#######################################################################
#######################################################################
allMDbio = function (input, factor, conditions, k = 0.5, norm = "rpkm", lc = 1,
r = 10, a0per = 0.9, nclust = 15, filter = 1, depth = NULL,
cv.cutoff = 0, cpm = 1)
# input: Set of data of type Input
# conditions: Levels of the factor to be compared (when the factor has more than 2 levels)
# k: When counts = 0, 0 will be changed to k. By default, k = 0.5.
# norm: Normalization method. It can be one of "rpkm" (default), "uqua"
# (upper quartile), "tmm" (trimmed mean of M) or "n" (no normalization).
# lc: Length correction in done by dividing expression by length^lc.
# By default, lc = 1.
# r: Number of permutations to compute null distribution (r=10).
# a0per: Percentile of S to compute a0. If NULL, a0 = 0. (a0per = 0.9)
{
# Check if the factor introduced is already defined
# If the factor introduced is defined and has more than 2 conditions,
# it will check if the conditions specified are defined too
condition_fac = FALSE
condition_lev = FALSE
datos1 <- datos2 <- matrix()
for (i in colnames(pData(input))) {
if (factor == i) {
condition_fac = TRUE
if (!is.factor(pData(input)[,i])) pData(input)[,i] = as.factor(pData(input)[,i])
if (length(levels(pData(input)[,i])) == 2) {
if (!is.null(assayData(input)$exprs)) {
datos1 <- assayData(input)$exprs[,which(pData(input)[,i] ==levels(pData(input)[,i])[1])]
datos2 <- assayData(input)$exprs[,which(pData(input)[,i] ==levels(pData(input)[,i])[2])]
} else {
datos1 <- assayData(input)$counts[,which(pData(input)[,i] ==levels(pData(input)[,i])[1])]
datos2 <- assayData(input)$counts[,which(pData(input)[,i] ==levels(pData(input)[,i])[2])]
}
# Define the comparison string
comparison <- paste(levels(pData(input)[,i])[1], levels(pData(input)[,i])[2], sep=" - ")
condition_lev = TRUE
if (!((ncol(datos1) > 1) && (ncol(datos2) > 1)))
stop("Error. NOISeqBIO needs at least 2 biological replicates per condition.\n")
}
else {
if (is.null(conditions))
stop("Error. You must specify which conditions you wish to compare when the factor has two or more conditions.\n")
if (length(conditions) != 2)
stop("Error. The argument conditions must contain the 2 conditions you wish to compare.")
l <- conditions %in% pData(input)[,i]
# If they are defined, they will be TRUE
if (l[1] == TRUE && l[2] == TRUE) {
if (!is.null(assayData(input)$exprs)) {
datos1 <- assayData(input)$exprs[,which(pData(input)[,i] == conditions[1])]
datos2 <- assayData(input)$exprs[,which(pData(input)[,i] == conditions[2])]
} else {
datos1 <- assayData(input)$counts[,which(pData(input)[,i] == conditions[1])]
datos2 <- assayData(input)$counts[,which(pData(input)[,i] == conditions[2])]
}
# Define the comparison string
comparison <- paste(conditions[1],conditions[2], sep=" - ")
condition_lev = TRUE
}
}
}
}
if (condition_fac == FALSE)
stop("The factor specified does not correspond with any of the ones you have defined.")
if (condition_lev == FALSE)
stop("The conditions specified don't exist for the factor specified.\n")
##-------------------------------------------------------------------------##
# Number of observations within each condition
n1 <- ncol(as.matrix(datos1))
n2 <- ncol(as.matrix(datos2))
if (max(n1,n2) == 1)
stop("There is only one replicate per condition. Please, use NOISeq instead of NOISeqBIO.\n")
# Rounding off data
if (norm == "n") { # no normalization
datos1 <- round(datos1, 10)
datos2 <- round(datos2, 10)
}
# Computing k
if (is.null(k)) {
m1 <- min(datos1[noceros(datos1, num = FALSE)], na.rm = TRUE)
m2 <- min(datos2[noceros(datos2, num = FALSE)], na.rm = TRUE)
k <- min(m1, m2)/2
}
# Total counts for each gene:
suma1 <- rowSums(as.matrix(datos1))
suma2 <- rowSums(as.matrix(datos2))
# Genes with counts in any condition
concounts <- names(which(suma1+suma2 > 0))
# All genes
todos <- rownames(as.matrix(datos1))
# Gene length
long <- 1000
g.sinL <- NULL # genes with no length defined
if (!is.null(featureData(input)@data$Length)) {
g.sinL <- names(which(is.na(featureData(input)@data$Length)))
if (any(!is.na(featureData(input)@data$Length)) == TRUE)
long <- featureData(input)@data[concounts, "Length"]
}
# Genes with counts and with length
gens.sin0 <- setdiff(concounts, g.sinL)
# cond1 and cond2 in the same matrix
datitos <- cbind(datos1, datos2)
datitos <- datitos[concounts,] # selecting only genes with counts
# Sequencing depth when filtering method = 3
if (filter == 3 && is.null(depth)) depth = colSums(datitos)
#-------------------------------------------------------------------------#
#-------------------------------------------------------------------------#
## Normalization
if (norm == "n") { # no normalization
datitos.0 <- sinceros(datitos, k = k)
datitos.norm <- datitos.0[gens.sin0, ]
}
if (norm == "rpkm") { # RPKM
datitos.0 <- rpkm(datitos, long = long, k = k, lc = lc)
datitos.norm <- datitos.0[gens.sin0, ]
}
if (norm == "uqua") { # Upper Quartile
datitos.0 <- uqua(datitos, long = long, lc = lc, k = k)
datitos.norm <- datitos.0[gens.sin0, ]
}
if (norm == "tmm") {
datitos.0 <- tmm(datitos, long = long, lc = lc, k = k)
datitos.norm <- datitos.0[gens.sin0, ]
}
#-------------------------------------------------------------------------#
## Filtering out low count features
if (filter != 0) {
datos.filt = filtered.data(dataset = datitos.norm, factor = c(rep("cond1", n1), rep("cond2", n2)),
norm = TRUE, depth = depth, method = filter, cv.cutoff = cv.cutoff, cpm = cpm)
} else {
datos.filt = datitos.norm
}
datos1.filt <- datos.filt[ ,1:n1]
datos2.filt <- datos.filt[ ,(n1+1):(n1+n2)]
#-------------------------------------------------------------------------#
## Noise distribution
Zr = NULL
if (n1+n2 <= 9) { # sharing information within clusters
Zr = share.info(mydata = datos.filt, n1 = n1, n2 = n2, r = r, nclust = nclust)
} else { # r permutations
for (i in 1:r) {
print(paste("r =", i))
mipermu = sample(1:(n1+n2))
mipermu = datos.filt[,mipermu]
mean1 = rowMeans(mipermu[,1:n1])
mean2 = rowMeans(mipermu[,(n1+1):(n1+n2)])
sd1 = apply(mipermu[,1:n1], 1, sd)
sd2 = apply(mipermu[,(n1+1):(n1+n2)], 1, sd)
myparam = list("n" = c(n1,n2), "sd" = cbind(sd1,sd2))
MDperm <- MDbio(dat = cbind(mean1, mean2), param = myparam, a0per = a0per)
Zr = cbind(Zr, myDfunction(mydif = MDperm$D, myrat = MDperm$M, stat = 1, coef = 0.5))
}
}
#-------------------------------------------------------------------------#
## Z-score for different experimental conditions (SIGNAL)
mean1 = rowMeans(as.matrix(datos1.filt))
mean2 = rowMeans(as.matrix(datos2.filt))
sd1 = apply(as.matrix(datos1.filt), 1, sd)
sd2 = apply(as.matrix(datos2.filt), 1, sd)
myparam = list("n" = c(n1,n2), "sd" = cbind(sd1,sd2))
MDs <- MDbio(dat = cbind(mean1, mean2), param = myparam, a0per = a0per)
Zs = myDfunction(mydif = MDs$D, myrat = MDs$M, stat = 1, coef = 0.5)
#-------------------------------------------------------------------------#
## Completing M and D (in signal)
lev1 <- mean1[todos]
lev2 <- mean2[todos]
names(lev1) <- names(lev2) <- todos
Zs <- as.numeric(Zs)
names(Zs) <- rownames(MDs$M)
Zs <- Zs[todos]
names(Zs) <- todos
## Computing Zn
Zn = as.numeric(Zr)
#-------------------------------------------------------------------------#
## Results
list("k" = k, "comp" = comparison, "Level1" = lev1, "Level2" = lev2, "Zs" = Zs, "Zn" = Zn)
}
##############################################################################
##############################################################################
## Function to summarize difference and ratio information (D and D0)
myDfunction <- function (mydif, myrat, stat, coef) {
if (stat == 1) { # linear combination of difference and ratio
myDvalues = coef*mydif + (1-coef)*myrat
}
if (stat == 2) { # distance to origin from (ratio, difference)
myDvalues = sign(mydif) * sqrt((mydif)^2 + (myrat)^2)
}
myDvalues
}
#######################################################################
SoniaTC/NOISeq documentation built on July 28, 2020, 6:31 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
allMD <- function (input, factor, conditions, k = 0.5, replicates, norm = "rpkm",
pnr = 0.2, nss = 5, v = 0.02, lc = 0)
# input: Set of data of type Input
# conditions: Levels of the factor to be compared (when the factor has more than 2 levels)
# k: When counts = 0, 0 will be changed to k. By default, k = 0.5.
# norm: Normalization method. It can be one of "rpkm" (default), "uqua"
# (upper quartile), "tmm" (trimmed mean of M) or "n" (no normalization).
# pnr: Percentage of total reads (seq.depth) for each simulated sample.
# Only needed when noise = "simul". By default, pnr = 1.
# nss: Number of simulated samples (>= 2). By default, nss = 5.
# If nss = 0, real samples are used to compute noise.
# v: Variability in sample total reads used to simulate samples.
# By default, v = 0.02. Sample total reads is computed as a
# random value from a uniform distribution in the interval
# [(pnr-v)*sum(counts), (pnr+v)*sum(counts)]
# lc: Length correction in done by dividing expression by length^lc.
# By default, lc = 1.
{
# n1 <- ncol(as.matrix(datos1))
# n2 <- ncol(as.matrix(datos2))
# g.sinL <- names(which(is.na(long)))
# Check if the factor introduced is already defined
# If the factor introduced is defined and has more than 2 conditions, it will check if the conditions specified are defined too
condition_fac = FALSE
condition_lev = FALSE
datos1 <- datos2 <- matrix()
for (i in colnames(pData(input))) {
if (factor == i) {
condition_fac = TRUE
if (!is.factor(pData(input)[,i])) pData(input)[,i] = as.factor(pData(input)[,i])
if (length(levels(pData(input)[,i])) == 2) {
if (!is.null(assayData(input)$exprs)) {
datos1 <- assayData(input)$exprs[,which(pData(input)[,i] ==levels(pData(input)[,i])[1]), drop = FALSE]
datos2 <- assayData(input)$exprs[,which(pData(input)[,i] ==levels(pData(input)[,i])[2]), drop = FALSE]
} else {
datos1 <- assayData(input)$counts[,which(pData(input)[,i] ==levels(pData(input)[,i])[1]), drop = FALSE]
datos2 <- assayData(input)$counts[,which(pData(input)[,i] ==levels(pData(input)[,i])[2]), drop = FALSE]
}
# Define the comparison string
comparison <- paste(levels(pData(input)[,i])[1], levels(pData(input)[,i])[2], sep=" - ")
condition_lev = TRUE
}
else {
if (is.null(conditions))
stop("Error. You must specify which conditions you wish to compare when the factor has two or more conditions.\n")
if (length(conditions) != 2)
stop("Error. The argument conditions must contain the 2 conditions you wish to compare.")
l <- conditions %in% pData(input)[,i]
# If they are defined, they will be TRUE
if (l[1] == TRUE && l[2] == TRUE) {
if (!is.null(assayData(input)$exprs)) {
datos1 <- assayData(input)$exprs[,which(pData(input)[,i] == conditions[1]), drop = FALSE]
datos2 <- assayData(input)$exprs[,which(pData(input)[,i] == conditions[2]), drop = FALSE]
} else {
datos1 <- assayData(input)$counts[,which(pData(input)[,i] == conditions[1]), drop = FALSE]
datos2 <- assayData(input)$counts[,which(pData(input)[,i] == conditions[2]), drop = FALSE]
}
# Define the comparison string
comparison <- paste(conditions[1],conditions[2], sep=" - ")
condition_lev = TRUE
}
}
}
}
if (condition_fac == FALSE)
stop("The factor you have written does not correspond with any of the ones you have defined.")
if (condition_lev == FALSE)
stop("The conditions you have written don't exist in the factor specified.\n")
# Correction to make it work when there are simulated samples
if (replicates == "no")
replicates = "technical"
# if (description(input)@samples[[1]] == "no")
# description(input)@samples[[1]] = "technical"
n1 <- ncol(as.matrix(datos1))
n2 <- ncol(as.matrix(datos2))
if (norm == "n") { # no normalization
datos1 <- round(datos1, 100)
datos2 <- round(datos2, 100)
}
if (is.null(k)) {
m1 <- min(datos1[noceros(datos1, num = FALSE)], na.rm = TRUE)
m2 <- min(datos2[noceros(datos2, num = FALSE)], na.rm = TRUE)
mm <- min(m1, m2)
k <- mm/2
}
# Total counts for each gene:
suma1 <- rowSums(as.matrix(datos1))
suma2 <- rowSums(as.matrix(datos2))
# All genes
todos <- rownames(as.matrix(datos1))
# Genes with counts in any condition
concounts <- names(which(suma1+suma2 > 0))
long <- 1000
g.sinL <- NULL
if (!is.null(featureData(input)@data$Length)) {
g.sinL <- names(which(is.na(featureData(input)@data$Length)))
if (any(!is.na(featureData(input)@data$Length)) == TRUE)
long <- featureData(input)@data[concounts, "Length"]
}
if (replicates == "technical") { ### technical replicates
suma1 <- suma1[concounts]
suma2 <- suma2[concounts]
#-------------------------------------------------------------------------#
# Normalization of counts for each condition (aggregating replicates)
if (norm == "rpkm") { # RPKM
suma1.norm <- rpkm(suma1, long = long, k = k, lc = lc)
suma2.norm <- rpkm(suma2, long = long, k = k, lc = lc)
}
if (norm == "uqua") {
suma.norm <- uqua(cbind(suma1, suma2), long = long, lc = lc, k = k)
suma1.norm <- as.matrix(suma.norm[ ,1])
suma2.norm <- as.matrix(suma.norm[ ,2])
}
if (norm == "tmm") {
suma.norm <- tmm(as.matrix(cbind(suma1, suma2)), long = long,
lc = lc, k = k)
suma1.norm <- as.matrix(suma.norm[ ,1])
suma2.norm <- as.matrix(suma.norm[ ,2])
}
}
#-------------------------------------------------------------------------#
## Noise distribution
if ((n1+n2)>2) { # with real samples
datitos <- cbind(datos1, datos2)
datitos <- datitos[concounts,]
gens.sin0 <- setdiff(concounts, g.sinL)
if (norm == "n") { # no normalization
datitos.0 <- sinceros(datitos, k = k)
datitos.norm <- datitos.0[gens.sin0, ]
}
if (norm == "rpkm") { # RPKM
datitos.0 <- rpkm(datitos, long = long, k = k, lc = lc)
datitos.norm <- datitos.0[gens.sin0, ]
}
if (norm == "uqua") { # Upper Quartile
datitos.0 <- uqua(datitos, long = long, lc = lc, k = k)
datitos.norm <- datitos.0[gens.sin0, ]
}
if (norm == "tmm") {
datitos.0 <- tmm(datitos, long = long, lc = lc, k = k)
datitos.norm <- datitos.0[gens.sin0, ]
}
datos1.norm <- datitos.norm[ ,1:n1]
datos2.norm <- datitos.norm[ ,(n1+1):(n1+n2)]
if (n1 > 1) {
MD1 <- MD(dat = datos1.norm)
} else { MD1 <- NULL }
if (n2 > 1) {
MD2 <- MD(dat = datos2.norm)
} else { MD2 <- NULL }
} else { # with simulated samples
if (nss == 0) {
nss <- 5
}
datos.sim <- sim.samples(counts1 = sinceros(suma1, k = k),
counts2 = sinceros(suma2, k = k),
pnr = pnr, nss = nss, v = v)
nn <- sapply(datos.sim, ncol)
dat.sim.norm <- vector("list", length = 2)
datitos <- cbind(datos.sim[[1]], datos.sim[[2]])
rownames(datitos) = names(suma1)
sumita <- rowSums(datitos)
g.sin0 <- names(which(sumita > 0))
gens.sin0 <- setdiff(g.sin0, g.sinL)
if (norm == "n") { # no normalization
datitos.0 <- sinceros(datitos, k = k)
datitos.norm <- datitos.0[gens.sin0, ]
}
if (norm == "rpkm") { # RPKM
datitos.0 <- rpkm(datitos, long = long, k = k, lc = lc)
datitos.norm <- datitos.0[gens.sin0, ]
}
if (norm == "uqua") { # Upper Quartile
datitos.0 <- uqua(datitos, long = long, lc = lc, k = k)
datitos.norm <- datitos.0[gens.sin0, ]
}
if (norm == "tmm") {
datitos.0 <- tmm(datitos, long = long, lc = lc, k = k)
datitos.norm <- datitos.0[gens.sin0, ]
}
dat.sim.norm[[1]] <- datitos.norm[ ,1:nn[1]]
dat.sim.norm[[2]] <- datitos.norm[ ,(nn[1]+1):sum(nn)]
MD1 <- MD(dat = dat.sim.norm[[1]])
MD2 <- MD(dat = dat.sim.norm[[2]])
}
Mr <- c(as.numeric(MD1$M), as.numeric(MD2$M))
Dr <- c(as.numeric(MD1$D), as.numeric(MD2$D))
#-------------------------------------------------------------------------#
## M and D for different experimental conditions
if (replicates == "technical" & norm != "n") {
MDs <- MD(dat = cbind(suma1.norm, suma2.norm))
lev1 <- suma1.norm[,1]
lev1 <- lev1[todos]
lev2 <- suma2.norm[,1]
lev2 <- lev2[todos]
} else {
if ((n1+n1) == 2) {
datos1.norm <- sinceros(as.matrix(datos1)[concounts,], k = k)
datos2.norm <- sinceros(as.matrix(datos2)[concounts,], k = k)
}
resum1.norm <- rowMeans(as.matrix(datos1.norm))
resum2.norm <- rowMeans(as.matrix(datos2.norm))
lev1 <- resum1.norm[todos]
lev2 <- resum2.norm[todos]
MDs <- MD(dat = cbind(resum1.norm, resum2.norm))
}
## Completing M and D
names(lev1) <- names(lev2) <- todos
Ms <- as.numeric(MDs$M)
names(Ms) <- rownames(MDs$M)
Ms <- Ms[todos]
names(Ms) <- todos
Ds <- as.numeric(MDs$D)
names(Ds) <- rownames(MDs$D)
Ds <- Ds[todos]
names(Ds) <- todos
## Results
list("k" = k, "comp" = comparison, "Level1" = lev1, "Level2" = lev2, "Ms" = Ms, "Ds" = Ds, "Mn" = Mr, "Dn" = Dr)
}
#######################################################################
#######################################################################
allMDbio = function (input, factor, conditions, k = 0.5, norm = "rpkm", lc = 1,
r = 10, a0per = 0.9, nclust = 15, filter = 1, depth = NULL,
cv.cutoff = 0, cpm = 1)
# input: Set of data of type Input
# conditions: Levels of the factor to be compared (when the factor has more than 2 levels)
# k: When counts = 0, 0 will be changed to k. By default, k = 0.5.
# norm: Normalization method. It can be one of "rpkm" (default), "uqua"
# (upper quartile), "tmm" (trimmed mean of M) or "n" (no normalization).
# lc: Length correction in done by dividing expression by length^lc.
# By default, lc = 1.
# r: Number of permutations to compute null distribution (r=10).
# a0per: Percentile of S to compute a0. If NULL, a0 = 0. (a0per = 0.9)
{
# Check if the factor introduced is already defined
# If the factor introduced is defined and has more than 2 conditions,
# it will check if the conditions specified are defined too
condition_fac = FALSE
condition_lev = FALSE
datos1 <- datos2 <- matrix()
for (i in colnames(pData(input))) {
if (factor == i) {
condition_fac = TRUE
if (!is.factor(pData(input)[,i])) pData(input)[,i] = as.factor(pData(input)[,i])
if (length(levels(pData(input)[,i])) == 2) {
if (!is.null(assayData(input)$exprs)) {
datos1 <- assayData(input)$exprs[,which(pData(input)[,i] ==levels(pData(input)[,i])[1])]
datos2 <- assayData(input)$exprs[,which(pData(input)[,i] ==levels(pData(input)[,i])[2])]
} else {
datos1 <- assayData(input)$counts[,which(pData(input)[,i] ==levels(pData(input)[,i])[1])]
datos2 <- assayData(input)$counts[,which(pData(input)[,i] ==levels(pData(input)[,i])[2])]
}
# Define the comparison string
comparison <- paste(levels(pData(input)[,i])[1], levels(pData(input)[,i])[2], sep=" - ")
condition_lev = TRUE
if (!((ncol(datos1) > 1) && (ncol(datos2) > 1)))
stop("Error. NOISeqBIO needs at least 2 biological replicates per condition.\n")
}
else {
if (is.null(conditions))
stop("Error. You must specify which conditions you wish to compare when the factor has two or more conditions.\n")
if (length(conditions) != 2)
stop("Error. The argument conditions must contain the 2 conditions you wish to compare.")
l <- conditions %in% pData(input)[,i]
# If they are defined, they will be TRUE
if (l[1] == TRUE && l[2] == TRUE) {
if (!is.null(assayData(input)$exprs)) {
datos1 <- assayData(input)$exprs[,which(pData(input)[,i] == conditions[1])]
datos2 <- assayData(input)$exprs[,which(pData(input)[,i] == conditions[2])]
} else {
datos1 <- assayData(input)$counts[,which(pData(input)[,i] == conditions[1])]
datos2 <- assayData(input)$counts[,which(pData(input)[,i] == conditions[2])]
}
# Define the comparison string
comparison <- paste(conditions[1],conditions[2], sep=" - ")
condition_lev = TRUE
}
}
}
}
if (condition_fac == FALSE)
stop("The factor specified does not correspond with any of the ones you have defined.")
if (condition_lev == FALSE)
stop("The conditions specified don't exist for the factor specified.\n")
##-------------------------------------------------------------------------##
# Number of observations within each condition
n1 <- ncol(as.matrix(datos1))
n2 <- ncol(as.matrix(datos2))
if (max(n1,n2) == 1)
stop("There is only one replicate per condition. Please, use NOISeq instead of NOISeqBIO.\n")
# Rounding off data
if (norm == "n") { # no normalization
datos1 <- round(datos1, 10)
datos2 <- round(datos2, 10)
}
# Computing k
if (is.null(k)) {
m1 <- min(datos1[noceros(datos1, num = FALSE)], na.rm = TRUE)
m2 <- min(datos2[noceros(datos2, num = FALSE)], na.rm = TRUE)
k <- min(m1, m2)/2
}
# Total counts for each gene:
suma1 <- rowSums(as.matrix(datos1))
suma2 <- rowSums(as.matrix(datos2))
# Genes with counts in any condition
concounts <- names(which(suma1+suma2 > 0))
# All genes
todos <- rownames(as.matrix(datos1))
# Gene length
long <- 1000
g.sinL <- NULL # genes with no length defined
if (!is.null(featureData(input)@data$Length)) {
g.sinL <- names(which(is.na(featureData(input)@data$Length)))
if (any(!is.na(featureData(input)@data$Length)) == TRUE)
long <- featureData(input)@data[concounts, "Length"]
}
# Genes with counts and with length
gens.sin0 <- setdiff(concounts, g.sinL)
# cond1 and cond2 in the same matrix
datitos <- cbind(datos1, datos2)
datitos <- datitos[concounts,] # selecting only genes with counts
# Sequencing depth when filtering method = 3
if (filter == 3 && is.null(depth)) depth = colSums(datitos)
#-------------------------------------------------------------------------#
#-------------------------------------------------------------------------#
## Normalization
if (norm == "n") { # no normalization
datitos.0 <- sinceros(datitos, k = k)
datitos.norm <- datitos.0[gens.sin0, ]
}
if (norm == "rpkm") { # RPKM
datitos.0 <- rpkm(datitos, long = long, k = k, lc = lc)
datitos.norm <- datitos.0[gens.sin0, ]
}
if (norm == "uqua") { # Upper Quartile
datitos.0 <- uqua(datitos, long = long, lc = lc, k = k)
datitos.norm <- datitos.0[gens.sin0, ]
}
if (norm == "tmm") {
datitos.0 <- tmm(datitos, long = long, lc = lc, k = k)
datitos.norm <- datitos.0[gens.sin0, ]
}
#-------------------------------------------------------------------------#
## Filtering out low count features
if (filter != 0) {
datos.filt = filtered.data(dataset = datitos.norm, factor = c(rep("cond1", n1), rep("cond2", n2)),
norm = TRUE, depth = depth, method = filter, cv.cutoff = cv.cutoff, cpm = cpm)
} else {
datos.filt = datitos.norm
}
datos1.filt <- datos.filt[ ,1:n1]
datos2.filt <- datos.filt[ ,(n1+1):(n1+n2)]
#-------------------------------------------------------------------------#
## Noise distribution
Zr = NULL
if (n1+n2 <= 9) { # sharing information within clusters
Zr = share.info(mydata = datos.filt, n1 = n1, n2 = n2, r = r, nclust = nclust)
} else { # r permutations
for (i in 1:r) {
print(paste("r =", i))
mipermu = sample(1:(n1+n2))
mipermu = datos.filt[,mipermu]
mean1 = rowMeans(mipermu[,1:n1])
mean2 = rowMeans(mipermu[,(n1+1):(n1+n2)])
sd1 = apply(mipermu[,1:n1], 1, sd)
sd2 = apply(mipermu[,(n1+1):(n1+n2)], 1, sd)
myparam = list("n" = c(n1,n2), "sd" = cbind(sd1,sd2))
MDperm <- MDbio(dat = cbind(mean1, mean2), param = myparam, a0per = a0per)
Zr = cbind(Zr, myDfunction(mydif = MDperm$D, myrat = MDperm$M, stat = 1, coef = 0.5))
}
}
#-------------------------------------------------------------------------#
## Z-score for different experimental conditions (SIGNAL)
mean1 = rowMeans(as.matrix(datos1.filt))
mean2 = rowMeans(as.matrix(datos2.filt))
sd1 = apply(as.matrix(datos1.filt), 1, sd)
sd2 = apply(as.matrix(datos2.filt), 1, sd)
myparam = list("n" = c(n1,n2), "sd" = cbind(sd1,sd2))
MDs <- MDbio(dat = cbind(mean1, mean2), param = myparam, a0per = a0per)
Zs = myDfunction(mydif = MDs$D, myrat = MDs$M, stat = 1, coef = 0.5)
#-------------------------------------------------------------------------#
## Completing M and D (in signal)
lev1 <- mean1[todos]
lev2 <- mean2[todos]
names(lev1) <- names(lev2) <- todos
Zs <- as.numeric(Zs)
names(Zs) <- rownames(MDs$M)
Zs <- Zs[todos]
names(Zs) <- todos
## Computing Zn
Zn = as.numeric(Zr)
#-------------------------------------------------------------------------#
## Results
list("k" = k, "comp" = comparison, "Level1" = lev1, "Level2" = lev2, "Zs" = Zs, "Zn" = Zn)
}
##############################################################################
##############################################################################
## Function to summarize difference and ratio information (D and D0)
myDfunction <- function (mydif, myrat, stat, coef) {
if (stat == 1) { # linear combination of difference and ratio
myDvalues = coef*mydif + (1-coef)*myrat
}
if (stat == 2) { # distance to origin from (ratio, difference)
myDvalues = sign(mydif) * sqrt((mydif)^2 + (myrat)^2)
}
myDvalues
}
#######################################################################
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.