How can I get normalized counts with ImpulseDE2?
0
0
Entering edit mode
MU • 0
@037ddc51
Last seen 7 months ago
Japan

Hi, I am trying RNA-seq analysis with ImpluseDE2. As a next step, I want to make a plot of normalized counts for each genes I am interested in, or make do co-expression analysis with metabolome data. How can I get the normalized count data for each gene? And, I don't understand how to use "Timecatg" in dfannotation. How should I set this parameter?

Any advice you could give would be much appreciated. I apologize if my English is difficult to understand.

ImpulseDE2 • 197 views
0
Entering edit mode

I also have similar issue. I think the best way to do this is use plotHeatmap function or rather its truncated version. I paste the code below:

objectImpulseDE2 <- objectImpulseDE2 # object, which is an output from runImpulseDE2 function
strCondition  <- "case"
boolIdentifyTransients <- TRUE # TRUE or FALSE, dependent on your analysis, it should work in both cases

#for some unknown reasons I have to define evalImpulse_comp function, which should be already defined

evalImpulse_comp <- function(vecImpulseParam, vecTimepoints)
{
vecImpulseValue <- sapply(vecTimepoints, function(t) {
(1/vecImpulseParam[3]) * (vecImpulseParam[2] + (vecImpulseParam[3] -
vecImpulseParam[2]) * (1/(1 + exp(-vecImpulseParam[1] *
(t - vecImpulseParam[5]))))) * (vecImpulseParam[4] +
(vecImpulseParam[3] - vecImpulseParam[4]) * (1/(1 +
exp(vecImpulseParam[1] * (t - vecImpulseParam[6])))))
})
vecImpulseValue[vecImpulseValue < 10^(-10)] <- 10^(-10)
return(vecImpulseValue)
}

#now its time to define our function of interest, which is a simple truncation of plotHeatmap
get_norm_count_table <- function(objectImpulseDE2, strCondition, boolIdentifyTransients, scaQThres)
{
dfAnnot <- get_dfAnnotationProc(obj = objectImpulseDE2)
scaNGenes <- dim(get_matCountDataProc(obj = objectImpulseDE2))[1]
vecSignificantIDs <- rownames(objectImpulseDE2$dfImpulseDE2Results[!is.na(objectImpulseDE2$dfImpulseDE2Results$padj) & objectImpulseDE2$dfImpulseDE2Results$padj < scaQThres, ]) vecTimePointsToEval <- sort(unique(dfAnnot$Time), decreasing = FALSE)
scaNTPtoEvaluate <- length(vecTimePointsToEval)
matImpulseValue <- do.call(rbind, lapply(vecSignificantIDs,
function(x) {
evalImpulse_comp(vecImpulseParam = get_lsModelFits(obj = objectImpulseDE2)[[strCondition]][[x]]$lsImpulseFit$vecImpulseParam,
vecTimepoints = vecTimePointsToEval)
}))
rownames(matImpulseValue) <- vecSignificantIDs
matidxMaxTimeSort <- t(apply(matImpulseValue, 1, function(genevalues) {
sort(genevalues, decreasing = TRUE, index.return = TRUE)$ix })) vecMaxTime <- vecTimePointsToEval[matidxMaxTimeSort[, 1]] matidxMinTimeSort <- t(apply(matImpulseValue, 1, function(genevalues) { sort(genevalues, decreasing = FALSE, index.return = TRUE)$ix
}))
if (boolIdentifyTransients) {
vecidxTransient <- which(objectImpulseDE2$dfImpulseDE2Results[vecSignificantIDs, ]$isTransient)
vecidxMonotonous <- which(objectImpulseDE2$dfImpulseDE2Results[vecSignificantIDs, ]$isMonotonous)
if (length(vecidxMonotonous) == 1) {
matImpulseMonot <- t(matImpulseValue[vecidxMonotonous,
])
}
else {
matImpulseMonot <- matImpulseValue[vecidxMonotonous,
]
}
vecboolMonotonousUp <- apply(matImpulseMonot, 1, function(gene) {
gene[1] < gene[scaNTPtoEvaluate]
})
vecboolMonotonousDown <- !vecboolMonotonousUp
vecidxMonotonousUp <- vecidxMonotonous[vecboolMonotonousUp]
vecidxMonotonousUpSort <- vecidxMonotonousUp[do.call(order,
as.data.frame(matidxMinTimeSort[vecidxMonotonousUp,
1]))]
vecidxMonotonousDown <- vecidxMonotonous[vecboolMonotonousDown]
vecidxMonotonousDownSort <- vecidxMonotonousDown[do.call(order,
as.data.frame(matidxMinTimeSort[vecidxMonotonousDown,
1]))]
if (length(vecidxTransient) == 1) {
matImpulseTransient <- t(matImpulseValue[vecidxTransient,
])
}
else {
matImpulseTransient <- matImpulseValue[vecidxTransient,
]
}
vecboolTransientValley <- apply(matImpulseTransient,
1, function(genevalues) {
boolValley <- any(genevalues[2:(scaNTPtoEvaluate -
1)] < genevalues[1] & genevalues[2:(scaNTPtoEvaluate -
1)] < genevalues[scaNTPtoEvaluate])
return(boolValley)
})
vecboolTransientPeak <- !vecboolTransientValley
vecidxTransientPeak <- vecidxTransient[vecboolTransientPeak]
vecidxTransientPeakSort <- vecidxTransientPeak[do.call(order,
as.data.frame(matidxMaxTimeSort[vecidxTransientPeak,
1]))]
vecidxTransientValley <- vecidxTransient[vecboolTransientValley]
vecidxTransientValleySort <- vecidxTransientValley[do.call(order,
as.data.frame(matidxMinTimeSort[vecidxTransientValley,
1]))]
vecidxAllSort <- c(vecidxMonotonousUpSort, vecidxMonotonousDownSort,
vecidxTransientPeakSort, vecidxTransientValleySort)
vecTrajectoryType <- c(rep("up", length(vecidxMonotonousUpSort)),
rep("down", length(vecidxMonotonousDownSort)), rep("*up",
length(vecidxTransientPeakSort)), rep("*down",
length(vecidxTransientValleySort)))
lsvecGeneGroups <- list(transition_up = vecSignificantIDs[vecidxMonotonousUpSort],
transition_down = vecSignificantIDs[vecidxMonotonousDownSort],
transient_up = vecSignificantIDs[vecidxTransientPeakSort],
transient_down = vecSignificantIDs[vecidxTransientValleySort])
}
else {
vecidxAllSort <- sort(vecMaxTime, decreasing = FALSE,
index.return = TRUE)$ix vecTrajectoryType <- rep(" ", length(vecidxAllSort)) lsvecGeneGroups <- list(all = vecSignificantIDs[vecidxAllSort]) } vecUniqueTP <- unique(dfAnnot$Time)
vecSizeFactors <- computeNormConst(matCountDataProc = get_matCountDataProc(obj = objectImpulseDE2),
vecSizeFactorsExternal = get_vecSizeFactors(obj = objectImpulseDE2))
matSizefactors <- matrix(vecSizeFactors, nrow = length(vecSignificantIDs),
ncol = dim(get_matCountDataProc(obj = objectImpulseDE2))[2],
byrow = TRUE)
matDataNorm <- get_matCountDataProc(obj = objectImpulseDE2)[vecSignificantIDs,
]/matSizefactors
matDataHeat <- do.call(cbind, lapply(vecUniqueTP, function(tp) {
vecidxCols <- which(dfAnnot\$Time %in% tp)
if (length(vecidxCols) > 1) {
return(rowMeans(matDataNorm[, vecidxCols], na.rm = TRUE))
}
else {
return(matDataNorm[, vecidxCols])
}
}))

}

#now you can run your analysis
get_norm_count_table(objectImpulseDE2, strCondition, boolIdentifyTransients, scaQThres)
#your results are in matDataNorm variable
my_norm_counts <-  matDataNorm


You can improve my solution in many ways (e.g. by adding 'return matDataNorm' at the end of function), but I see that you post the issue ~half year ago, so I put here first working solution. Hopefuly it will help.