calcFDRByExperiment: FDR correction per experiment

In Carldeboer/MAUDE: Mean Alterations Using Discrete Expression

FDR correction per experiment

Description

Perform Benjamini-Hochberg FDR correction of p-values within each experiment. The returned values correspond to Q values. This is run automatically within getTilingElementwiseStats and getElementwiseStats, and doesn't generally need to be used directly.

Usage

calcFDRByExperiment(experiments, x, tails)

Arguments

experiments	a data.frame containing the headers that demarcate the screen ID, which are all also present in normNBSummaries
x	data.frame of element-level statistics, including columns for every column in 'experiments' and a column named 'p.value'
tails	whether to test for increased expression ("upper"), decreased ("lower"), or both ("both"); (defaults to "both")

Value

a numerical vector containing B-H Q values corrected separately for every experiment.

Examples

fakeReadData = data.frame(id=rep(1:10000,2), expt=c(rep("e1",10000), rep("e2",10000)), 
                          A=rpois(20000, lambda = 100), B=rpois(20000, lambda = 100),
                          C=rpois(20000, lambda = 100), D=rpois(20000, lambda = 100),
                          E=rpois(20000, lambda = 100), F=rpois(20000, lambda = 100),
                          NotSorted=rpois(20000, lambda = 100), 
                          position = rep(c(rep(NA, 1000), (1:9000)*10 + 5E7),2), 
                          chr=rep(c(rep(NA, 1000), rep("chr1", 9000)),2), 
                          negControl = rep(c(rep(TRUE,1000),rep(FALSE,9000)),2), 
                          stringsAsFactors = FALSE)
#make one region an "enhancer" and "repressor" by skewing the reads 
enhancers = data.frame(name = c("enh","repr"), start = c(40000, 70000) + 5E7, 
  end = c(40500, 70500) + 5E7, chr="chr1")
enhancerData = findOverlappingElements(fakeReadData[!is.na(fakeReadData$position),], enhancers, 
  guides.pos = "position",elements.start = "start", elements.end = "end")
readSkew=1.2 # we will scale up/down the reads in ABC and DEF by this amount
enhancerData[enhancerData$name=="enh", c("D","E","F")] = 
  floor(readSkew* enhancerData[enhancerData$name=="enh", c("D","E","F")]);
enhancerData[enhancerData$name=="repr", c("D","E","F")] = 
  floor(enhancerData[enhancerData$name=="repr", c("D","E","F")]/readSkew);
enhancerData[enhancerData$name=="repr", c("A","B","C")] = 
  floor(readSkew* enhancerData[enhancerData$name=="repr", c("A","B","C")]);
enhancerData[enhancerData$name=="enh", c("A","B","C")] = 
  floor(enhancerData[enhancerData$name=="enh", c("A","B","C")]/readSkew);
#replace the original data for these elements
fakeReadData = rbind(fakeReadData[!(fakeReadData$id %in% enhancerData$id), ],
  enhancerData[names(fakeReadData)])
#make experiments and sorting strategy
expts = unique(fakeReadData["expt"]);
curSortBins = makeBinModel(data.frame(Bin = c("A","B","C","D","E","F"), fraction = rep(0.1,6)))
curSortBins = rbind(curSortBins, curSortBins) # duplicate and use same for both expts
curSortBins$expt = c(rep(expts$expt[1],6),rep(expts$expt[2],6))
guideHits = findGuideHitsAllScreens(experiments = expts, countDataFrame=fakeReadData,
                                    binStats = curSortBins, unsortedBin = "NotSorted",
                                    negativeControl="negControl")
tilingElementStats = getTilingElementwiseStats(experiments = expts, normNBSummaries = guideHits, 
  tails = "both", chr = "chr", location="position", negativeControl = "negControl")
tilingElementStats$Q = calcFDRByExperiment(expts, tilingElementStats,"both") 
if(require("ggplot2")){
  p=ggplot(tilingElementStats, aes(x=FDR, y=Q)) +geom_point()+geom_abline(intercept=0, slope=1)+
    scale_y_log10()+scale_x_log10(); print(p)
}

Carldeboer/MAUDE documentation built on March 27, 2022, 8:50 p.m.