Hot answers tagged

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


12

The normalized counts themselves can be accessed with counts(dds, normalized=T). Now as to what the baseMean actually means, that will depend upon whether an "expanded model matrix" is in use or not. Given your previous question, we can see that geno_treat has a bunch of levels, which means that expanded models are not in use. In such cases, the baseMean ...


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


8

You need to specify the number without the version. Instead of "ENSMUST00000178862.1" just "ENSMUST00000178862": You can do this with one more line: g <- gsub("\\..*", "", rownames(txi.kallisto$counts)) (hgnc_symbols <- getBM(attributes = c("hgnc_symbol", "chromosome_name", "ensembl_transcript_id"), filters = "ensembl_transcript_id", values = g, mart ...


8

You only have 4 samples total. I think it would be difficult to not have the PCA show big differences between the groups with so few points. On the other hand, for differential expression, it is hard to get something to be statistically significant with only 2 replicates.


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.


7

The counts are "reads" for single-end datasets and "fragments" for paired-end datasets. In other words, they're what featureCounts and htseq-count produce. The "normalized counts" that you'll be able to access in DESeq2 are, as aptly named, simply normalized versions such that values are comparable across samples.


7

You originally had asked a very broad question, so I'll try to demonstrate why that is such a hard question to answer. I've done two fairly large differential analysis studies (and a few smaller ones) covering very different areas of research, and the approaches that other researchers used subsequent to my differential expression calculations were ...


6

It depends what you mean by “normalised”. As Devon said, the normalized = TRUE argument to the count function gives you normalised counts. However, these are “only” library-size normalised (i.e. divided by the sizeFactors(dds)). However, as the vignette explains, downstream processing generally requires more advanced normalisation, to account for the ...


6

There are a variety of reasons people use gene-level quantitations. Transcript-level differences are difficult to biologically interpret. Let's be honest, few groups are likely to put in the work required to determine what these might mean. Most genes have at least some level of characterization, so it's much more biologically tractable to think in terms of ...


6

Ultimately my colleague helped me to solve the issue by following the steps: Created environment: conda create --name myenv Activated it: source activate myenv Edited .condarc to set the priority of the channels to be the following (top to bottom): conda-forge, bioconda, r, default Checked that the channels are set in the correct order: conda config --get ...


6

You can't (easily) use a single design, since clone is nested within diabetic status, so you were correct using two separate designs. Your designs are correct. If you really wanted to use a single design for some reason, it'd be ~clone+structure and for the latter comparison you'd end up using a rather complicated contrast (something like c(rep(c(1,-1), 3), ...


6

Transcript abundance quantification is a tricky topic since a read often could belong to several transcripts, so any "count" is a best guess as to which transcript it actually originates from. That being said, there are tools that can help you here: salmon (as you mentioned) to quanitfy. Run it with --numGibbsSamples 50 (or higher if your computer ...


5

If they're truly technical replicates, then there's no way to model them using DESeq2*, as you've alluded to with the collapseReplicates function. DESeq2/ Mike Love's general recommendation with collapseReplicates is to just add the reads together for technical replicates. If you want to model them instead of collapsing them down, you can voom transform ...


5

Practically speaking, there's no way to include the technical replicates in that design (in DESeq2 at least). Your concern regarding inflating the power is exactly correct and the only way to combat that would be to add a pairing factor like one might do with case-control or tumor-normal studies. That is, something like: group libraryPrep sample 1 WT ...


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

You are right, samples.txt is not generated by Salmon (or any other transcript abundance quantifiers). From documentation, you can find a link to an example of how samples.txt should look like. Note: At the moment, the server is unavailable (I got an Error 503: Service Temporarily Unavailable) - I'd suggest you get directly in contact with the owners / ...


5

I think the problem is that you have time as a linear value rather than a factor. While this naively makes sense (after all, you have to go through 3 hours awake before you can get to 9 hours awake), what you end up doing is fitting a single coefficient to each gene of "linearly changes with time". What you want instead is factor, so you can have a system ...


4

The first method (resSig$log2FoldChange) is for getting the most differentially expressed genes. These are the genes with the biggest differences between specific groups you pre-define. The second method (rowVars(assay(rld))) is for getting the highly variable genes. These are the genes most variable across all samples regardless of which samples they are.


4

It's likely that the fold-change of the topGene is greater than 5. Consequently, the circle is being drawn and the text produced, but they're outside the bounds of the plot (the bounds don't update if you add a circle and text label). Look at the log2FoldChange column of resLFC1[topGene, ] and increase the bounds to encompass that.


4

Any IDs (e.g. replicate number) don't need to be included in the design formula if it's not a repeated factor in multiple samples. In this case, you have said that mouse number is related to the number of mice used in the pool (and is a variable factor that is consistent across multiple samples), in which case it should be included in the design matrix. ...


4

See the second example from the vignette, which is exactly your situation. The software cannot differentiate between the batch effects vs. condition effects because the groups are exactly the same.


4

Provide rank sufficient design to DESeqDataSetFromMatrix and then use your custom model matrix in DESeq. In essence: dds = DESeqDataSetFromMatric(counts, s2c, design=~batch) design <- model.matrix(~strain+batch, s2c) design = design[, -9] DESeq(dds, full=design) See this thread on the bioconductor site for details.


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

To add to the list that Devon Ryan outlined (or perhaps to elaborate on point 2?): Although Salmon/Kallisto/RSEM are the more accurate in their transcript quantification than the methods they superseded, transcript level quantification is still not as accurate as gene level quantification (see the tximport paper, which would also be the tool i'd recommend ...


4

First I installed this libraries: sudo apt-get install libcurl4-openssl-dev libxml2-dev I then installed RCurl: BiocManager::install("RCurl") and lastly: BiocManager::install("DESeq2", version = "3.8")


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


4

In the DESeq2 manual there is a section titled "How can I include a continuous covariate in the design formula?" that deals with your question. Basically the process is no different from using a discrete covariate. DESeq <- DESeqDataSetFromMatrix(countData = counts, colData = metadata, design = ~ RIN + condition) In this example condition is the ...


4

From what I can gather, you want to account for the effect of individual, nested within region. That is, you want to see after accounting for these, is there a consistent effect for Injection:Social across all conditions. So you set up the model like this: m1 <- model.matrix(~ ind.n*Region + Injection + Social + Injection:Social,data=..) The last term ...


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