Entering edit mode
Dear all, I ran a DESEQ2 script to evaluate gene ratio changes. When I have a look on the raw data it seems that the ratio is higher in the PET14 samples when in the Soil/ NoMP controls. But thr Deseq analysis shows me a positive log2fold change for otu 2 which corresponds in a lower ratio in my case. Does anyone have an explanation for this?
Thanks a lot Stephan Below are the original data
T14NoMP1 T14NoMP2 T14NoMP3 T14NoMP4 T14NoMP5 T14PE1 T14PE2 T14PE3 T14PE4 T14PE5 T14Soil1 T14Soil2 T14Soil3 T14Soil4 T14Soil5
otu001 15.49600219 0.71672122 1.04939165 7.25112383 12.41321836 4.235161e-01 20.60396068 5.080065432 15.990734851 20.145296115 0.31653702 0.55708149 0.92024388 6.41080776 23.406969736
otu002 1.71866085 2.43154776 2.02328512 1.90285255 1.76140920 2.153202e+00 2.00962395 3.776349930 2.463601552 2.107220045 0.94157150 1.21138013 1.38092935 1.13129426 1.098722348
otu003 0.04139551 0.02164607 0.01839983 0.02538947 0.02674946 9.308081e-05 0.00458436 0.006081188 0.007444051 0.002374211 0.00717658 0.01013808 0.00558299 0.00509914 0.004146075
otu004 0.33576741 0.50464516 0.40462458 0.32373456 0.34749654 8.001034e-01 0.08304846 0.064430695 0.096548372 0.123210245 0.42297635 0.25579105 0.19918321 0.32468280 0.235927196
otu005 4.77548214 3.29980457 2.59519889 7.42613315 3.54733314 3.703788e+01 52.07337724 36.632964717 18.614310414 40.830865803 39.61681424 15.55572315 31.69627507 16.42143724 18.240599756
otu006 1.14832861 1.25183746 1.14993670 1.33360559 1.07910937 8.022288e-01 3.60871444 3.704028002 2.365334311 2.730540930 0.57839931 1.33972200 1.60898705 1.00387325 1.287353742
diagdds = phyloseq_to_deseq2(psPE_T14, ~ Polymer)
calculate geometric means prior to estimate size factors
gm_mean = function(x, na.rm=TRUE){
exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
}
geoMeans = apply(counts(diagdds), 1, gm_mean)
diagdds = estimateSizeFactors(diagdds, geoMeans = geoMeans)
diagdds = DESeq(diagdds, fitType="mean")
res<-results(diagdds)
res
baseMean log2FoldChange lfcSE stat pvalue padj
<numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
otu001 7.285483 0.4001557 0.837073 0.4780415 0.632621 0.925354
otu002 1.891746 0.6738062 0.766467 0.8791071 0.379343 0.925354
otu003 0.000000 NA NA NA NA NA
otu004 0.162719 1.0391874 1.371087 0.7579296 0.448493 0.925354
otu005 19.543096 -0.0685950 0.634957 -0.1080309 0.913971 0.925354
otu006 1.474963 0.0734311 0.783753 0.0936916 0.925354 0.925354
Sorry the code and data has to bad format to be understood correctly but are you talking about normalized data ? Because you should never interpret raw counts Have you checked what is your reference when you performed differential analysis ? (this is very important)
Excuse me for the formatting, it was displayed correctly in the preview but not after submitting. Thanks for improving the format!
Regarding your comment, the reference are the NoMP_T14 and Soil_T14 samples. The data are not normalized, because the values are gene ratios. So we want to check if the gene ratios are different between treatment and control, Its just surprising, since for every other treatment the DESEQ analysis correlates with the ratio changes, but not for this one and I could not find any errors.