14

The simplest manner is to not use a wald test, but rather an LRT with a reduced model lacking the factor of interest: dds = DESeq(dds, test="LRT" reduced=~geno+geno:Treatment) The above would give you results for Treatment regardless of level while still accounting for a possible interaction (i.e., a "main effect of treatment, regardless of the type of ...


9

For normal RNA-seq PCR duplicates are normally kept in, but the duplication rate can be used as a quality control: The higher the duplication rate, the lower the quality. For expression analysis, it is probably best to discard high duplication rate samples, rather than deduplicate them. In general, the smaller the amount of RNA input into the library ...


9

FDR stands for False Discovery Rate. It is a statistic tool used in multiple hypothesis testing. As you may know, when you use a p-value cutoff (usually 0.05) for your experiments, it means strictly speaking that "is there was actually no signal, there would be a probability of 0.05 to observe this kind of extreme values". This can be understood as "there'...


8

You have scale="col" in the code. What you are plotting is the z-score calculated based on the value distribution per column. Try changing it to "row" or "none".


8

It seems that the "combining factors" trick described in part 3.3 of DESeq2 current "vignette" (as of may 2017) under the title "Interaction" is a way to access to the desired contrasts. It seems possible do do it directly when building the colData and when calling DESeqDataSetFromMatrix: Let's add a combined "geno" and "treat" factors to the future ...


7

Mike Love is right. If the response you are looking for is linear in terms of change per minute, the most productive approach is likely to be fitting a linear model. You might get something out of this because the difference between successive time points represents multiple measurements of the rate of change over time. The biggest problem is that the time ...


7

Although it is not recommended to use no replicates, in the edgeR manual they give some advice on how to go on with no replicates design. See page 21 of user guide. You can e.g. estimate a BCV value. Of course it is still a trick and not sound statistics.


6

Yes, that blog post does represent just one guy's opinion (hi!) and it does date all the way back to 2014, which is, like, decades in genomics years. :-) By the way, there is quite a bit of literature discussing the improvements that expected read counts derived from an Expectation Maximization algorithm provide over raw read counts. I'd suggest reading the ...


5

This happen when the variables (strain +batch) create a design matrix like this: batch strain 1 1 # 1 1 # 1 2 2 2 3 3 4 3 ... 16 72 Which means that some of the covariates are not linearly independent (ie batch 1 and strain 1), all the strain 1 is in batch 1. You can correct for batch effects, but not in this design of the linear model (if you want to ...


5

Try the Gene Expression Omnibus - it looks like they have some datasets.


5

I don't think that the issue is the low counts, but rather the number of features without any real variance (the black dots at the bottom). So what the heck is the dispersion plot and why does one need to fit it anyway? In a typical RNAseq experiment, one measures many thousands of genes with only a few replicates per biological group. This then leads to ...


5

I would be very hesitant to blindly exclude those samples based on the clustering. Check to see if the clusters actually denote some sort of batch effect, since it's not like all of the TCGA datasets were processed at the same place or time. Check for clinical covariates too. Also, do a PCA and check the genes with high loading on the relevant PC. If it ...


5

In the linked article the authors formalize microarray analysis as the study of the joint distributions of $\overrightarrow{X}_i$ and $Y_i$, where $\overrightarrow{X}_i$ is a vector of random variables, the distribution of each of which (i.e. each of $X_{ji}$) is determined by the level of expression of gene j in sample i, and $Y$ is some response or ...


4

Remove low-count features in advance. This is standard for most tools including DESeq2 and edgeR (see section 2.6). This will keep you from testing a lot of features that cannot be differentially expressed (using NB model you need at least four reads - below that the uncertainty is too high to ever be DE) and you will make the false discovery correction ...


4

If you're interested in looking at the differences between two treatments then you'll end up wanting to do both a direct contrast as well as the individual comparisons to baseline. The direct contrast will give you the genes actually differentially expressed between the two conditions. In practice, you may want to filter this a bit so you only have genes ...


4

I assume that you are talking about the implementation of these methods in the limma package. Otherwise this answer does not apply. I think that your questions can be answered with some simulations where we can test with some "genes" with a known relationship: library("limma") set.seed(123) # Create some genes and samples nGenes <- 40 nSamples <- ...


4

You can simply apply() t.test() to your matrix. In general, though, I expect lmFit() from the limma package will be simpler to use. Here is an example with t.test(): m = matrix(runif(100, max=c(rep(1, 50), rep(1.5, 50))), nrow=10) # make some data apply(m, 1, function(x) t.test(x[1:5], x[6:10])$p.value) Neither t.test() nor lmFit() should be expected to ...


4

The edgeR authors recommend that you use a relatively low logFC threshold for glmTreat such as lfc=log2(1.2). A lfc value as high as lfc=log2(2) is seldom required and only would only be appropriate for datasets with really large amounts of DE. The first threshold you used of lfc=2 is equivalent to lfc=log2(4) and is too high for any dataset. There's a few ...


3

You can specify the exact comparisons you want in the results() function. So: dds = DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = directory, design= ~ condition) dds = DESeq(dds) res = results(dds, contrast("condition", "treatment1", "control1")) Note that condition should not be, "stim001-control1, stim001-control2, etc.", but instead, ...


3

Generally you should just leave them as is. One does remove/mark duplicates in DNA seq. For further read check this Nature paper


3

Include all genes/transcripts in your analysis. A transcript that is not detected could be undetected through sampling error (i.e. the sequencer / library prep just happened to miss that transcript), or it could be because the transcript isn't generated in a particular sample. It's not uncommon for genes to be switched off in response to different ...


3

Depending upon the type of treatment used the set of DEGs will change.If the treatments have similar kind of effect you will get a small list (less variable genes will have higher p-val) using a cutoff of p<0.05. So, it's better to start with control vs treated comparisons then T1 vs T2. Compare the lists form CTRL vs T1 and CTRL vs T2 you will get the ...


3

You're misinterpreting the moderated T-statistic, it's basically the fold-change divided by its variance. The p-value comes directly from that, so if you filter by moderated fold-change you're just setting an unknown (unless you go through the trouble of figuring it out) p-value threshold. Instead, do exactly as you've been doing and set appropriate (...


3

I'm not sure if nowadays matlab is still a good option for informatics. I would rather go for a R script to make a volcano plot. If you don't want to have fancy code but functional one, you can even make the foldchange calculations and log2 transformation in excel and just plot the log2 fold change in x axis and the -log10 transform of the p/q values on the ...


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 ...


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

Fold-change >= 2 is the same as logFC (log2(fold-change)) >= 1, so your example is doing exactly what you want. logFC is generally easier to think in and work with than fold-change, since (aside from computational reasons) then increases and decreases in expression differ only in sign and are, thus, easier to compare in magnitude (e.g., it's easier to tell ...


3

You can usually get away with FDR < 0.1, but that's as high as you can go. This all presumes you're doing follow-up experiments of some sort, of course. I guess you could increase the FDR more, but you're then really increasing the odds that your follow-up experiments will fail. Obviously increasing the FDR will decrease the confidence in the results, it'...


3

I would suggest using a likelihood ratio test for differential expression using logistic regression with batch as a latent variable. In Seurat you can do: markers <- FindAllMarkers(object, test.use = 'LR', latent.vars = 'batch') (change "object" and "batch" accordingly) See https://www.biorxiv.org/content/early/2018/02/14/258566 from Lior Pachter's ...


3

No. If the only difference between Tumour_batch1 and Tumour_batch2 was the library prep/sequencing run, then all this will tell you is if the particular sample of RNA taken from the tumour is different from the particular sample of RNA taken from normal. It almost certainly is. Indeed if you did technical replicates of two RNA samples, both taken from ...


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