This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Can I use pseudocounts for differential abundance of OTUs?

Hi !!! I am working with WGS metagenome data profiled by MetaPhlAn software which gives relative abundance of taxa (no read count data of the taxa is provided). Now, I want to find out which are the significantly differentially present taxa between the test and the control samples with DESEQ package. So, my question is should I multiply the relative abundance data with some constant (e.g. 1 Million) converting the relative abundance to pseudo counts for Deseq analysis? Will it give correct result?

Thanks, dpc

deseq
• 3,177 views • link
updated by antonioggsousa antonioggsousa • written by dpc dpc

1 answer

Hi,

I don't think that you can do it, because DESeq2 uses raw counts or estimated counts: http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#input-data

This post might help: A: DESeq2 of log2(count+1) transformed TCGA htseq counts from Xenabrowser

António

• 0 views • link
written by antonioggsousa antonioggsousa

Do you have any idea how can I find out (statistically) the differentially abundant taxa from metaphlan output?

• 0 views • link
written by dpc dpc

Hi,

I never use it myself, but LEFSe seems to do what you're interested in. It uses relative abundances to find differential abundant features (OTUs/taxa/genes) ranked by effect size: https://genomebiology.biomedcentral.com/articles/10.1186/gb-2011-12-6-r60

Actually this method was developed by the same lab that developed metaphlan, so I think they are compatible, but as I said I never used them, therefore I can't be sure.

I hope this helps,

António

• 0 views • link
written by antonioggsousa antonioggsousa

Thanks, Antonio. Yes, I have already used it. I asked just to have a cross-check my data with another method if available.

• 0 views • link
written by dpc dpc

Sir, I am directly writing to you with a problem regarding correlation heatmap generation. I want a correlation heatmap showing the significance level within the heatmap cells. I want to show the P-value <0.05 as "" and P-value <0.001 as "*". For that I am following the script from here. You can find the script under the header "Stars instead of p values:" (also, I have added the codes at the end of this query). When I run with my data ("my_data") I am facing an error message

"Error in FUN(X[[i]], ...) : object 'Method' not found".

Can you please tell me how should I run the code and where am I doing a mistake? Please bear with me if it is a silly question as I am not an expert with R.

Thanks, dpc

CODE:

    library(Hmisc)
    library(car)



  mydf=df
    set.seed(12345)
    mydf[,2:5] = sapply(mydf[,2:5],jitter,amount=15)
    mydf=mydf[c(1:10,20:29,39:48,58:67),]

    # calculate r
    c = rcorr(as.matrix(mydf[sapply(mydf,is.numeric)]))$r

    # calculate p values
    p = rcorr(as.matrix(mydf[sapply(mydf,is.numeric)]))$P

    plots <- dlply(mydf, .(Method), function (x1) {
      ggplot(data.frame(subset(melt(rcorr(as.matrix(x1[sapply(x1,is.numeric)]))$r)[lower.tri(c),],Var1 != Var2),
                        pvalue=Recode(subset(melt(rcorr(as.matrix(x1[sapply(x1,is.numeric)]))$P)[lower.tri(p),],Var1 != Var2)$value , "lo:0.01 = '***'; 0.01:0.05 = '*'; else = ' ';")),
             aes(x=Var1,y=Var2,fill=value)) +
        geom_tile(aes(fill = value),colour = "white") +
        geom_text(aes(label = sprintf("%1.2f",value)), vjust = 0) + 
        geom_text(aes(label = pvalue), vjust = 1) +
        theme_bw() +
        scale_fill_gradient2(name="R^2",midpoint=0.25,low = "blue", high = "red") + 
        xlab(NULL) + 
        ylab(NULL) + 
        theme(axis.text.x=element_blank(),
              axis.text.y=element_blank(),
              axis.ticks=element_blank(),
              panel.border=element_blank()) + 
        ggtitle(x1$Method) + theme(plot.title = element_text(lineheight=1,face="bold")) + 
        geom_text(data = subset(melt(rcorr(as.matrix(x1[sapply(x1,is.numeric)]))$r),Var1==Var2),
                  aes(label=Var1),vjust=1 ) 
    })

    grid.arrange(plots$Single_ROI + theme(legend.position='none'), 
                 plots$Simple_2_ROI + theme(legend.position='none'),
                 plots$WIG_Method + theme(legend.position='none'), 
                 plots$WIG_drawn_bg + theme(legend.position='none'),
                 ncol=2, 
                 nrow=2)
• 0 views • link
written by dpc dpc

Hi @dpc,

Please, post this as another independent post/question to keep things tight. In this way others can look to your post/question and help you, because the question is not related with gene expression or differential abundance, but instead plotting. Therefore, you might get better help.

António

• 0 views • link
written by antonioggsousa antonioggsousa

Thanks Antonio. I am deleting the question here. I have posted it as a new question here . I will like to draw your attention to the query.

Thanks, dpc

• 0 views • link
written by dpc dpc

Log in to answer this question.

More posts like this