Fun with qvalue()

Hollie asked me to use the bioconductor package qvalue() to try and identify the false detection rate for the Geoduck DMR/Ls we’ve been identifying. I attempted to read the paper behind the q-value located here. It sort of, almost made sense, so I figured I’d progress to actually playing with the package.

First, we need to install the package. It comes from the Bioconductor repo. There were a lot of other things that neeed updated, so I went ahead and did that. I won’t make you suffer through the output, because it’s long and very boring.

library(readr)
library(BiocInstaller)
#biocLite("qvalue")
library(qvalue)

before we explore our Geoduck data, I want to explore their test dataset “hedenfalk” that they use in their documentation to see what an “expected” result looks like. Hedenfalk is a gene expression study on breastcancer with 3226 genes and 3170 two sample t-tests run between the genes.

stat is a vector of the t-test t-statistics, and stat0 is simulated t-tests from permutations on the labels, which generates a null distribution based on the data.

empPvals calculates p-values from supplied obsered test-statistics and simulated null distribution

data(hedenfalk)
stat <- hedenfalk$stat
stat0 <- hedenfalk$stat0
p.pooled <- empPvals(stat = stat, stat0 = stat0)
p.testspecific <- empPvals(stat = stat, stat0=stat0, pool= FALSE)

next we use a quantile-quantile plot to make sure that our observed and our simulated null distributions come from the same distribution.

qqplot(p.pooled, p.testspecific);abline(0,1)

And they do! so hooray. Lets make some q scores.

qobj <- qvalue(p.pooled)
summary(qobj)

Call:
qvalue(p = p.pooled)

pi0:    0.669926    

Cumulative number of significant calls:

          <1e-04 <0.001 <0.01 <0.025 <0.05 <0.1   <1
p-value       15     76   265    424   605  868 3170
q-value        0      0     1     73   162  319 3170
local FDR      0      0     3     30    85  167 2241

The summary above gives us cumulative sigificant results for different q-value cutoffs, as well as a theoretical number of false detectections at a given q-value cutoff.

hist(qobj)

plot(qobj)

Lets look at our data now

results_10x <- read_delim("~/Documents/Genewiz_hdd/Day10/macau/output/result.assoc.txt", 
    "\t", escape_double = FALSE, trim_ws = TRUE)

Parsed with column specification:
cols(
  id = col_character(),
  n = col_integer(),
  acpt_rate = col_double(),
  beta = col_double(),
  se_beta = col_double(),
  pvalue = col_double(),
  h = col_double(),
  se_h = col_double(),
  sigma2 = col_double(),
  se_sigma2 = col_double(),
  alpha0 = col_double(),
  se_alpha0 = col_double()
)

head(results_10x)

We care about the pvalue column, it represents the p-value that the beta parameter adds meaningful information to the overall model developed by program. Beta is the parameter for the treatment variable, so if it’s significant, that means the treatment is likely meaningful.

test.q <- qvalue(results_10x$pvalue)
summary(test.q)

Call:
qvalue(p = results_10x$pvalue)

pi0:    1   

Cumulative number of significant calls:

          <1e-04 <0.001 <0.01 <0.025 <0.05 <0.1   <1
p-value        0      1    11     22    41  103 1974
q-value        0      0     0      0     0    0 1974
local FDR      0      0     0      0     0    0    1

hist(test.q)

plot(test.q)

Those… don’t look anything like the example data.

Lets try it again, with a more liberal cutoff point applied to trimming the data.

results_5x <- read_delim("~/Documents/Genewiz_hdd/Day10/macauv2/output/5x.assoc.txt", 
    "\t", escape_double = FALSE, trim_ws = TRUE)

Parsed with column specification:
cols(
  id = col_character(),
  n = col_integer(),
  acpt_rate = col_double(),
  beta = col_double(),
  se_beta = col_double(),
  pvalue = col_double(),
  h = col_double(),
  se_h = col_double(),
  sigma2 = col_double(),
  se_sigma2 = col_double(),
  alpha0 = col_double(),
  se_alpha0 = col_double()
)

nrow(results_5x)

[1] 8936

head(results_5x)

qval.5x <- qvalue(results_5x$pvalue)
summary(qval.5x)

Call:
qvalue(p = results_5x$pvalue)

pi0:    1   

Cumulative number of significant calls:

          <1e-04 <0.001 <0.01 <0.025 <0.05 <0.1   <1
p-value        3     16    64    144   300  625 8936
q-value        0      0     1      1     2    2 8936
local FDR      0      0     1      1     1    1   18

That’s… a little more heartening. Though still it seems like a very low significance rate via q-value.

hist(qval.5x)

plot(qval.5x)

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).