Hi Noa, On Feb 11, 2014 7:52 AM, "Noa Henig" <ntzunz@gmail.com> wrote: > > Hi Michael, > I am working on an a new setup that follows the design you suggested, should have the results soon. Meanwhile I wanted to get a quick estimate for the differential expression only in one of the cell types but I still have a problem with high fold change for samples having zero or very low counts: > I used a partial matrix that includes only the replicates of one of the cell types, so there is a single factor with 3 possible levels in this design. The code was the same as I sent last Sunday, where betaPrior was set to FALSE, except the counts matrix that included this time only 6 replicates of U_V0, 4 replicates of U_K1 and 2 replicates of U_p5: > > "U_V0", "U_V0","U_V0", "U_V0", "U_V0", "U_V0", "U_K1","U_K1", "U_K1", "U_K1", "U_p5", "U_p5" > > The normalized counts for a certain gene were the following: > 0, 0, 0, 0, 0, 0.83, 0, 1.12, 0, 0, 0, 0 > > The log2 fold change for that gene for the comparison U_V0 vs. U_p5 was -13.38 > For the comparison of U_V0 vs. U_K1 it was 1.01. > > Can you tell what is the problem in this case? See my last email. This is not a problem. U_p5 has zero counts, so expected value zero. The base level has small, positive expected value. The fold change is from a small positive value to zero. So log(0). Note that U_V0 is the base level in your design. So I assume you are actually reporting log fold changes over the base level. Mike > Thanks again, > Noa > > > On Sun, Feb 9, 2014 at 8:27 PM, Michael Love <michaelisaiahlove@gmail.com> wrote: >> >> Sorry, I left out the interaction term in my earlier response, which is important in order to test if the effect of condition is different in the second cell type as in the first cell type. >> >> On Sun, Feb 9, 2014 at 12:31 PM, Michael Love < michaelisaiahlove@gmail.com> wrote: >>> >>> hi Noa, >>> >>> >>> On Sun, Feb 9, 2014 at 10:25 AM, Noa Henig <ntzunz@gmail.com> wrote: >>>> >>>> Hi DESeq developers, >>>> I am running DESeq2 1.2.8 to receive differential expression in 2 different >>>> cell types. For each cell type I have control samples, over- expression of >>>> gene A, and over-expression of gene B, while the second cell type is used >>>> in order to support the findings from the first cell type. >>>> >>>> I have 2 questions: >>>> First, I found some genes that had zero normalized counts in 7 out of 8 >>>> replicates (the last one had 1 count or 1.6 counts). The fold change for >>>> comparison for such two genes was -13.4 and 2.47. Since the 'condition' has >>>> more than 3 levels I set the betaPrior to False as Michael indicated in the >>>> earlier correspondence (the code appears below). >>>> Can you tell the reason for this this and how to correct this? >>> >>> >>> I don't understand what is unexpected here. It looks like you are not following my recommendation in a previous email to use a design formula with celltype and condition, e.g.: >>> >>> ~ celltype + condition >> >> >> >> I would recommend you use for this experiment the design formula: >> >> ~ celltype + condition + celltype:condition >> >> Let's say that celltype has two levels "Cell1" and "Cell2", and condition has three levels "Level1", "Level2", "Level3". Then to examine the effect of one of the levels of condition you would use, e.g.: >> >> results(dds, name="condition_Level2_vs_Level1") >> ...or equivalently: >> results(dds, contrast=c("condition","Level2","Level1")) >> >> You can compare any pairs of levels using the contrast argument. >> >> In order to test if the second cell type had a different effect of condition as the first cell-type, you would use: >> >> results(dds, name="celltypeCell2:conditionLevel2") >> >> ...for the level2 effect and, >> >> results(dds, name="celltypeCell2:conditionLevel3") >> >> ...for the level3 effect. >> >> >> Mike >> >> >>> >>> >>> where celltype would be a column in colData with levels U or M, and condition would be a column in colData with levels V0, K1 or p5. >>> >>> When one of the levels specified in the contrast has zero counts, and we do not use a prior on log fold change, the log fold change estimate will slope toward + or - infinity, but the algorithm will stop if the increase in the likelihood is less than a threshold (as does the glm() function in base R), or if the estimate passes fold change of 1024 or 1/1024. >>> >>> You don't need to correct these estimates, they will either be significant or not by adjusted p-value (as can be seen in the MA-plot) depending on the estimate of dispersion. There is typically a curved boundary region in the MA-plot where fold changes from low count genes are not significant by adjusted p-value, while fold change from high count genes are significant. >>> >>> >>> Mike >>> >>> >>>> >>>> >>>> In addition, I received the following warning: >>>> Warning message: >>>> In parametricDispersionFit(mcols(objectNZ)$baseMean[useForFit], : >>>> the parametric dispersion fit did not converge, >>>> which occurs when the dispersion estimates over the mean do not fit to the >>>> curve y = a/x + b. You should re-run the analysis using fitType='local' >>>> or 'mean' (DESeq2 > v1.2 will re-run automatically). >>>> When using local regression fit, the user should examine plotDispEsts(dds) >>>> to make sure the fitted line is not sharply curving up or down based on >>>> the position of individual points. >>>> >>>> I guess that by running version 1.2.8 the additional fit types were ran >>>> automatically but can I see the actual code that was ran ? >>>> >>> >>> This error message is admittedly a bit confusing, and I've improved this in the devel branch. Better would be to say: fitType="parametric" didn't work. The dispersion estimates do not fit well to the parametrization. So fitType="local" was automatically substituted. You can avoid this warning if you rerun the analysis by specifying fitType="local" or "mean" from the start. >>> >>> >>>> >>>> This is the sessionInfo: >>>> R version 3.0.2 (2013-09-25) >>>> Platform: i386-w64-mingw32/i386 (32-bit) >>>> >>>> locale: >>>> [1] LC_COLLATE=Hebrew_Israel.1255 LC_CTYPE=Hebrew_Israel.1255 >>>> LC_MONETARY=Hebrew_Israel.1255 >>>> [4] LC_NUMERIC=C LC_TIME=Hebrew_Israel.1255 >>>> >>>> attached base packages: >>>> [1] parallel stats graphics grDevices utils datasets methods >>>> base >>>> >>>> other attached packages: >>>> [1] DESeq2_1.2.8 RcppArmadillo_0.4.000 Rcpp_0.10.6 >>>> GenomicRanges_1.14.4 XVector_0.2.0 >>>> [6] IRanges_1.20.6 BiocGenerics_0.8.0 >>>> >>>> loaded via a namespace (and not attached): >>>> [1] annotate_1.40.0 AnnotationDbi_1.24.0 Biobase_2.22.0 >>>> DBI_0.2-7 genefilter_1.44.0 >>>> [6] grid_3.0.2 lattice_0.20-23 locfit_1.5-9.1 >>>> RColorBrewer_1.0-5 RSQLite_0.11.4 >>>> [11] splines_3.0.2 stats4_3.0.2 survival_2.37-4 >>>> XML_3.98-1.1 xtable_1.7-1 >>>> >>>> And this is the code I was running: >>>> >>>> tmpcountData = read.table(countsTable, header=TRUE, row.names=1) >>>> >>>> numOfConds=6 >>>> >>>> samplesList = c("U_V0", "U_V0","U_V0", "U_V0", "U_V0", "U_V0", "U_K1", >>>> "U_K1", "U_K1", "U_K1", "U_p5", "U_p5", "M_V0", "M_V0", "M_V0", >>>> "M_K1","M_K1", "M_p5", "M_p5") >>>> >>>> countData = tmpcountData[,3:21] # counts table contains 2 more columns >>>> before the counts >>>> >>>> colData = data.frame(row.names=colnames(countData), condition = samplesList) >>>> >>>> conditionsList = c("U_V0","U_K1", "U_p5","M_V0","M_K1","M_p5") >>>> >>>> colData$condition = factor( colData$condition, levels = >>>> conditionsList) >>>> >>>> dds <- DESeqDataSetFromMatrix(countData = countData, colData = colData, >>>> design = ~condition) >>>> >>>> dds <- estimateSizeFactors(dds) >>>> >>>> dds <- estimateDispersions(dds, maxit = 100) >>>> >>>> dds <- nbinomWaldTest(dds, betaPrior = FALSE) >>>> >>>> png("plotDispEsts_final.png") >>>> >>>> plotDispEsts(dds) >>>> >>>> dev.off() >>>> >>>> tmp_counts = counts(dds, normalized=TRUE) >>>> >>>> normalized_counts = data.frame(tmp_counts) >>>> >>>> >>>> ... reading the columns of counts for comparisons and writing the results to >>>> a table.... >>>> >>>> >>>> Thanks! >>>> Noa Henig >>>> >>>> >>>> >>>> >>>> >>>> On Tue, Jan 7, 2014 at 4:21 PM, Michael Love < michaelisaiahlove@gmail.com>wrote: >>>> >>>> > Hi Noa, >>>> > >>>> > Thanks for reporting in. This sounds like the same problem that a user >>>> > had in November: >>>> > >>>> > >>>> > http://permalink.gmane.org/gmane.science.biology.informatics.conductor /51749 >>>> > >>>> > >>>> > If you are using the updated release version (>= 1.2.6), there should >>>> > have been a warning printed when you build the DESeqDataSet that you should >>>> > set betaPrior=FALSE because factors are present in the design with 3 or >>>> > more levels. >>>> > >>>> > In the devel branch of DESeq2 (v1.3), we have implemented a solution that >>>> > allows using a prior on log fold changes for factors with 3 or more levels >>>> > results. Then DESeq() or nbinomWaldTest() will provide symmetric results >>>> > regardless of the base level and it takes care of this situation with >>>> > nonzero log fold changes from contrasts in which both conditions have zeros. >>>> > >>>> > With future questions please >>>> > post your full code, the output of sessionInfo() >>>> > , as this helps us provide better answers. >>>> > >>>> > Mike >>>> > On Jan 7, 2014 8:49 AM, "Noa Henig" <ntzunz@gmail.com> wrote: >>>> > >>>> >> Hi, >>>> >> I am using DESeq2 to get differential expression of 45 samples which >>>> >> represent 26 different conditions. >>>> >> I read that the calculation of the fold change in changed DESeq2 and >>>> >> affects the genes having low counts more than the highly expressed genes. >>>> >> However, I found genes that had zero counts in series of 4 samples (2 >>>> >> replicates in each of 2 conditions), while their log2fold change was 3. >>>> >> what is the reason for that? or how can this be prevented? >>>> >> >>>> >> I used the nbinomWaldTest and did not replace outliers with the trimmed >>>> >> mean since I have only 2 replicates per condition. >>>> >> >>>> >> I'd appreciate your help very much, >>>> >> Noa Henig >>>> >> >>>> >> [[alternative HTML version deleted]] >>>> >> >>>> >> _______________________________________________ >>>> >> Bioconductor mailing list >>>> >> Bioconductor@r-project.org >>>> >> https://stat.ethz.ch/mailman/listinfo/bioconductor >>>> >> Search the archives: >>>> >> http://news.gmane.org/gmane.science.biology.informatics.conductor >>>> >> >>>> > >>>> >>>> _______________________________________________ >>>> Bioconductor mailing list >>>> Bioconductor@r-project.org >>>> https://stat.ethz.ch/mailman/listinfo/bioconductor >>>> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor >>> >>> >> > [[alternative HTML version deleted]]