Hot answers tagged

5

It is indeed true that for the DE analysis one should include batch into the formula to avoid changing the original counts. Still, for everything else such as plotting heatmaps use of removeBatchEffects is perfectly fine and (at least to me) a standard and well-accepted procedure. It essentially does not matter what you use to correct for the batch effect ...


4

Not really an answer but an extended comment... and most likely something you don't like to hear I guess by technical replicates, it means taking the same biological sample and making 2 methylation libraries. If this is used as replicates in deseq2 or edger, the variance you are estimating is the technical variation that comes with preparing the library, ...


3

It looks good to me. One think I could suggest is to see the distribution of your data/model (with limma/voom). After you get the list of your genes with: top.table <- topTable(fit.eBayes, sort.by = "P", n = Inf) Maybe you can check your multidimensional scaling (MDS): plotMDS(your_data, col = as.numeric(group)) Then you can also check with ...


3

For comparing the counts of different samples from DESeq2, Michael Love recommends using the variance-stabilized transform. It'd be great if you could provide some specific code examples in your question, but without that here's something that should work with the DESeq2 workflow as mentioned in the package documentation: ## [assuming "dds <- DESeq(dds)" ...


3

While I share Llopis' concern about estimating variance from 2 samples, the statement you quoted is about avoiding false positives from genes that are only expressed in a few samples. It's fairly common to exclude genes that expressed in fewer samples than the smallest group even if the number of samples per group is much higher than 3.


3

The problem with most of the methods is that use the gene's variance for each group, which can't be calculated (reliable) when the sample is <= 2. Also, statistically, it would have extremely low power, so the conclusions couldn't be trusted much. You could also calculate the "raw" fold change by yourself (ie without the variance estimation and ...


2

If you just want the change with time then use time as a continuous variable. The log2FC will then be the change per hour (you can use the default Wald test to extract it). For time-point specific comparisons you'll instead need to treat the time points as groups and use contrasts as appropriate.


2

The mroast function has an argument to specify which contrast do you want to test, quoting from the help page: contrast contrast for which the test is required. Can be an integer specifying a column of design, or the name of a column of design, or a numeric contrast vector of length equal to the number of columns of design. So if you have the contrast ...


2

You can use either, but lmFit has the benefit of returning an object that can be used with eBayes() so you can pool information across genes/probes/whatever. lm() is a base R function applicable basically everywhere. lmFit() if from the limma package, so originally intended for microarray data, though these days pretty much everything omics is analyzed with ...


2

To add to what @gringer, when you do use TPM, the normalization done is for both library size and gene length. When you use rlog, the normalization is done via median normalization (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2864565/). See below for a quick example about the rlog, and size factor normalization. dds=makeExampleDESeqDataSet(betaSD=0.5,...


2

You can use the following code to calculate the coefficient of variation: # expr is your expression matrix. SD <- apply(expr, 1, sd) CV <- sqrt(exp(SD^2) - 1) It might be implemented in some package but it is so brief that you can write it again yourself. Then you can filter out those that are below certain percentage of the distribution of CV (like ...


2

You're basically subtracting a constant per-gene per-level. The relevant portion of the code is: fit <- lmFit(x,cbind(design,X.batch),...) beta <- fit$coefficients[,-(1:ncol(design)),drop=FALSE] beta[is.na(beta)] <- 0 as.matrix(x) - beta %*% t(X.batch) where x is your input matrix, design is the design matrix and X.batch is the matrix of batch and ...


2

August 2021 To date, it seems that the response to this question is "No". However, there is a GitHub repository called edgePy, aiming to "become an implementation of edgeR for differential expression analysis in the Python language" (edgeR is somewhat similar to limma). PS: Many thanks to @ATpoint90 (the Twitter "hive mind" ...


2

The error is telling you that there isn't any Infected column on the model.matrix. Check the column names of mm. Also not sure that having two ~ will work well.


1

You may run edgeR for methylation analysis without replicates (https://f1000research.com/articles/6-2055). But I also recommend you to have a look at the R DSS package (http://bioconductor.org/packages/release/bioc/html/DSS.html). It has a smoothing step which allows to modelize the biological variability by taking information from neighboring sites.


1

That is not trivial. GEO2R uses the Bayesian linear model-based framework limma https://bioconductor.org/packages/release/bioc/html/limma.html for the analysis. I doubt you can easily and meaningfully implement the relevant code in Python without a major effort and notable background knowledge. The user guide https://bioconductor.org/packages/release/bioc/...


1

DESeq2 uses the batch information (and everything else in the design) to produce offsets for its GLM. For a background on that please check how linear models work, e.g. using the StatQuest series of statistics videos over at YouTube. It still operates on the raw counts. The same goes for the normalization factors. removeBatchEffect fits a linear model to the ...


1

If you got DEGs with statistically significant p-values after multiple testing following a valid pipeline, I'd believe that, even if looking at the first few PCs doesn't look promising.


1

Try to use DElist() function before you transform, and also make rownames first. counts <- prac_count_10[,-1] rownames(counts) <- prac_count_10[,1] DGE1 <- DGEList(counts) And then continue with your voom...


1

Treat your disease samples as individual groups and then follow the normal routine in limma to use contrasts to compare two group (the control group versus individual disease samples). Note that you need to think long and hard about what these results then mean. Normally we compare groups because then the results should generalize to other samples. That will ...


1

I would use the p-value/FDR which each method returns to rank the gene list for that method in order from 'most likely to be DE' to least likely. You would then have three ranked lists of genes - you would put the actual rankings in for each gene - 1 for the top ranked, 2 for the next, etc. You could then carry out, say Spearman's rank correlation on the ...


1

You'll have a bit more power if you use all of the samples at once. Nothing is useful for this, don't waste your time with it. You're not looking for outliers, you're looking to see how samples cluster. You can just use correlation for that.


Only top voted, non community-wiki answers of a minimum length are eligible