In Carldeboer/MAUDE: Mean Alterations Using Discrete Expression
knitr::opts_chunk$set(echo = TRUE)
CRISPR sort-screen evaluation test
This document includes both the code and results for the evaluation metrics used in the MAUDE paper. This should allow fairly easy extension to allow new approaches and datasets.
library("ggplot2")
library(reshape)
library(cowplot)
library(MAUDE)
library(openxlsx)
library(edgeR)
library(DESeq2)
I forgot to set a seed when doing this for the publication, and there is an element of stochasticity in the evaluation, and so the final results are likely to differ a little from the publication, especially for the datasets with less signal/noise; if you run the stochastic steps multiple times, set a different seed, or don't set the seed at all, they may also differ from this tutorial
set.seed(35263377)
CD69 screen
Load the CD69 data and process it.
#a mapping to unify bin names from Simeonov data
binmapBack = list("baseline" = "baseline", "low"="low", "medium"="medium","high"="high","back_" = "NS",
"baseline_" = "baseline", "low_"="low", "medium_"="medium", "high_"="high",
"A"="baseline", "B"="low", "E" = "medium", "F"="high")
#this comes from manually reconstructing the CD69 density curve from extended data figure 1a (Simeonov et al)
binBoundsCD69 = data.frame(Bin = c("A","F","B","E"),
fraction = c(0.65747100, 0.02792824, 0.25146688, 0.06313389),
stringsAsFactors = FALSE)
fractionalBinBounds = makeBinModel(binBoundsCD69[c("Bin","fraction")])
fractionalBinBounds = rbind(fractionalBinBounds, fractionalBinBounds)
fractionalBinBounds$screen = c(rep("1",6),rep("2",6));
#only keep bins A,B,E,F
fractionalBinBounds = fractionalBinBounds[fractionalBinBounds$Bin %in% c("A","B","E","F"),]
fractionalBinBounds$Bin = unlist(binmapBack[fractionalBinBounds$Bin]);
#load data
cd69OriginalResults = read.xlsx('https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5675716/bin/NIHMS913084-supplement-supplementary_table_1.xlsx')
cd69OriginalResults$NT = grepl("negative_control", cd69OriginalResults$gRNA_systematic_name)
cd69OriginalResults$pos = cd69OriginalResults$PAM_3primeEnd_coord;
cd69OriginalResults = unique(cd69OriginalResults)
cd69CountData = melt(cd69OriginalResults, id.vars = c("pos","NT","gRNA_systematic_name"))
cd69CountData = cd69CountData[grepl(".count$",cd69CountData$variable),]
cd69CountData$theirBin = gsub("CD69(.*)([12]).count","\\1",cd69CountData$variable)
cd69CountData$screen = gsub("CD69(.*)([12]).count","\\2",cd69CountData$variable)
cd69CountData$reads= as.numeric(cd69CountData$value); cd69CountData$value=NULL;
# convert their bin to one that is consistent
cd69CountData$Bin = unlist(binmapBack[cd69CountData$theirBin]);
binReadMat = data.frame(cast(cd69CountData[!is.na(cd69CountData$pos) | cd69CountData$NT,], pos+gRNA_systematic_name+NT+screen ~ Bin, value="reads"))
I wanted to confirm that this is how they calculated logFC in Simeonov et al. The plot below should form a straight line on y=x.
#confirm how to calc log2FC:
cd69OriginalResults$l2fc.vsbg1 = log2((1+ cd69OriginalResults$CD69high_1.count.norm)/(1+cd69OriginalResults$CD69back_1.count.norm))
cd69OriginalResults$l2fc.vsbg2 = log2((1+ cd69OriginalResults$CD69high_2.count.norm)/(1+cd69OriginalResults$CD69back_2.count.norm))
cd69OriginalResults$l2fc.hilo1 = log2((1+ cd69OriginalResults$CD69high_1.count.norm)/(1+cd69OriginalResults$CD69baseline_1.count.norm))
cd69OriginalResults$l2fc.hilo2 = log2((1+ cd69OriginalResults$CD69high_2.count.norm)/(1+cd69OriginalResults$CD69baseline_2.count.norm))
#confirm that how we just calculated log2 fc is the same as originally calculated
p=ggplot(cd69OriginalResults, aes(x=high2.l2fc,y= l2fc.vsbg2)) + geom_point(); print(p)
Run the various analysis methods on the CD69 data.
#edgeR
x <- unique(cd69OriginalResults[c("gRNA_systematic_name","CD69baseline_1.count","CD69baseline_2.count",
"CD69high_1.count", "CD69high2.count")])
row.names(x) = x$gRNA_systematic_name; x$gRNA_systematic_name=NULL;
x=as.matrix(x)
group <- factor(c(1,1,2,2))
y <- DGEList(counts=x,group=group)
y <- calcNormFactors(y)
design <- model.matrix(~group)
y <- estimateDisp(y,design)
#To perform likelihood ratio tests:
fit <- glmFit(y,design)
lrt <- glmLRT(fit,coef=2)
edgeRGuideLevel = topTags(lrt, n=nrow(x))@.Data[[1]]
edgeRGuideLevel$gRNA_systematic_name = row.names(edgeRGuideLevel);
edgeRGuideLevel$metric = edgeRGuideLevel$logFC;
edgeRGuideLevel$significance = -log(edgeRGuideLevel$PValue);
edgeRGuideLevel$method="edgeR";
#DESeq2
deseqGroups = data.frame(bin=factor(c(1,1,2,2)));
row.names(deseqGroups) = c("CD69baseline_1.count","CD69baseline_2.count",
"CD69high_1.count", "CD69high2.count");
dds <- DESeqDataSetFromMatrix(countData = x,colData = deseqGroups, design= ~ bin)
dds <- DESeq(dds)
res <- results(dds, name=resultsNames(dds)[2])
deseqGuideLevel = as.data.frame(res@listData)
deseqGuideLevel$gRNA_systematic_name =res@rownames;
deseqGuideLevel$metric = deseqGuideLevel$log2FoldChange;
deseqGuideLevel$significance = -log(deseqGuideLevel$pvalue);
deseqGuideLevel$method="DESeq2";
#MAUDE
guideLevelStatsCD69 = findGuideHitsAllScreens(unique(binReadMat["screen"]), binReadMat, fractionalBinBounds, sortBins = c("baseline","high","low","medium"), unsortedBin = "NS")
guideLevelStatsCD69$chr="chr12"
guideLevelStatsCastCD69 = cast(guideLevelStatsCD69, gRNA_systematic_name + pos+NT ~ screen, value="Z")
names(guideLevelStatsCastCD69)[ncol(guideLevelStatsCastCD69)-1:0]=c("s1","s2")
guideLevelStatsCastCD69$significance = apply(guideLevelStatsCastCD69[c("s1","s2")],1, combineZStouffer)
guideLevelStatsCastCD69$metric=apply(guideLevelStatsCastCD69[c("s1","s2")],1, mean)
guideLevelStatsCastCD69$method = "MAUDE"
#Two log fold change methods
cd69OriginalResultsHiLow = cd69OriginalResults[c("gRNA_systematic_name","l2fc.hilo1","l2fc.hilo2")]
cd69OriginalResultsVsBG = cd69OriginalResults[c("gRNA_systematic_name","l2fc.vsbg1","l2fc.vsbg2")]
cd69OriginalResultsHiLow$significance = apply(cd69OriginalResultsHiLow[2:3], 1, FUN = mean)
cd69OriginalResultsHiLow$metric = apply(cd69OriginalResultsHiLow[2:3], 1, FUN = mean)
cd69OriginalResultsVsBG$significance = apply(cd69OriginalResultsVsBG[2:3], 1, FUN = mean)
cd69OriginalResultsVsBG$metric = apply(cd69OriginalResultsVsBG[2:3], 1, FUN = mean)
cd69OriginalResultsHiLow$method="logHivsLow"
cd69OriginalResultsVsBG$method="logVsUnsorted"
Compile the results from the various methods.
# predictions
allResults = rbind(cd69OriginalResultsVsBG[c("method","gRNA_systematic_name","significance","metric")],
cd69OriginalResultsHiLow[c("method","gRNA_systematic_name","significance","metric")],
deseqGuideLevel[c("method","gRNA_systematic_name","significance","metric")],
edgeRGuideLevel[c("method","gRNA_systematic_name","significance","metric")],
guideLevelStatsCastCD69[c("method","gRNA_systematic_name","significance","metric")])
allResults = merge(allResults, cd69OriginalResults[c("gRNA_systematic_name","NT","pos")],
by="gRNA_systematic_name")
allResults = allResults[!is.na(allResults$pos) | allResults$NT,]
allResults$promoter = allResults$pos <= 9913996 & allResults$pos >= 9912997
allResults$gID = allResults$gRNA_systematic_name; allResults$gRNA_systematic_name=NULL;
allResults$locus ="CD69"
allResults$type ="CRISPRa"
allResults$celltype ="Jurkat"
TNFAIP3
Load and parse the TNFAIP3 data from Ray et al.
##read in TNFAIP3 data
binFractionsA20 = read.table(textConnection(readLines(gzcon(url("https://ftp.ncbi.nlm.nih.gov/geo/series/GSE136nnn/GSE136693/suppl/GSE136693_20190828_CRISPR_FF_bin_fractions.txt.gz")))),
sep="\t", stringsAsFactors = FALSE, header = TRUE)
CRISPRaCountsA20 = read.table(textConnection(readLines(gzcon(url("https://ftp.ncbi.nlm.nih.gov/geo/series/GSE136nnn/GSE136693/suppl/GSE136693_20190828_CRISPRa_FF_countMatrix.txt.gz")))),
sep="\t", stringsAsFactors = FALSE, header = TRUE)
CRISPRiCountsA20 = read.table(textConnection(readLines(gzcon(url("https://ftp.ncbi.nlm.nih.gov/geo/series/GSE136nnn/GSE136693/suppl/GSE136693_20190828_CRISPRi_FF_countMatrix.txt.gz")))),
sep="\t", stringsAsFactors = FALSE, header = TRUE)
crispraGuides = read.table(textConnection(readLines(gzcon(url("https://ftp.ncbi.nlm.nih.gov/geo/series/GSE136nnn/GSE136693/suppl/GSE136693_20180205_selected_CRISPRa_guides_seq.txt.gz")))),
sep="\t", stringsAsFactors = FALSE, header = TRUE)
crispraGuides$seq = gsub("(.*)(.{3})","\\1",crispraGuides$seq.w.PAM)
crispraGuides$pos = round((as.numeric(gsub("([0-9]+):([0-9]+)-([0-9]+):([+-])","\\2",
crispraGuides$guideID))+
as.numeric(gsub("([0-9]+):([0-9]+)-([0-9]+):([+-])","\\3",
crispraGuides$guideID)))/2)
CRISPRaCountsA20 = merge(CRISPRaCountsA20, unique(crispraGuides[c("seq","pos")]), by=c("seq"), all.x=TRUE)
CRISPRiCountsA20$pos = round((as.numeric(gsub("([0-9]+):([0-9]+)-([0-9]+):([+-])","\\2",
CRISPRiCountsA20$gID))+
as.numeric(gsub("([0-9]+):([0-9]+)-([0-9]+):([+-])","\\3",
CRISPRiCountsA20$gID)))/2)
binFractionsA20$expt = paste(binFractionsA20$celltype, binFractionsA20$CRISPRType,sep="_")
#merge CRISPRi and CRISPRa
A20CountData = melt(CRISPRaCountsA20, id.vars=c("seq","elementClass","element","NT","gID","pos"))
A20CountData$CRISPRType="CRISPRa";
A20CountData2 = melt(CRISPRiCountsA20, id.vars=c("seq","elementClass","element","NT","gID","pos"))
A20CountData2$CRISPRType="CRISPRi";
A20CountData = rbind(A20CountData, A20CountData2)
rm('A20CountData2');
A20CountData$count = A20CountData$value; A20CountData$value=NULL;
A20CountData$sample = A20CountData$variable; A20CountData$variable=NULL;
A20CountData$bin = gsub("(.*)_(.*)_(.*)", "\\3", A20CountData$sample);
A20CountData$screen = gsub("(.*)_(.*)_(.*)", "\\2", A20CountData$sample);
A20CountData$celltype = gsub("(.*)_(.*)_(.*)", "\\1", A20CountData$sample);
A20CountData$expt = paste(A20CountData$celltype, A20CountData$CRISPRType,sep="_")
Run each method on each experiment from Ray et al.
#combine CD69 metrics
# Pearson's r between replicates ## metricsBoth will contain all our evaluation metrics
metricsBoth =
data.frame(method=c("MAUDE","logHivsLow","logVsUnsorted"), metric="r",which="replicate_correl",
value=c(cor(guideLevelStatsCastCD69$s1[!guideLevelStatsCastCD69$NT],
guideLevelStatsCastCD69$s2[!guideLevelStatsCastCD69$NT]),
cor(cd69OriginalResults$l2fc.hilo1[!cd69OriginalResults$NT],
cd69OriginalResults$l2fc.hilo2[!cd69OriginalResults$NT]),
cor(cd69OriginalResults$l2fc.vsbg1[!cd69OriginalResults$NT],
cd69OriginalResults$l2fc.vsbg2[!cd69OriginalResults$NT])),
sig=c(cor.test(guideLevelStatsCastCD69$s1[!guideLevelStatsCastCD69$NT],
guideLevelStatsCastCD69$s2[!guideLevelStatsCastCD69$NT])$p.value,
cor.test(cd69OriginalResults$l2fc.hilo1[!cd69OriginalResults$NT],
cd69OriginalResults$l2fc.hilo2[!cd69OriginalResults$NT])$p.value,
cor.test(cd69OriginalResults$l2fc.vsbg1[!cd69OriginalResults$NT],
cd69OriginalResults$l2fc.vsbg2[!cd69OriginalResults$NT])$p.value),
locus="CD69",type="CRISPRa",celltype="Jurkat",stringsAsFactors = FALSE)
allResultsBoth = allResults; #combined results (significance, effect sizes) for CD69 and A20 screens
for (e in unique(binFractionsA20$expt)){
curCelltype = gsub("(.*)_(.*)", "\\1", e);
curtype = gsub("(.*)_(.*)", "\\2", e);
curA20CountData = unique(A20CountData[A20CountData$expt==e,
c("seq","pos", "NT","gID","count","screen","bin")])
curA20CountDataTotals = cast(curA20CountData, screen +bin~ ., fun.aggregate = sum, value="count")
names(curA20CountDataTotals)[3] = "total";
curA20CountData = merge(curA20CountData, curA20CountDataTotals, by=c("screen","bin"))
curA20CountData$CPM = curA20CountData$count/curA20CountData$total * 1E6;
curCPMMat = cast(curA20CountData, seq + NT + gID + screen + pos ~ bin, value="CPM")
curCPMMat$l2fc_hilo = log2((1+curCPMMat$F)/(1+curCPMMat$A))
if(curtype=="CRISRPi"){
curCPMMat$l2fc_vsbg = log2((1+curCPMMat$NS)/(1+curCPMMat$A))
}else{ #CRISPRa
curCPMMat$l2fc_vsbg = log2((1+curCPMMat$F)/(1+curCPMMat$NS))
}
curBins = as.data.frame(melt(binFractionsA20[binFractionsA20$expt==e,],
id.vars = c("celltype","screen","CRISPRType","expt")))
names(curBins)[ncol(curBins) - (1:0)] = c("Bin","fraction")
curBins2 = data.frame();
for (s in unique(curBins$screen)){
curBins3 = makeBinModel(curBins[curBins$screen==s,c("Bin","fraction")])
curBins3$screen = s;
curBins2 = rbind(curBins2, curBins3)
}
curBins2$Bin = as.character(curBins2$Bin);
curCountMat = cast(curA20CountData, seq + NT + gID +pos + screen ~ bin, value="count")
guideLevelStats = findGuideHitsAllScreens(experiments = unique(curCountMat["screen"]),
countDataFrame = curCountMat, binStats = curBins2,
sortBins = c("A","B","C","D","E","F"), unsortedBin = "NS",
negativeControl="NT")
guideLevelStatsCast = cast(guideLevelStats, gID + pos+NT ~ screen, value="Z")
#names(guideLevelStatsCast)[4:ncol(guideLevelStatsCast)]=sprintf("s%i", 1:(ncol(guideLevelStatsCast)-3))
maudeZs = guideLevelStatsCast;
guideLevelStatsCast$significance = apply(maudeZs[unique(curA20CountData$screen)],1, combineZStouffer)
guideLevelStatsCast$metric=apply(maudeZs[unique(curA20CountData$screen)],1, mean)
guideLevelStatsCast$method = "MAUDE"
### EdgeR
library(edgeR)
x= cast(unique(curA20CountData[curA20CountData$bin %in% c("A","F"), c("bin","gID","screen","count")]), gID ~ screen + bin, value="count")
row.names(x) = x$gID; x$gID=NULL;
x=as.matrix(x)
group = grepl("_F",colnames(x))+1
group <- factor(group)
y <- DGEList(counts=x,group=group)
y <- calcNormFactors(y)
design <- model.matrix(~group)
y <- estimateDisp(y,design)
#To perform likelihood ratio tests:
fit <- glmFit(y,design)
lrt <- glmLRT(fit,coef=2)
edgeRGuideLevel = topTags(lrt, n=nrow(x))@.Data[[1]]
edgeRGuideLevel$gID = row.names(edgeRGuideLevel);
edgeRGuideLevel$metric = edgeRGuideLevel$logFC;
edgeRGuideLevel$significance = -log(edgeRGuideLevel$PValue);
edgeRGuideLevel$method="edgeR";
### DEseq
library(DESeq2)
deseqGroups = data.frame(bin=group);
row.names(deseqGroups) = colnames(x);
dds <- DESeqDataSetFromMatrix(countData = x,colData = deseqGroups, design= ~ bin)
dds <- DESeq(dds)
#resultsNames(dds) # lists the coefficients
res <- results(dds, name=resultsNames(dds)[2])
stopifnot(resultsNames(dds)[1]=="Intercept")
deseqGuideLevel = as.data.frame(res@listData)
deseqGuideLevel$gID =res@rownames;
deseqGuideLevel$metric = deseqGuideLevel$log2FoldChange;
deseqGuideLevel$significance = -log(deseqGuideLevel$pvalue);
deseqGuideLevel$method="DESeq2";
curLRHiLow = cast(unique(curCPMMat[c("gID","NT","screen","l2fc_hilo")]), gID + NT ~ screen, value="l2fc_hilo")
curLRVsBG = cast(unique(curCPMMat[c("gID","NT","screen","l2fc_vsbg")]), gID + NT ~ screen, value="l2fc_vsbg")
numSamples = ncol(curLRHiLow)-2;
sampleNames = unique(curCPMMat$screen)
curLRVsBG$significance = apply(curLRVsBG[3:(numSamples+2)], MARGIN = 1, FUN = mean);
curLRVsBG$metric = apply(curLRVsBG[3:(numSamples+2)], MARGIN = 1, FUN = mean);
curLRVsBG$method="logVsUnsorted"
curLRHiLow$significance = apply(curLRHiLow[3:(numSamples+2)], MARGIN = 1, FUN = mean);
curLRHiLow$metric = apply(curLRHiLow[3:(numSamples+2)], MARGIN = 1, FUN = mean);
curLRHiLow$method="logHivsLow"
#compile results for A20
curResults = rbind(unique(curLRHiLow[c("method","gID","significance","metric")]),
unique(curLRVsBG[c("method","gID","significance","metric")]),
deseqGuideLevel[c("method","gID","significance","metric")],
edgeRGuideLevel[c("method","gID","significance","metric")],
unique(guideLevelStatsCast[c("method","gID","significance","metric")]))
curResults = merge(curResults, unique(curCPMMat[c("gID","NT","pos")]), by="gID")
curResults = curResults[!is.na(curResults$pos) | curResults$NT,]
curResults$promoter = grepl("TNFAIP3", curResults$gID) |
(curResults$pos <= 138189439 & curResults$pos >= 138187040) # chr6:138188077-138188379;138187040
#append the current results to all
curResults$locus ="TNFAIP3"
curResults$type =curtype
curResults$celltype =curCelltype
allResultsBoth = rbind(allResultsBoth, curResults);
# (1) similarity between the effect sizes estimated per replicate,
corLRHiLow = cor(curLRHiLow[!curLRHiLow$NT, 3:(3+numSamples-1)])
corLRVsBG = cor(curLRVsBG[!curLRVsBG$NT, 3:(3+numSamples-1)])
maudeZCors = cor(maudeZs[!maudeZs$NT, 4:ncol(maudeZs)])
maudeCorP=1
maudeCorR=-1
corLRHiLowP=1
corLRHiLowR=-1
corLRVsBGP=1;
corLRVsBGR=-1
#select the best inter-replicate correlation for each of the three approaches for which this is possible
for(i in 1:(length(sampleNames)-1)){
for(j in (i+1):length(sampleNames)){
curR = cor(maudeZs[!maudeZs$NT, sampleNames[i]], maudeZs[!maudeZs$NT, sampleNames[j]])
curP = cor.test(maudeZs[!maudeZs$NT, sampleNames[i]], maudeZs[!maudeZs$NT, sampleNames[j]])$p.value
if (maudeCorR < curR){
maudeCorR = curR;
maudeCorP = curP;
}
curR = cor(curLRVsBG[!curLRVsBG$NT, sampleNames[i]], curLRVsBG[!curLRVsBG$NT, sampleNames[j]])
curP = cor.test(curLRVsBG[!curLRVsBG$NT, sampleNames[i]],
curLRVsBG[!curLRVsBG$NT, sampleNames[j]])$p.value
if (corLRVsBGR < curR){
corLRVsBGR = curR;
corLRVsBGP = curP;
}
curR = cor(curLRHiLow[!curLRHiLow$NT, sampleNames[i]], curLRHiLow[!curLRHiLow$NT, sampleNames[j]])
curP = cor.test(curLRHiLow[!curLRHiLow$NT, sampleNames[i]],
curLRHiLow[!curLRHiLow$NT, sampleNames[j]])$p.value
if (corLRHiLowR < curR){
corLRHiLowR = curR;
corLRHiLowP = curP;
}
}
}
metricsBoth = rbind(metricsBoth,
data.frame(method=c("MAUDE","logHivsLow","logVsUnsorted"),
metric="r",which="replicate_correl",
value=c(maudeCorR, corLRHiLowR, corLRVsBGR),
sig=c(maudeCorP, corLRHiLowP, corLRVsBGP),
locus ="TNFAIP3", type =curtype, celltype =curCelltype,
stringsAsFactors = FALSE))
}
Here are some functions that I will make use of below.
#Some useful functions
ranksumROC = function(x,y,na.rm=TRUE,...){
if (na.rm){
x=na.rm(x);
y=na.rm(y);
}
curTest = wilcox.test(x,y,...);
curTest$AUROC = (curTest$statistic/length(x))/length(y)
return(curTest)
}
na.rm = function(x){ x[!is.na(x)]}
Evaluation
The inter-replicate correlations were calculated in the sections above and are stored in metricsBoth. Below, the two remaining metrics are calculated (distinguishing promoter vs other targeting guides, and adjacent vs randomly paired guides).
# (1) similarity between the effect sizes estimated per replicate,
# (above)
metricsBoth$significant= metricsBoth$sig < 0.01;
# Other evaluation metrics
allExpts = unique(allResultsBoth[c("celltype","locus","type")])
for (ei in 1:nrow(allExpts)){
curCelltype = allExpts$celltype[ei]
curtype = allExpts$type[ei];
curLocus = allExpts$locus[ei];
curResults = allResultsBoth[allResultsBoth$celltype==curCelltype &
allResultsBoth$type==curtype & allResultsBoth$locus==curLocus,]
for(m in unique(curResults$method)){
curData = curResults[curResults$method==m & !curResults$NT,]
# (2) similarity in effect size between adjacent guides
curData = curData[order(curData$pos),]
guideEffectDistances =
data.frame(method = m, random=FALSE,
difference = abs(curData$metric[2:nrow(curData)] - curData$metric[1:(nrow(curData)-1)]),
dist =abs(curData$pos[2:nrow(curData)] - curData$pos[1:(nrow(curData)-1)]),
stringsAsFactors = FALSE)
guideEffectDistances = guideEffectDistances[guideEffectDistances$dist < 100,]
guideEffectDistances$dist=NULL;
### changed 10 in next line to 100 to make this more robust
curData = curData[sample(nrow(curData), size = nrow(curData)*100, replace = TRUE),]
guideEffectDistances =
rbind(guideEffectDistances,
data.frame(method = m, random=TRUE,
difference = abs(curData$metric[2:nrow(curData)] -
curData$metric[1:(nrow(curData)-1)]), stringsAsFactors = FALSE));
# random should have more different than adjacent
curRS = ranksumROC(guideEffectDistances$difference[guideEffectDistances$method==m &
guideEffectDistances$random],
guideEffectDistances$difference[guideEffectDistances$method==m &
!guideEffectDistances$random])
metricsBoth = rbind(metricsBoth, data.frame(method=m, metric="AUROC-0.5",
which="adjacent_vs_random",value = curRS$AUROC-0.5,
locus=curLocus,celltype=curCelltype, type=curtype,
sig=curRS$p.value, significant = curRS$p.value < 0.01))
# (3) ability to distinguish promoter-targeting guides from other guides.
if (curtype=="CRISPRi" & !(m %in% c("edgeR","DESeq2"))){
# edgeR and DESeq2 are reversed for CRISPRi
# non promoter should have higher effect than promoter (more -ve)
curRS = ranksumROC(curResults$significance[curResults$method==m & !curResults$NT &
!curResults$promoter],
curResults$significance[curResults$method==m & !curResults$NT &
curResults$promoter])
}else{
# promoter should have larger effect than non-promoter
curRS = ranksumROC(curResults$significance[curResults$method==m & !curResults$NT &
curResults$promoter],
curResults$significance[curResults$method==m & !curResults$NT &
!curResults$promoter])
}
metricsBoth = rbind(metricsBoth, data.frame(method=m, metric="AUROC-0.5", which="promoter_vs_T",
value = curRS$AUROC-0.5, locus=curLocus,
celltype=curCelltype, type=curtype, sig=curRS$p.value,
significant = curRS$p.value < 0.01))
}
}
##compile all metrics; label the best in each test and whether any tests were significant (P<0.01)
metricsBoth2 = metricsBoth;
metricsBoth2$method = factor(as.character(metricsBoth2$method),
levels = c("logVsUnsorted","logHivsLow","DESeq2","edgeR","MAUDE"))
metricsBoth2Best = cast(metricsBoth2, which + locus + type+celltype ~ ., value="value",
fun.aggregate = max)
names(metricsBoth2Best)[ncol(metricsBoth2Best)] = "best"
metricsBoth2 = merge(metricsBoth2, metricsBoth2Best, by = c("which","locus","type","celltype"))
metricsBoth2AnySig = cast(metricsBoth2[colnames(metricsBoth2)!="value"],
which + locus + type+celltype ~ ., value="significant", fun.aggregate = any)
names(metricsBoth2AnySig)[ncol(metricsBoth2AnySig)] = "anySig"
metricsBoth2 = merge(metricsBoth2, metricsBoth2AnySig, by = c("which","locus","type","celltype"))
metricsBoth2$isBest = metricsBoth2$value==metricsBoth2$best;
metricsBoth2$isBestNA = metricsBoth2$isBest;
metricsBoth2$isBestNA[!metricsBoth2$isBestNA]=NA;
metricsBoth2$pctOfMax = metricsBoth2$value/metricsBoth2$best * 100;
#fill in NAs for edgeR and DESeq2 which cannot have inter-replicate correlations
temp = metricsBoth2[metricsBoth2$metric=="r",]
temp = temp[1:2,];
temp$method = c("edgeR","DESeq2")
temp$value=NA; temp$isBest=NA; temp$significant=FALSE; temp$pctOfMax=NA; temp$isBestNA=NA;
metricsBoth2 = rbind(metricsBoth2, temp)
Finally make the graph with all evaluation metrics as shown in the MAUDE paper.
#make the final evaluation graph
p1 = ggplot(metricsBoth2[metricsBoth2$which=="adjacent_vs_random",],
aes(x=method, fill=value, y=paste(locus,type,celltype))) + geom_tile() +
geom_text(data=metricsBoth2[metricsBoth2$which=="adjacent_vs_random" & metricsBoth2$anySig,],
aes(label="*",colour=isBestNA),show.legend = FALSE) +theme_bw() +
scale_fill_gradient2(high="red", low="blue", mid="black") +
theme(legend.position="top", axis.text.x = element_text(hjust=1, angle=45),
axis.title.y = element_blank())+scale_colour_manual(values = c("green"), na.value=NA) +
scale_y_discrete(expand=c(0,0))+scale_x_discrete(expand=c(0,0))+ggtitle("Adjacent vs\nrandom guides");
p2 = ggplot(metricsBoth2[metricsBoth2$which=="promoter_vs_T",],
aes(x=method, fill=value, y=paste(locus,type,celltype))) + geom_tile() +
geom_text(data=metricsBoth2[metricsBoth2$which=="promoter_vs_T" & metricsBoth2$anySig,],
aes(label="*",colour=isBestNA),show.legend = FALSE) +theme_bw() +
scale_fill_gradient2(high="red", low="blue", mid="black") +
theme(legend.position="top", axis.text.x = element_text(hjust=1, angle=45),
axis.title.y = element_blank())+scale_colour_manual(values = c("green"), na.value=NA)+
scale_y_discrete(expand=c(0,0))+scale_x_discrete(expand=c(0,0))+
ggtitle("Promoter vs other\ntargeting guides");
p3 = ggplot(metricsBoth2[metricsBoth2$which=="replicate_correl",],
aes(x=method, fill=value, y=paste(locus,type,celltype))) + geom_tile() +
geom_text(data=metricsBoth2[metricsBoth2$which=="replicate_correl" & metricsBoth2$anySig,],
aes(label="*",colour=isBestNA),show.legend = FALSE) +theme_bw() +
scale_fill_gradient2(high="red", low="blue", mid="black") +
theme(legend.position="top", axis.text.x = element_text(hjust=1, angle=45),
axis.title.y = element_blank())+scale_colour_manual(values = c("green"), na.value=NA)+
scale_y_discrete(expand=c(0,0))+scale_x_discrete(expand=c(0,0))+ggtitle("Replicate\ncorrelations");
g= plot_grid(p1,p2,p3, align = 'h', nrow = 1); print(g)
Session info
sessionInfo()
Carldeboer/MAUDE documentation built on March 27, 2022, 8:50 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = TRUE)
CRISPR sort-screen evaluation test
This document includes both the code and results for the evaluation metrics used in the MAUDE paper. This should allow fairly easy extension to allow new approaches and datasets.
library("ggplot2") library(reshape) library(cowplot) library(MAUDE) library(openxlsx) library(edgeR) library(DESeq2)
I forgot to set a seed when doing this for the publication, and there is an element of stochasticity in the evaluation, and so the final results are likely to differ a little from the publication, especially for the datasets with less signal/noise; if you run the stochastic steps multiple times, set a different seed, or don't set the seed at all, they may also differ from this tutorial
set.seed(35263377)
CD69 screen
Load the CD69 data and process it.
#a mapping to unify bin names from Simeonov data binmapBack = list("baseline" = "baseline", "low"="low", "medium"="medium","high"="high","back_" = "NS", "baseline_" = "baseline", "low_"="low", "medium_"="medium", "high_"="high", "A"="baseline", "B"="low", "E" = "medium", "F"="high") #this comes from manually reconstructing the CD69 density curve from extended data figure 1a (Simeonov et al) binBoundsCD69 = data.frame(Bin = c("A","F","B","E"), fraction = c(0.65747100, 0.02792824, 0.25146688, 0.06313389), stringsAsFactors = FALSE) fractionalBinBounds = makeBinModel(binBoundsCD69[c("Bin","fraction")]) fractionalBinBounds = rbind(fractionalBinBounds, fractionalBinBounds) fractionalBinBounds$screen = c(rep("1",6),rep("2",6)); #only keep bins A,B,E,F fractionalBinBounds = fractionalBinBounds[fractionalBinBounds$Bin %in% c("A","B","E","F"),] fractionalBinBounds$Bin = unlist(binmapBack[fractionalBinBounds$Bin]); #load data cd69OriginalResults = read.xlsx('https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5675716/bin/NIHMS913084-supplement-supplementary_table_1.xlsx') cd69OriginalResults$NT = grepl("negative_control", cd69OriginalResults$gRNA_systematic_name) cd69OriginalResults$pos = cd69OriginalResults$PAM_3primeEnd_coord; cd69OriginalResults = unique(cd69OriginalResults) cd69CountData = melt(cd69OriginalResults, id.vars = c("pos","NT","gRNA_systematic_name")) cd69CountData = cd69CountData[grepl(".count$",cd69CountData$variable),] cd69CountData$theirBin = gsub("CD69(.*)([12]).count","\\1",cd69CountData$variable) cd69CountData$screen = gsub("CD69(.*)([12]).count","\\2",cd69CountData$variable) cd69CountData$reads= as.numeric(cd69CountData$value); cd69CountData$value=NULL; # convert their bin to one that is consistent cd69CountData$Bin = unlist(binmapBack[cd69CountData$theirBin]); binReadMat = data.frame(cast(cd69CountData[!is.na(cd69CountData$pos) | cd69CountData$NT,], pos+gRNA_systematic_name+NT+screen ~ Bin, value="reads"))
I wanted to confirm that this is how they calculated logFC in Simeonov et al. The plot below should form a straight line on y=x.
#confirm how to calc log2FC: cd69OriginalResults$l2fc.vsbg1 = log2((1+ cd69OriginalResults$CD69high_1.count.norm)/(1+cd69OriginalResults$CD69back_1.count.norm)) cd69OriginalResults$l2fc.vsbg2 = log2((1+ cd69OriginalResults$CD69high_2.count.norm)/(1+cd69OriginalResults$CD69back_2.count.norm)) cd69OriginalResults$l2fc.hilo1 = log2((1+ cd69OriginalResults$CD69high_1.count.norm)/(1+cd69OriginalResults$CD69baseline_1.count.norm)) cd69OriginalResults$l2fc.hilo2 = log2((1+ cd69OriginalResults$CD69high_2.count.norm)/(1+cd69OriginalResults$CD69baseline_2.count.norm)) #confirm that how we just calculated log2 fc is the same as originally calculated p=ggplot(cd69OriginalResults, aes(x=high2.l2fc,y= l2fc.vsbg2)) + geom_point(); print(p)
Run the various analysis methods on the CD69 data.
#edgeR x <- unique(cd69OriginalResults[c("gRNA_systematic_name","CD69baseline_1.count","CD69baseline_2.count", "CD69high_1.count", "CD69high2.count")]) row.names(x) = x$gRNA_systematic_name; x$gRNA_systematic_name=NULL; x=as.matrix(x) group <- factor(c(1,1,2,2)) y <- DGEList(counts=x,group=group) y <- calcNormFactors(y) design <- model.matrix(~group) y <- estimateDisp(y,design) #To perform likelihood ratio tests: fit <- glmFit(y,design) lrt <- glmLRT(fit,coef=2) edgeRGuideLevel = topTags(lrt, n=nrow(x))@.Data[[1]] edgeRGuideLevel$gRNA_systematic_name = row.names(edgeRGuideLevel); edgeRGuideLevel$metric = edgeRGuideLevel$logFC; edgeRGuideLevel$significance = -log(edgeRGuideLevel$PValue); edgeRGuideLevel$method="edgeR"; #DESeq2 deseqGroups = data.frame(bin=factor(c(1,1,2,2))); row.names(deseqGroups) = c("CD69baseline_1.count","CD69baseline_2.count", "CD69high_1.count", "CD69high2.count"); dds <- DESeqDataSetFromMatrix(countData = x,colData = deseqGroups, design= ~ bin) dds <- DESeq(dds) res <- results(dds, name=resultsNames(dds)[2]) deseqGuideLevel = as.data.frame(res@listData) deseqGuideLevel$gRNA_systematic_name =res@rownames; deseqGuideLevel$metric = deseqGuideLevel$log2FoldChange; deseqGuideLevel$significance = -log(deseqGuideLevel$pvalue); deseqGuideLevel$method="DESeq2"; #MAUDE guideLevelStatsCD69 = findGuideHitsAllScreens(unique(binReadMat["screen"]), binReadMat, fractionalBinBounds, sortBins = c("baseline","high","low","medium"), unsortedBin = "NS") guideLevelStatsCD69$chr="chr12" guideLevelStatsCastCD69 = cast(guideLevelStatsCD69, gRNA_systematic_name + pos+NT ~ screen, value="Z") names(guideLevelStatsCastCD69)[ncol(guideLevelStatsCastCD69)-1:0]=c("s1","s2") guideLevelStatsCastCD69$significance = apply(guideLevelStatsCastCD69[c("s1","s2")],1, combineZStouffer) guideLevelStatsCastCD69$metric=apply(guideLevelStatsCastCD69[c("s1","s2")],1, mean) guideLevelStatsCastCD69$method = "MAUDE" #Two log fold change methods cd69OriginalResultsHiLow = cd69OriginalResults[c("gRNA_systematic_name","l2fc.hilo1","l2fc.hilo2")] cd69OriginalResultsVsBG = cd69OriginalResults[c("gRNA_systematic_name","l2fc.vsbg1","l2fc.vsbg2")] cd69OriginalResultsHiLow$significance = apply(cd69OriginalResultsHiLow[2:3], 1, FUN = mean) cd69OriginalResultsHiLow$metric = apply(cd69OriginalResultsHiLow[2:3], 1, FUN = mean) cd69OriginalResultsVsBG$significance = apply(cd69OriginalResultsVsBG[2:3], 1, FUN = mean) cd69OriginalResultsVsBG$metric = apply(cd69OriginalResultsVsBG[2:3], 1, FUN = mean) cd69OriginalResultsHiLow$method="logHivsLow" cd69OriginalResultsVsBG$method="logVsUnsorted"
Compile the results from the various methods.
# predictions allResults = rbind(cd69OriginalResultsVsBG[c("method","gRNA_systematic_name","significance","metric")], cd69OriginalResultsHiLow[c("method","gRNA_systematic_name","significance","metric")], deseqGuideLevel[c("method","gRNA_systematic_name","significance","metric")], edgeRGuideLevel[c("method","gRNA_systematic_name","significance","metric")], guideLevelStatsCastCD69[c("method","gRNA_systematic_name","significance","metric")]) allResults = merge(allResults, cd69OriginalResults[c("gRNA_systematic_name","NT","pos")], by="gRNA_systematic_name") allResults = allResults[!is.na(allResults$pos) | allResults$NT,] allResults$promoter = allResults$pos <= 9913996 & allResults$pos >= 9912997 allResults$gID = allResults$gRNA_systematic_name; allResults$gRNA_systematic_name=NULL; allResults$locus ="CD69" allResults$type ="CRISPRa" allResults$celltype ="Jurkat"
TNFAIP3
Load and parse the TNFAIP3 data from Ray et al.
##read in TNFAIP3 data binFractionsA20 = read.table(textConnection(readLines(gzcon(url("https://ftp.ncbi.nlm.nih.gov/geo/series/GSE136nnn/GSE136693/suppl/GSE136693_20190828_CRISPR_FF_bin_fractions.txt.gz")))), sep="\t", stringsAsFactors = FALSE, header = TRUE) CRISPRaCountsA20 = read.table(textConnection(readLines(gzcon(url("https://ftp.ncbi.nlm.nih.gov/geo/series/GSE136nnn/GSE136693/suppl/GSE136693_20190828_CRISPRa_FF_countMatrix.txt.gz")))), sep="\t", stringsAsFactors = FALSE, header = TRUE) CRISPRiCountsA20 = read.table(textConnection(readLines(gzcon(url("https://ftp.ncbi.nlm.nih.gov/geo/series/GSE136nnn/GSE136693/suppl/GSE136693_20190828_CRISPRi_FF_countMatrix.txt.gz")))), sep="\t", stringsAsFactors = FALSE, header = TRUE) crispraGuides = read.table(textConnection(readLines(gzcon(url("https://ftp.ncbi.nlm.nih.gov/geo/series/GSE136nnn/GSE136693/suppl/GSE136693_20180205_selected_CRISPRa_guides_seq.txt.gz")))), sep="\t", stringsAsFactors = FALSE, header = TRUE) crispraGuides$seq = gsub("(.*)(.{3})","\\1",crispraGuides$seq.w.PAM) crispraGuides$pos = round((as.numeric(gsub("([0-9]+):([0-9]+)-([0-9]+):([+-])","\\2", crispraGuides$guideID))+ as.numeric(gsub("([0-9]+):([0-9]+)-([0-9]+):([+-])","\\3", crispraGuides$guideID)))/2) CRISPRaCountsA20 = merge(CRISPRaCountsA20, unique(crispraGuides[c("seq","pos")]), by=c("seq"), all.x=TRUE) CRISPRiCountsA20$pos = round((as.numeric(gsub("([0-9]+):([0-9]+)-([0-9]+):([+-])","\\2", CRISPRiCountsA20$gID))+ as.numeric(gsub("([0-9]+):([0-9]+)-([0-9]+):([+-])","\\3", CRISPRiCountsA20$gID)))/2) binFractionsA20$expt = paste(binFractionsA20$celltype, binFractionsA20$CRISPRType,sep="_") #merge CRISPRi and CRISPRa A20CountData = melt(CRISPRaCountsA20, id.vars=c("seq","elementClass","element","NT","gID","pos")) A20CountData$CRISPRType="CRISPRa"; A20CountData2 = melt(CRISPRiCountsA20, id.vars=c("seq","elementClass","element","NT","gID","pos")) A20CountData2$CRISPRType="CRISPRi"; A20CountData = rbind(A20CountData, A20CountData2) rm('A20CountData2'); A20CountData$count = A20CountData$value; A20CountData$value=NULL; A20CountData$sample = A20CountData$variable; A20CountData$variable=NULL; A20CountData$bin = gsub("(.*)_(.*)_(.*)", "\\3", A20CountData$sample); A20CountData$screen = gsub("(.*)_(.*)_(.*)", "\\2", A20CountData$sample); A20CountData$celltype = gsub("(.*)_(.*)_(.*)", "\\1", A20CountData$sample); A20CountData$expt = paste(A20CountData$celltype, A20CountData$CRISPRType,sep="_")
Run each method on each experiment from Ray et al.
#combine CD69 metrics # Pearson's r between replicates ## metricsBoth will contain all our evaluation metrics metricsBoth = data.frame(method=c("MAUDE","logHivsLow","logVsUnsorted"), metric="r",which="replicate_correl", value=c(cor(guideLevelStatsCastCD69$s1[!guideLevelStatsCastCD69$NT], guideLevelStatsCastCD69$s2[!guideLevelStatsCastCD69$NT]), cor(cd69OriginalResults$l2fc.hilo1[!cd69OriginalResults$NT], cd69OriginalResults$l2fc.hilo2[!cd69OriginalResults$NT]), cor(cd69OriginalResults$l2fc.vsbg1[!cd69OriginalResults$NT], cd69OriginalResults$l2fc.vsbg2[!cd69OriginalResults$NT])), sig=c(cor.test(guideLevelStatsCastCD69$s1[!guideLevelStatsCastCD69$NT], guideLevelStatsCastCD69$s2[!guideLevelStatsCastCD69$NT])$p.value, cor.test(cd69OriginalResults$l2fc.hilo1[!cd69OriginalResults$NT], cd69OriginalResults$l2fc.hilo2[!cd69OriginalResults$NT])$p.value, cor.test(cd69OriginalResults$l2fc.vsbg1[!cd69OriginalResults$NT], cd69OriginalResults$l2fc.vsbg2[!cd69OriginalResults$NT])$p.value), locus="CD69",type="CRISPRa",celltype="Jurkat",stringsAsFactors = FALSE) allResultsBoth = allResults; #combined results (significance, effect sizes) for CD69 and A20 screens for (e in unique(binFractionsA20$expt)){ curCelltype = gsub("(.*)_(.*)", "\\1", e); curtype = gsub("(.*)_(.*)", "\\2", e); curA20CountData = unique(A20CountData[A20CountData$expt==e, c("seq","pos", "NT","gID","count","screen","bin")]) curA20CountDataTotals = cast(curA20CountData, screen +bin~ ., fun.aggregate = sum, value="count") names(curA20CountDataTotals)[3] = "total"; curA20CountData = merge(curA20CountData, curA20CountDataTotals, by=c("screen","bin")) curA20CountData$CPM = curA20CountData$count/curA20CountData$total * 1E6; curCPMMat = cast(curA20CountData, seq + NT + gID + screen + pos ~ bin, value="CPM") curCPMMat$l2fc_hilo = log2((1+curCPMMat$F)/(1+curCPMMat$A)) if(curtype=="CRISRPi"){ curCPMMat$l2fc_vsbg = log2((1+curCPMMat$NS)/(1+curCPMMat$A)) }else{ #CRISPRa curCPMMat$l2fc_vsbg = log2((1+curCPMMat$F)/(1+curCPMMat$NS)) } curBins = as.data.frame(melt(binFractionsA20[binFractionsA20$expt==e,], id.vars = c("celltype","screen","CRISPRType","expt"))) names(curBins)[ncol(curBins) - (1:0)] = c("Bin","fraction") curBins2 = data.frame(); for (s in unique(curBins$screen)){ curBins3 = makeBinModel(curBins[curBins$screen==s,c("Bin","fraction")]) curBins3$screen = s; curBins2 = rbind(curBins2, curBins3) } curBins2$Bin = as.character(curBins2$Bin); curCountMat = cast(curA20CountData, seq + NT + gID +pos + screen ~ bin, value="count") guideLevelStats = findGuideHitsAllScreens(experiments = unique(curCountMat["screen"]), countDataFrame = curCountMat, binStats = curBins2, sortBins = c("A","B","C","D","E","F"), unsortedBin = "NS", negativeControl="NT") guideLevelStatsCast = cast(guideLevelStats, gID + pos+NT ~ screen, value="Z") #names(guideLevelStatsCast)[4:ncol(guideLevelStatsCast)]=sprintf("s%i", 1:(ncol(guideLevelStatsCast)-3)) maudeZs = guideLevelStatsCast; guideLevelStatsCast$significance = apply(maudeZs[unique(curA20CountData$screen)],1, combineZStouffer) guideLevelStatsCast$metric=apply(maudeZs[unique(curA20CountData$screen)],1, mean) guideLevelStatsCast$method = "MAUDE" ### EdgeR library(edgeR) x= cast(unique(curA20CountData[curA20CountData$bin %in% c("A","F"), c("bin","gID","screen","count")]), gID ~ screen + bin, value="count") row.names(x) = x$gID; x$gID=NULL; x=as.matrix(x) group = grepl("_F",colnames(x))+1 group <- factor(group) y <- DGEList(counts=x,group=group) y <- calcNormFactors(y) design <- model.matrix(~group) y <- estimateDisp(y,design) #To perform likelihood ratio tests: fit <- glmFit(y,design) lrt <- glmLRT(fit,coef=2) edgeRGuideLevel = topTags(lrt, n=nrow(x))@.Data[[1]] edgeRGuideLevel$gID = row.names(edgeRGuideLevel); edgeRGuideLevel$metric = edgeRGuideLevel$logFC; edgeRGuideLevel$significance = -log(edgeRGuideLevel$PValue); edgeRGuideLevel$method="edgeR"; ### DEseq library(DESeq2) deseqGroups = data.frame(bin=group); row.names(deseqGroups) = colnames(x); dds <- DESeqDataSetFromMatrix(countData = x,colData = deseqGroups, design= ~ bin) dds <- DESeq(dds) #resultsNames(dds) # lists the coefficients res <- results(dds, name=resultsNames(dds)[2]) stopifnot(resultsNames(dds)[1]=="Intercept") deseqGuideLevel = as.data.frame(res@listData) deseqGuideLevel$gID =res@rownames; deseqGuideLevel$metric = deseqGuideLevel$log2FoldChange; deseqGuideLevel$significance = -log(deseqGuideLevel$pvalue); deseqGuideLevel$method="DESeq2"; curLRHiLow = cast(unique(curCPMMat[c("gID","NT","screen","l2fc_hilo")]), gID + NT ~ screen, value="l2fc_hilo") curLRVsBG = cast(unique(curCPMMat[c("gID","NT","screen","l2fc_vsbg")]), gID + NT ~ screen, value="l2fc_vsbg") numSamples = ncol(curLRHiLow)-2; sampleNames = unique(curCPMMat$screen) curLRVsBG$significance = apply(curLRVsBG[3:(numSamples+2)], MARGIN = 1, FUN = mean); curLRVsBG$metric = apply(curLRVsBG[3:(numSamples+2)], MARGIN = 1, FUN = mean); curLRVsBG$method="logVsUnsorted" curLRHiLow$significance = apply(curLRHiLow[3:(numSamples+2)], MARGIN = 1, FUN = mean); curLRHiLow$metric = apply(curLRHiLow[3:(numSamples+2)], MARGIN = 1, FUN = mean); curLRHiLow$method="logHivsLow" #compile results for A20 curResults = rbind(unique(curLRHiLow[c("method","gID","significance","metric")]), unique(curLRVsBG[c("method","gID","significance","metric")]), deseqGuideLevel[c("method","gID","significance","metric")], edgeRGuideLevel[c("method","gID","significance","metric")], unique(guideLevelStatsCast[c("method","gID","significance","metric")])) curResults = merge(curResults, unique(curCPMMat[c("gID","NT","pos")]), by="gID") curResults = curResults[!is.na(curResults$pos) | curResults$NT,] curResults$promoter = grepl("TNFAIP3", curResults$gID) | (curResults$pos <= 138189439 & curResults$pos >= 138187040) # chr6:138188077-138188379;138187040 #append the current results to all curResults$locus ="TNFAIP3" curResults$type =curtype curResults$celltype =curCelltype allResultsBoth = rbind(allResultsBoth, curResults); # (1) similarity between the effect sizes estimated per replicate, corLRHiLow = cor(curLRHiLow[!curLRHiLow$NT, 3:(3+numSamples-1)]) corLRVsBG = cor(curLRVsBG[!curLRVsBG$NT, 3:(3+numSamples-1)]) maudeZCors = cor(maudeZs[!maudeZs$NT, 4:ncol(maudeZs)]) maudeCorP=1 maudeCorR=-1 corLRHiLowP=1 corLRHiLowR=-1 corLRVsBGP=1; corLRVsBGR=-1 #select the best inter-replicate correlation for each of the three approaches for which this is possible for(i in 1:(length(sampleNames)-1)){ for(j in (i+1):length(sampleNames)){ curR = cor(maudeZs[!maudeZs$NT, sampleNames[i]], maudeZs[!maudeZs$NT, sampleNames[j]]) curP = cor.test(maudeZs[!maudeZs$NT, sampleNames[i]], maudeZs[!maudeZs$NT, sampleNames[j]])$p.value if (maudeCorR < curR){ maudeCorR = curR; maudeCorP = curP; } curR = cor(curLRVsBG[!curLRVsBG$NT, sampleNames[i]], curLRVsBG[!curLRVsBG$NT, sampleNames[j]]) curP = cor.test(curLRVsBG[!curLRVsBG$NT, sampleNames[i]], curLRVsBG[!curLRVsBG$NT, sampleNames[j]])$p.value if (corLRVsBGR < curR){ corLRVsBGR = curR; corLRVsBGP = curP; } curR = cor(curLRHiLow[!curLRHiLow$NT, sampleNames[i]], curLRHiLow[!curLRHiLow$NT, sampleNames[j]]) curP = cor.test(curLRHiLow[!curLRHiLow$NT, sampleNames[i]], curLRHiLow[!curLRHiLow$NT, sampleNames[j]])$p.value if (corLRHiLowR < curR){ corLRHiLowR = curR; corLRHiLowP = curP; } } } metricsBoth = rbind(metricsBoth, data.frame(method=c("MAUDE","logHivsLow","logVsUnsorted"), metric="r",which="replicate_correl", value=c(maudeCorR, corLRHiLowR, corLRVsBGR), sig=c(maudeCorP, corLRHiLowP, corLRVsBGP), locus ="TNFAIP3", type =curtype, celltype =curCelltype, stringsAsFactors = FALSE)) }
Here are some functions that I will make use of below.
#Some useful functions ranksumROC = function(x,y,na.rm=TRUE,...){ if (na.rm){ x=na.rm(x); y=na.rm(y); } curTest = wilcox.test(x,y,...); curTest$AUROC = (curTest$statistic/length(x))/length(y) return(curTest) } na.rm = function(x){ x[!is.na(x)]}
Evaluation
The inter-replicate correlations were calculated in the sections above and are stored in metricsBoth. Below, the two remaining metrics are calculated (distinguishing promoter vs other targeting guides, and adjacent vs randomly paired guides).
# (1) similarity between the effect sizes estimated per replicate, # (above) metricsBoth$significant= metricsBoth$sig < 0.01; # Other evaluation metrics allExpts = unique(allResultsBoth[c("celltype","locus","type")]) for (ei in 1:nrow(allExpts)){ curCelltype = allExpts$celltype[ei] curtype = allExpts$type[ei]; curLocus = allExpts$locus[ei]; curResults = allResultsBoth[allResultsBoth$celltype==curCelltype & allResultsBoth$type==curtype & allResultsBoth$locus==curLocus,] for(m in unique(curResults$method)){ curData = curResults[curResults$method==m & !curResults$NT,] # (2) similarity in effect size between adjacent guides curData = curData[order(curData$pos),] guideEffectDistances = data.frame(method = m, random=FALSE, difference = abs(curData$metric[2:nrow(curData)] - curData$metric[1:(nrow(curData)-1)]), dist =abs(curData$pos[2:nrow(curData)] - curData$pos[1:(nrow(curData)-1)]), stringsAsFactors = FALSE) guideEffectDistances = guideEffectDistances[guideEffectDistances$dist < 100,] guideEffectDistances$dist=NULL; ### changed 10 in next line to 100 to make this more robust curData = curData[sample(nrow(curData), size = nrow(curData)*100, replace = TRUE),] guideEffectDistances = rbind(guideEffectDistances, data.frame(method = m, random=TRUE, difference = abs(curData$metric[2:nrow(curData)] - curData$metric[1:(nrow(curData)-1)]), stringsAsFactors = FALSE)); # random should have more different than adjacent curRS = ranksumROC(guideEffectDistances$difference[guideEffectDistances$method==m & guideEffectDistances$random], guideEffectDistances$difference[guideEffectDistances$method==m & !guideEffectDistances$random]) metricsBoth = rbind(metricsBoth, data.frame(method=m, metric="AUROC-0.5", which="adjacent_vs_random",value = curRS$AUROC-0.5, locus=curLocus,celltype=curCelltype, type=curtype, sig=curRS$p.value, significant = curRS$p.value < 0.01)) # (3) ability to distinguish promoter-targeting guides from other guides. if (curtype=="CRISPRi" & !(m %in% c("edgeR","DESeq2"))){ # edgeR and DESeq2 are reversed for CRISPRi # non promoter should have higher effect than promoter (more -ve) curRS = ranksumROC(curResults$significance[curResults$method==m & !curResults$NT & !curResults$promoter], curResults$significance[curResults$method==m & !curResults$NT & curResults$promoter]) }else{ # promoter should have larger effect than non-promoter curRS = ranksumROC(curResults$significance[curResults$method==m & !curResults$NT & curResults$promoter], curResults$significance[curResults$method==m & !curResults$NT & !curResults$promoter]) } metricsBoth = rbind(metricsBoth, data.frame(method=m, metric="AUROC-0.5", which="promoter_vs_T", value = curRS$AUROC-0.5, locus=curLocus, celltype=curCelltype, type=curtype, sig=curRS$p.value, significant = curRS$p.value < 0.01)) } } ##compile all metrics; label the best in each test and whether any tests were significant (P<0.01) metricsBoth2 = metricsBoth; metricsBoth2$method = factor(as.character(metricsBoth2$method), levels = c("logVsUnsorted","logHivsLow","DESeq2","edgeR","MAUDE")) metricsBoth2Best = cast(metricsBoth2, which + locus + type+celltype ~ ., value="value", fun.aggregate = max) names(metricsBoth2Best)[ncol(metricsBoth2Best)] = "best" metricsBoth2 = merge(metricsBoth2, metricsBoth2Best, by = c("which","locus","type","celltype")) metricsBoth2AnySig = cast(metricsBoth2[colnames(metricsBoth2)!="value"], which + locus + type+celltype ~ ., value="significant", fun.aggregate = any) names(metricsBoth2AnySig)[ncol(metricsBoth2AnySig)] = "anySig" metricsBoth2 = merge(metricsBoth2, metricsBoth2AnySig, by = c("which","locus","type","celltype")) metricsBoth2$isBest = metricsBoth2$value==metricsBoth2$best; metricsBoth2$isBestNA = metricsBoth2$isBest; metricsBoth2$isBestNA[!metricsBoth2$isBestNA]=NA; metricsBoth2$pctOfMax = metricsBoth2$value/metricsBoth2$best * 100; #fill in NAs for edgeR and DESeq2 which cannot have inter-replicate correlations temp = metricsBoth2[metricsBoth2$metric=="r",] temp = temp[1:2,]; temp$method = c("edgeR","DESeq2") temp$value=NA; temp$isBest=NA; temp$significant=FALSE; temp$pctOfMax=NA; temp$isBestNA=NA; metricsBoth2 = rbind(metricsBoth2, temp)
Finally make the graph with all evaluation metrics as shown in the MAUDE paper.
#make the final evaluation graph p1 = ggplot(metricsBoth2[metricsBoth2$which=="adjacent_vs_random",], aes(x=method, fill=value, y=paste(locus,type,celltype))) + geom_tile() + geom_text(data=metricsBoth2[metricsBoth2$which=="adjacent_vs_random" & metricsBoth2$anySig,], aes(label="*",colour=isBestNA),show.legend = FALSE) +theme_bw() + scale_fill_gradient2(high="red", low="blue", mid="black") + theme(legend.position="top", axis.text.x = element_text(hjust=1, angle=45), axis.title.y = element_blank())+scale_colour_manual(values = c("green"), na.value=NA) + scale_y_discrete(expand=c(0,0))+scale_x_discrete(expand=c(0,0))+ggtitle("Adjacent vs\nrandom guides"); p2 = ggplot(metricsBoth2[metricsBoth2$which=="promoter_vs_T",], aes(x=method, fill=value, y=paste(locus,type,celltype))) + geom_tile() + geom_text(data=metricsBoth2[metricsBoth2$which=="promoter_vs_T" & metricsBoth2$anySig,], aes(label="*",colour=isBestNA),show.legend = FALSE) +theme_bw() + scale_fill_gradient2(high="red", low="blue", mid="black") + theme(legend.position="top", axis.text.x = element_text(hjust=1, angle=45), axis.title.y = element_blank())+scale_colour_manual(values = c("green"), na.value=NA)+ scale_y_discrete(expand=c(0,0))+scale_x_discrete(expand=c(0,0))+ ggtitle("Promoter vs other\ntargeting guides"); p3 = ggplot(metricsBoth2[metricsBoth2$which=="replicate_correl",], aes(x=method, fill=value, y=paste(locus,type,celltype))) + geom_tile() + geom_text(data=metricsBoth2[metricsBoth2$which=="replicate_correl" & metricsBoth2$anySig,], aes(label="*",colour=isBestNA),show.legend = FALSE) +theme_bw() + scale_fill_gradient2(high="red", low="blue", mid="black") + theme(legend.position="top", axis.text.x = element_text(hjust=1, angle=45), axis.title.y = element_blank())+scale_colour_manual(values = c("green"), na.value=NA)+ scale_y_discrete(expand=c(0,0))+scale_x_discrete(expand=c(0,0))+ggtitle("Replicate\ncorrelations"); g= plot_grid(p1,p2,p3, align = 'h', nrow = 1); print(g)
Session info
sessionInfo()
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.