Question: MAST - not getting Hurdle p values
gravatar for jka8119
13 months ago by
jka811910 wrote:


I am trying to follow analysis steps as outlined in

I run into an issue when running

summaryDt <- summaryCond$datatable 

fcHurdle <- merge(summaryDt[contrast=='conditionStim' & component=='H',.(primerid, `Pr(>Chisq)`)], #hurdle P values summaryDt[contrast=='conditionStim' & component=='logFC', .(primerid, coef, ci.hi, ci.lo)], by='primerid') #logFC coefficients

My summaryCond$datatable only has 'C', 'D', 'S' and 'logFC' components but not 'H'. Because of this issue, I cannot perform subsequent MAST analysis steps. Not sure what I am doing incorrectly. I am comparing single cell RNA-seq data from two groups of cells.

Any help would be appreciated!


R version 3.4.1 (2017-06-30)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] MAST_1.4.1                              SummarizedExperiment_1.8.1              DelayedArray_0.4.1                      matrixStats_0.53.1                     
 [5] RColorBrewer_1.1-2                      rsvd_0.9                                NMF_0.20.6                              cluster_2.0.6                          
 [9] rngtools_1.2.4                          pkgmaker_0.22                           registry_0.5                            stringr_1.3.0                          
[13] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2 GenomicFeatures_1.30.3                  GenomicRanges_1.30.2                    GenomeInfoDb_1.14.0                    
[17] knitr_1.20                              data.table_1.10.4-3                     reshape2_1.4.3                          limma_3.34.8                           
[21] GSEABase_1.40.1                         graph_1.56.0                            annotate_1.56.1                         XML_3.98-1.10                          
[25] AnnotationDbi_1.40.0                    IRanges_2.12.0                          S4Vectors_0.16.0                        Biobase_2.38.0                         
[29] BiocGenerics_0.24.0                     GGally_1.3.2                            ggplot2_2.2.1                          

loaded via a namespace (and not attached):
 [1] httr_1.3.1               RMySQL_0.10.13           bit64_0.9-7              foreach_1.4.4            assertthat_0.2.0         blob_1.1.0              
 [7] GenomeInfoDbData_1.0.0   Rsamtools_1.30.0         progress_1.1.2           pillar_1.1.0             RSQLite_2.0              lattice_0.20-35         
[13] digest_0.6.15            XVector_0.18.0           colorspace_1.3-2         Matrix_1.2-12            plyr_1.8.4               pkgconfig_2.0.1         
[19] biomaRt_2.34.2           zlibbioc_1.24.0          xtable_1.8-2             scales_0.5.0             BiocParallel_1.12.0      tibble_1.4.2            
[25] lazyeval_0.2.1           magrittr_1.5             memoise_1.1.0            doParallel_1.0.11        tools_3.4.1              prettyunits_1.0.2       
[31] gridBase_0.4-7           munsell_0.4.3            Biostrings_2.46.0        compiler_3.4.1           rlang_0.2.0              grid_3.4.1              
[37] RCurl_1.95-4.10          iterators_1.0.9          bitops_1.0-6             gtable_0.2.0             codetools_0.2-15         abind_1.4-5             
[43] DBI_0.7                  reshape_0.8.7            R6_2.2.2                 GenomicAlignments_1.14.1 rtracklayer_1.38.3       bit_1.1-12              
[49] stringi_1.1.6            Rcpp_0.12.15            

mast • 306 views
ADD COMMENTlink modified 13 months ago by Andrew_McDavid190 • written 13 months ago by jka811910

Are you directly cutting and pasting code from vignette and seeing this issue?

ADD REPLYlink written 13 months ago by Andrew_McDavid190


Correct. The interesting things is that MAST correctly identifies differentially expressed genes but does not give me the Hurdle p value. 


ADD REPLYlink written 13 months ago by jka811910

I am unable to replicate.  The vignette compiles fine and summaryCond$datatable contains all the expected components:


C D H logFC S 6900 6900 2300 4600 6900

Are you trying to adapt the code in the vignette to your own data and it's not working?

ADD REPLYlink modified 13 months ago • written 13 months ago by Andrew_McDavid190

Yes. Please see below. My starting data is 'normmatrix', a numeric matrix with genes as rows and cells as columns.


>options(stringsAsFactors = F)
>genes <-, stringsAsFactors = F)
>colnames(genes) <- "gene"
>rownames(genes) <- genes$gene
>colnames(cellidentities) <- "celltype"
>cellidentities$celltype <- as.character(cellidentities$celltype)

>sle <- normmatrix[,row.names(subset(cellidentities, cellidentities$celltype %in% c("Th_conv", "Plasma_cells")))]
>cellsubset <- subset(cellidentities, cellidentities$celltype %in% c("Th_conv", "Plasma_cells"))
>scaRaw <- FromMatrix(as.matrix(sle), cellsubset, subset(genes, genes$gene %in% row.names(sle)))

>cdr2 <-colSums(assay(scaRaw)>0)

>colData(scaRaw)$cngeneson <- scale(cdr2)

>zlmCond <- zlm(~condition+cngeneson, scaRaw)

>summaryCond <- summary(zlmCond, doRLT = 'conditionPlasma_cells', logFC = TRUE, level = 0.95) 

>summaryDt <- summaryCond$datatable


ADD REPLYlink written 13 months ago by jka811910

What is the output of 

summaryCond <- summary(zlmCond, doRLT = TRUE, logFC = TRUE, level = 0.95)
table(summaryCond$datatable$component, summaryCond$datatable$contrast, exclude = NULL)


ADD REPLYlink written 13 months ago by Andrew_McDavid190

Please see below.

> table(summaryCond$datatable$component, summaryCond$datatable$contrast, exclude = NULL)
        (Intercept) cngeneson conditionPlasma_cells
  C           20922     20922                 20922
  D           20922     20922                 20922
  logFC           0     20922                 20922
  S           20922     20922                 20922

ADD REPLYlink written 13 months ago by jka811910
Answer: MAST - not getting Hurdle p values
gravatar for Andrew_McDavid
13 months ago by
Andrew_McDavid190 wrote:

There is a typo in the arguments to `summary`: `doRLT` should read `doLRT`.  No error is triggered because the generic signature includes `...`.

Perhaps the function should check that all of the dot args match and throw an error manually?

ADD COMMENTlink written 13 months ago by Andrew_McDavid190

It was the typo. Sometimes it is the smallest things... I appreciate your help.

ADD REPLYlink written 13 months ago by jka811910
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 156 users visited in the last hour