filtCortExpr2 = cort2
#Remove outliers (quite few genes)
filtGeneLength1 = filtCortExpr1$geneLength < 10000
filtGeneLength2 = filtCortExpr2$geneLength < 10000
filtGCFull1 = (filtCortExpr1$gcFullLength >= 0.33) & (filtCortExpr1$gcFullLength <= 0.67)
filtGCFull2 = (filtCortExpr2$gcFullLength >= 0.33) & (filtCortExpr2$gcFullLength <= 0.67)
filtGCTail1 = (filtCortExpr1$gcTail >= 0.2) & (filtCortExpr1$gcTail <= 0.65)
filtGCTail2 = (filtCortExpr2$gcTail >= 0.2) & (filtCortExpr2$gcTail <= 0.65)
filtAll1 = filtGeneLength1 & filtGCFull1 & filtGCTail1
filtAll2 = filtGeneLength2 & filtGCFull2 & filtGCTail2
sum(!filtAll1) #number of discarded outliers cortex 1, 232
sum(!filtAll2) #number of discarded outliers cortex 2, 232
sum(filtAll1) #number of genes left cortex 1
sum(filtAll2) #number of genes left cortex 2
filtCort1 = filtCortExpr1[filtAll1,]
filtCort2 = filtCortExpr2[filtAll2,]
sum(is.na(filtCort1$remUMIFracOtherSample))
sum(is.na(filtCort2$remUMIFracOtherSample))
#experiment with not using umi copy fraction for genes with very few counts - it seems 5 UMIs or less gives only noise
filtCortRemRUFExpr1 = filtCortExpr1;
filtCortRemRUFExpr1$remUMIFracOtherSample[filtCortExpr2$UMI < 6] = NA
#filtCortRemRUFExpr1$remUMIFracOtherSample[filtCortExpr1$UMITPM < 0.001 | filtCortExpr1$bulkTPM < 0.001] = NA
filtCortRemRUF1 = filtCortRemRUFExpr1[filtAll1,]
#resCort1RemUMIFracSpec = genScToBulkCovGraphs(filtCortRemRUF1, logUMITMM ~ remUMIFracOtherSample, LogUMIDivBulk ~ remUMIFracOtherSample, 7, "UMI copy fraction")
filtCortRemRUFExpr2 = filtCortExpr2;
filtCortRemRUFExpr2$remUMIFracOtherSample[filtCortExpr1$UMI < 6] = NA
#filtCortRemRUFExpr1$remUMIFracOtherSample[filtCortExpr1$UMITPM < 0.001 | filtCortExpr1$bulkTPM < 0.001] = NA
filtCortRemRUF2 = filtCortRemRUFExpr2[filtAll2,]
#resCort2RemUMIFracSpec = genScToBulkCovGraphs(filtCortRemRUF2, logUMITMM ~ remUMIFracOtherSample, LogUMIDivBulk ~ remUMIFracOtherSample, 7, "UMI copy fraction")
#0 Copies per UMI
###########################
#resCort1CopiesPerUMI = genScToBulkCovGraphs(filtCort1, logUMITMM ~ CopiesPerUMIOtherSample, LogUMIDivBulk ~ CopiesPerUMIOtherSample, 14, "Copies per UMI")
#resCort2RemUMIFrac = genScToBulkCovGraphs(filtCort2, logUMITMM ~ remUMIFracOtherSample, LogUMIDivBulk ~ remUMIFracOtherSample, 7, "UMI copy fraction")
#1 Removed counts' fraction
###########################
#resCort1RemUMIFrac = genScToBulkCovGraphs(filtCort1, logUMITMM ~ remUMIFracOtherSample, LogUMIDivBulk ~ remUMIFracOtherSample, 7, "UMI copy fraction")
#resCort2RemUMIFrac = genScToBulkCovGraphs(filtCort2, logUMITMM ~ remUMIFracOtherSample, LogUMIDivBulk ~ remUMIFracOtherSample, 7, "UMI copy fraction")
#use the data where UMI copy fraction is removed where it is based on too few molecules
resCort1RemUMIFrac = genScToBulkCovGraphs(filtCortRemRUF1, logUMITMM ~ remUMIFracOtherSample, LogUMIDivBulk ~ remUMIFracOtherSample, 7, "UMI copy fraction")
resCort2RemUMIFrac = genScToBulkCovGraphs(filtCortRemRUF2, logUMITMM ~ remUMIFracOtherSample, LogUMIDivBulk ~ remUMIFracOtherSample, 7, "UMI copy fraction")
head(filtCortRemRUF1)
resCort1RemUMIFrac2 = genScToBulkCovGraphs(filtCortRemRUF1, logUMITMM ~ remUMIFrac, LogUMIDivBulk ~ remUMIFrac, 6, "UMI copy fraction")
#EVAL no corr
tr2gMouse <- transcript2gene(species = "Mus musculus", type = "vertebrate",
kallisto_out_path = paste0(dataPath, "EVAL/bus_output"),
ensembl_version = 96, verbose=T)
busFileEVALNoCorr = read.table(paste0(dataPath,"EVAL/evalnocorr.txt"))
resEVALNoCorr = plotCopiesPerGene(busFileEVALNoCorr, whitelist_v2, paste0(dataPath,"EVAL/bus_output/"), tr2gMouse)
saveRDS(resEVALNoCorr, paste0(cachedDataPath,"resEVALNoCorr.Rds"))
install.packages("bootstrap")
setwd("C:\Work\R\RNASeqCTProfileEval")
setwd("C:/Work/R/RNASeqCTProfileEval")
load("saves/2019-12-19-rf-bulk-lab-5-celltype-tissue.RData")
load("saves/2019-12-19-rf-bulk-lab-5-only-b-subcelltype-tissue.RData")
load("saves/2019-12-19-rf-bulk-lab-5-only-t-subcelltype-tissue.RData")
load("saves/2019-12-19-rf-bulk-celltype-lab.RData")
load("saves/2019-12-19-rf-sc-celltype-lab.RData")
load("saves/2019-12-19-rf-lab-5-6-sc_or_bulk-celltype.RData")
create_subplot <- function(res, var_groups, title, ylab = FALSE) {
oob_mean_var_groups_imp <- lapply(
res, function(r) colMeans(r$oob_var_groups_imp))
p <- tibble(
variable = factor(rep(var_groups, times = 3), levels = var_groups),
var_imp = do.call(c,
lapply(oob_mean_var_groups_imp, function(vi) vi / sum(vi))),
norm = factor(
rep(c("TPM", "TMM", "Quantile"), each = length(var_groups)),
levels = c("TPM", "TMM", "Quantile"))) %>%
ggplot() +
geom_bar(
aes(x = variable, y = var_imp, fill = norm),
stat = "identity",
position = "dodge") +
labs(x = NULL, y = NULL, title = title) +
scale_y_continuous(limits = c(0, 1)) +
scale_fill_manual("Normalisation", values = cbPalette) +
theme_minimal() +
theme(plot.title = element_text(size = 10, face = "bold"))
if (ylab) {
p <- p +
labs(x = NULL, y = "Variable Importance (in %)") +
theme(plot.margin = unit(c(0.2, 0, 0.2, 0.55), "cm"),
axis.title = element_text(size = 8))
} else {
p <- p + theme(plot.margin = unit(c(0.2, 0, 0.2, 1), "cm"))
}
# Return final plot
p
}
p1 <- create_subplot(
res_sc_or_bulk_celltype_lab_5_6,
c("SC vs Bulk", "Celltype"),
"Tissue 1 and labs (5, 6)")
# Colour-blind friendly palette
cbPalette <- c(
"#999999", "#E69F00", "#56B4E9", "#009E73",
"#F0E442", "#0072B2", "#D55E00", "#CC79A7")
p1 <- create_subplot(
res_sc_or_bulk_celltype_lab_5_6,
c("SC vs Bulk", "Celltype"),
"Tissue 1 and labs (5, 6)")
p2 <- create_subplot(
res_celltype_tissue_lab_5_bulk,
c("Celltype", "Tissue"),
"Lab 5 and bulk samples")
p3 <- create_subplot(
res_subcelltype_tissue_lab_5_bulk_only_b,
c("Sub-celltype", "Tissue"),
"Lab 5, bulk samples, and B cells",
ylab = TRUE)
p4 <- create_subplot(
res_subcelltype_tissue_lab_5_bulk_only_t,
c("Sub-celltype", "Tissue"),
"Lab 5, bulk samples, and T cells")
p5 <- create_subplot(
res_celltype_lab_bulk,
c("Celltype", "Lab"),
"Tissue 1, mixed sub-celltype and bulk samples")
p6 <- create_subplot(
res_celltype_lab_sc,
c("Celltype", "Lab"),
"Tissue 1, mixed sub-celltype and SC samples")
ggpubr::ggarrange(
p1, p2, p3, p4, p5, p6,
ncol = 2, nrow = 3,
common.legend = TRUE, legend = "bottom",
labels = c("A", "B", "C", "D", "E", "F"),
label.x = 0.03)
p1 <- create_subplot(
res_sc_or_bulk_celltype_lab_5_6,
c("SC vs Bulk", "Celltype"),
"Samples from Blood, Lab Bulk5 and SC HCA CB")
ggpubr::ggarrange(
p1, p2, p3, p4, p5, p6,
ncol = 2, nrow = 3,
common.legend = TRUE, legend = "bottom",
labels = c("A", "B", "C", "D", "E", "F"),
label.x = 0.03)
p1 <- create_subplot(
res_sc_or_bulk_celltype_lab_5_6,
c("SC vs Bulk", "Celltype"),
"Lab Bulk5 and SC HCA CB, Samples from Blood")
ggpubr::ggarrange(
p1, p2, p3, p4, p5, p6,
ncol = 2, nrow = 3,
common.legend = TRUE, legend = "bottom",
labels = c("A", "B", "C", "D", "E", "F"),
label.x = 0.03)
source('C:/Work/R/RNASeqCTProfileEval/GenFigSourcesOfVarPairwise.R', echo=TRUE)
boxFac = factor(boxes, 1:15, c("Techn. repl. bulk 3",
"Diff. indiv. bulk 3",
"Same indiv. bulk 4",
"Diff. indiv. bulk 4",
"Diff. Tissue bulk 5",
"Cell subtype bulk 5",
"Diff. indiv. bulk 5",
"CT bulk 5",
"CT Bulk",
"CT Sc",
"Lab Bulk",
"Lab Sc",
"Sc vs Bulk",
"Bulk diff. indiv.",
"Sc diff. indiv."))
p = boxPlotWithDots(joinedData, boxFac, 1)
p<-p + facet_wrap( ~ boxes, nrow=1, strip.position = "bottom") +
theme(strip.text.x = element_text(size=10, angle=90, hjust=1),
strip.background = element_rect(colour="transparent", fill="transparent"))
p
boxFac = factor(boxes, 1:15, c("Techn. repl. Bulk3",
"Diff. indiv. Bulk3",
"Same indiv. Bulk4",
"Diff. indiv. Bulk4",
"Diff. Tissue Bulk5",
"Cell subtype Bulk5",
"Diff. indiv. Bulk5",
"CT Bulk5",
"CT bulk",
"CT sc",
"Lab bulk",
"Lab sc",
"Sc vs bulk",
"Bulk diff. indiv.",
"Sc diff. indiv."))
p = boxPlotWithDots(joinedData, boxFac, 1)
p<-p + facet_wrap( ~ boxes, nrow=1, strip.position = "bottom") +
theme(strip.text.x = element_text(size=10, angle=90, hjust=1),
strip.background = element_rect(colour="transparent", fill="transparent"))
p
boxFac = factor(boxes, 1:15, c("Techn. repl. Bulk3",
"Diff. indiv. Bulk3",
"Same indiv. Bulk4",
"Diff. indiv. Bulk4",
"Diff. Tissue Bulk5",
"Cell subtype Bulk5",
"Diff. indiv. Bulk5",
"Cell type Bulk5",
"Cell type all bulk",
"Cell type all sc",
"Lab all bulk",
"Lab all sc",
"Sc vs bulk",
"Bulk diff. indiv.",
"Sc diff. indiv."))
p = boxPlotWithDots(joinedData, boxFac, 1)
p<-p + facet_wrap( ~ boxes, nrow=1, strip.position = "bottom") +
theme(strip.text.x = element_text(size=10, angle=90, hjust=1),
strip.background = element_rect(colour="transparent", fill="transparent"))
p
ggpubr::ggarrange(
p1, p2, p3, p4, p5, p6,
ncol = 2, nrow = 3,
common.legend = TRUE, legend = "bottom",
labels = c("A", "B", "C", "D", "E", "F"),
label.x = 0.03)
p1 <- create_subplot(
res_sc_or_bulk_celltype_lab_5_6,
c("Mix of sc and bulk", "Cell type"),
"Mix of Sc and Bulk vs Cell Type")
p2 <- create_subplot(
res_celltype_tissue_lab_5_bulk,
c("Cell type", "Tissue"),
"Cell Type vs Tissue")
p3 <- create_subplot(
res_subcelltype_tissue_lab_5_bulk_only_b,
c("Cell subtype", "Tissue"),
"Cell Subtype vs Tissue, B cells",
ylab = TRUE)
p4 <- create_subplot(
res_subcelltype_tissue_lab_5_bulk_only_t,
c("Cell subtype", "Tissue"),
"Cell Subtype vs Tissue, T cells")
p5 <- create_subplot(
res_celltype_lab_bulk,
c("Cell type", "Lab"),
"Cell Type vs Lab, Bulk")
p6 <- create_subplot(
res_celltype_lab_sc,
c("Cell type", "Lab"),
"Cell Type vs Lab, Sc")
ggpubr::ggarrange(
p1, p2, p3, p4, p5, p6,
ncol = 2, nrow = 3,
common.legend = TRUE, legend = "bottom",
labels = c("A", "B", "C", "D", "E", "F"),
label.x = 0.03)
ggpubr::ggarrange(
p1, p2, p3, p4, p5, p6,
ncol = 2, nrow = 3,
common.legend = TRUE, legend = "bottom",
labels = c("A", "B", "C", "D", "E", "F"),
label.x = 0.03)
ggsave("plots/rf-var-group-importance.png",
width = 20.32, height = 18, unit = "cm", dpi = 100)
ggsave("plots/rf-var-group-importance.png",
width = 20.32, height = 18, unit = "cm", dpi = 100)
p1 <- create_subplot(
res_sc_or_bulk_celltype_lab_5_6,
c("Mix of sc and bulk", "Cell type"),
"Mix of Single-Cell and Bulk vs Cell Type")
p6 <- create_subplot(
res_celltype_lab_sc,
c("Cell type", "Lab"),
"Cell Type vs Lab, Single-Cell")
ggpubr::ggarrange(
p1, p2, p3, p4, p5, p6,
ncol = 2, nrow = 3,
common.legend = TRUE, legend = "bottom",
labels = c("A", "B", "C", "D", "E", "F"),
label.x = 0.03)
ggsave("plots/rf-var-group-importance.png",
width = 20.32, height = 18, unit = "cm", dpi = 100)
p = boxPlotWithDots(joinedData, boxFac, 1)
p
ggpubr::ggarrange(
p1, p2, p3, p4, p5, p6,
ncol = 2, nrow = 3,
common.legend = TRUE, legend = "bottom", title="Test"
labels = c("A", "B", "C", "D", "E", "F"),
label.x = 0.03)
ggpubr::ggarrange(
p1, p2, p3, p4, p5, p6,
ncol = 2, nrow = 3,
common.legend = TRUE, legend = "bottom", title="Test",
labels = c("A", "B", "C", "D", "E", "F"),
label.x = 0.03)
ggpubr::ggarrange(
p1, p2, p3, p4, p5, p6,
ncol = 2, nrow = 3,
common.legend = TRUE, legend = "bottom", title="Test",
labels = c("A", "B", "C", "D", "E", "F"),
label.x = 0.03)
fig4 = ggpubr::ggarrange(
p1, p2, p3, p4, p5, p6,
ncol = 2, nrow = 3,
common.legend = TRUE, legend = "bottom", title="Test",
labels = c("A", "B", "C", "D", "E", "F"),
label.x = 0.03)
fig4 = ggpubr::ggarrange(
p1, p2, p3, p4, p5, p6,
ncol = 2, nrow = 3,
common.legend = TRUE, legend = "bottom",
labels = c("A", "B", "C", "D", "E", "F"),
label.x = 0.03)
#check that the title is shown on the graph, it sometimes randomly disappears.
annotate_figure(fig4,
top = text_grob("Estimating Variation Factors Using Random Forest Regression", face = "bold", size = 14))
library(ggpubr)
#check that the title is shown on the graph, it sometimes randomly disappears.
annotate_figure(fig4,
top = text_grob("Estimating Variation Factors Using Random Forest Regression", face = "bold", size = 14))
source('C:/Work/R/RNASeqCTProfileEval/GenFigScVsBulk.R', echo=TRUE)
fig1 = ggplot(data= df , aes(x=plotSampleOrder, y=value, fill=Dataset)) +
geom_boxplot(outlier.shape = NA, coef = 0) +
labs(y="Relative Log Expression", x="Samples") +
scale_fill_manual(values=colors) +
coord_cartesian(ylim=c(-2, 2)) +
theme( axis.text.x=element_blank(), axis.ticks.x=element_blank()) +
geom_line(data = dfline, aes(x=xline, y=yline, fill=NA)) +
facet_wrap( ~ nm, nrow=3)
#Run GenerateMixedDataMatrices before running this file!
library("ggplot2")
library("Seurat")
dataFolder = "C:/Work/R/RNASeqCTProfileEval/"
source(paste0(dataFolder, "FigureHelpFunc.R"))
#############################
# load data
#############################
library(dplyr)
library(tibble)
tpmScAndBulk = readRDS(paste0(dataFolder, "data/tpmScAndBulk.RDS"))
tmmScAndBulk = readRDS(paste0(dataFolder, "data/tmmScAndBulk.RDS"))
quantileScAndBulk = readRDS(paste0(dataFolder, "data/quantileScAndBulk.RDS"))
bcScAndBulk = readRDS(paste0(dataFolder, "data/bcScAndBulk.RDS"))
cellTypes = readRDS(paste0(dataFolder, "data/cellTypes.RDS"))
labs = readRDS(paste0(dataFolder, "data/labs.RDS"))
subCellTypes = readRDS(paste0(dataFolder, "data/subCellTypes.RDS"))
tissues = readRDS(paste0(dataFolder, "data/tissues.RDS"))
###########################
## Figure 1
###########################
#for comparison:
#plotRLE(as.matrix(bulk_counts), outline=FALSE, ylim=c(-4, 4), col=as.numeric(cellTypes))
#plotRLE(as.matrix(bulk_tpm), outline=FALSE, ylim=c(-4, 4), col=as.numeric(cellTypes))
#plotRLE(as.matrix(bulk_tmm), outline=FALSE, ylim=c(-4, 4), col=as.numeric(cellTypes))
#so, instead of plotRLE we want to plot this using ggplot so we can put all of them together in a figure
#first: prepare data
library("matrixStats")
genRLEData <- function(df) {
df = as.matrix(df)
med = rowMedians(df)
output = df
for (i in 1:dim(output)[2]) {
output[,i] = log2((df[,i] + 1) / (med + 1))
}
return (output)
}
a = genRLEData(tpmScAndBulk)
b = genRLEData(tmmScAndBulk)
c = genRLEData(quantileScAndBulk)
library(dplyr)
library(tidyr)
labvals = c(3,2,4,5,1,6,7,8,9,10)
labnames = c("Bulk 1", "Bulk 2", "Bulk 3", "Bulk 4", "Bulk 5", "SC HCA CB", "SC LC", "SC PBMC68k", "SC Mixed 10x data", "SC Melanoma")
plotSampleOrder = c(2:7,1,8:12,25,13:24,26:64,65:74,75:90,91:98,104:105,101:103,99:100)
plotSampleOrderFac = factor(1:length(plotSampleOrder), levels=plotSampleOrder)
#now with ggplot:
dfTPM <- a %>%
as.data.frame %>%
gather(., key= sample, value="value") %>%
mutate(Dataset = rep(factor(labs, labvals, labnames),each=dim(a)[1]),
plotSampleOrder = rep(plotSampleOrderFac, each=dim(a)[1]))
dfTMM <- b %>%
as.data.frame %>%
gather(., key= sample, value="value") %>%
mutate(Dataset = rep(factor(labs, labvals, labnames),each=dim(b)[1]),
plotSampleOrder = rep(plotSampleOrderFac, each=dim(b)[1]))
dfQ <- c %>%
as.data.frame %>%
gather(., key= sample, value="value") %>%
mutate(Dataset = rep(factor(labs, labvals, labnames),each=dim(c)[1]),
plotSampleOrder = rep(plotSampleOrderFac, each=dim(c)[1]))
xline = c(0,dim(a)[2])
yline = c(0,0)
dfline = data.frame(xline,yline)
#Combine all three normalization methods:
df = rbind(dfTPM,dfTMM, dfQ)
len = dim(dfTPM)[1]
nm = rep(c(1,2,3), each=len)
nm = factor(nm, c(1,2,3), c("TPM/CPM", "TMM", "Quantile"))
df = cbind(df,nm)
colors = c("#a6cee3", "#1f78b4", "#b2df8a", "#33a02c", "#fb9a99",
"#e31a1c", "#fdbf6f", "#ff7f00", "#cab2d6", "#6a3d9a") #created with colorbrewer
fig1 = ggplot(data= df , aes(x=plotSampleOrder, y=value, fill=Dataset)) +
geom_boxplot(outlier.shape = NA, coef = 0) +
labs(y="Relative Log Expression", x="Samples") +
scale_fill_manual(values=colors) +
coord_cartesian(ylim=c(-2, 2)) +
theme( axis.text.x=element_blank(), axis.ticks.x=element_blank()) +
geom_line(data = dfline, aes(x=xline, y=yline, fill=NA)) +
facet_wrap( ~ nm, nrow=3)
#check that the title is shown on the graph, it sometimes randomly disappears.
annotate_figure(fig1,
top = text_grob("Relative Log Expression After Normalization", face = "bold", size = 14))
#check that the title is shown on the graph, it sometimes randomly disappears.
annotate_figure(fig1,
top = text_grob("Relative Log Expression After Normalization", face = "bold", size = 14))
a2 = prcomp(t(log2(tpmScAndBulk + 1)))
b2 = prcomp(t(log2(tmmScAndBulk + 1)))
c2 = prcomp(t(log2(quantileScAndBulk + 1)))
d2 = prcomp(t(log2(bcScAndBulk + 1)))
#calculate the explained variance for principle component 1 and 2 for each
explVar <- function(pcaData, pc) {
eigs <- pcaData$sdev^2
return (eigs[pc] / sum(eigs))
}
x=c(1,1,1,1)
explPC1=c(explVar(a2, 1), explVar(b2, 1), explVar(c2, 1), explVar(d2, 1))
explPC2=c(explVar(a2, 2), explVar(b2, 2), explVar(c2, 2), explVar(d2, 2))
textPC1=c(paste0(round(100*explPC1[1],2), "% - Technical"),
paste0(round(100*explPC1[2],2), "% - Technical"),
paste0(round(100*explPC1[3],2), "% - Technical"),
paste0(round(100*explPC1[4],2), "% - Cell type"))
textPC2=c(paste0(round(100*explPC2[1],2), "% - Cell type"),
paste0(round(100*explPC2[2],2), "% - Cell type"),
paste0(round(100*explPC2[3],2), "% - Cell type"),
paste0(round(100*explPC2[4],2), "% - Unknown"))
nm = factor(c(1,2,3,4), c(1,2,3,4), c("TPM/CPM", "TMM", "Quantile", "TMM + ComBat"))
dat2 = data.frame(x,explPC1,explPC2,textPC1,nm)
summary(d2) #check that the PC1 explained variance fraction was correctly calculated -> OK!
#scale all PCs to be between 0 and 1
scalePCs <- function(d) {
for (i in 1:(dim(d)[2])) {
minval = min(d[,i])
maxval = max(d[,i])
d[,i] = (d[,i] - minval) / (maxval - minval)
}
return (d)
}
a2$x = scalePCs(a2$x)
b2$x = scalePCs(b2$x)
c2$x = scalePCs(c2$x)
d2$x = scalePCs(d2$x)
#flip PC2 upside down for c so the graphs are more similar. Same for PC1 for b
c2$x[,2] = -c2$x[,2] + 1
b2$x[,1] = -b2$x[,1] + 1
#plot(a2$x[,1], a2$x[,2])
dfTPM <- as.data.frame(a2$x)
dfTPM$labs <- factor(labs, labvals, labnames)
dfTPM$bulkAndCt <- factor((cellTypes - 1) + (labs > 5) * 2, c(0,1,2,3), c("Bulk B Cell", "Bulk T Cell", "Sc B Cell", "Sc T Cell"))
dfTPM$pc1Expl = explVar(a2, 1)
dfTPM$pc2Expl = explVar(a2, 2)
dfTMM <- as.data.frame(b2$x)
dfTMM$labs <- factor(labs, labvals, labnames)
dfTMM$bulkAndCt <- factor((cellTypes - 1) + (labs > 5) * 2, c(0,1,2,3), c("Bulk B Cell", "Bulk T Cell", "Sc B Cell", "Sc T Cell"))
dfTMM$pc1Expl = explVar(b2, 1)
dfTMM$pc2Expl = explVar(b2, 2)
dfQ <- as.data.frame(c2$x)
dfQ$labs <- factor(labs, labvals, labnames)
dfQ$bulkAndCt <- factor((cellTypes - 1) + (labs > 5) * 2, c(0,1,2,3), c("Bulk B Cell", "Bulk T Cell", "Sc B Cell", "Sc T Cell"))
dfQ$pc1Expl = explVar(c2, 1)
dfQ$pc2Expl = explVar(c2, 2)
dfBC <- as.data.frame(d2$x)
dfBC$labs <- factor(labs, labvals, labnames)
dfBC$bulkAndCt <- factor((cellTypes - 1) + (labs > 5) * 2, c(0,1,2,3), c("Bulk B Cell", "Bulk T Cell", "Sc B Cell", "Sc T Cell"))
dfBC$pc1Expl = explVar(d2, 1)
dfBC$pc2Expl = explVar(d2, 2)
#Combine all methods:
df = rbind(dfTPM,dfTMM, dfQ, dfBC)
len = dim(dfTPM)[1]
nm = rep(c(1,2,3,4), each=len)
nm = factor(nm, c(1,2,3,4), c("TPM/CPM", "TMM", "Quantile", "TMM + ComBat"))
df = cbind(df,nm)
colors = c("#a6cee3", "#1f78b4", "#b2df8a", "#33a02c", "#fb9a99",
"#e31a1c", "#fdbf6f", "#ff7f00", "#cab2d6", "#6a3d9a") #created with colorbrewer
p<-ggplot(df,aes(x=PC1,y=PC2,color=labs, shape=bulkAndCt))
p<-p + geom_point()
p<-p + scale_color_manual(values=colors)
p$labels$shape = "Sample properties"
p$labels$colour = "Dataset"
p<-p + coord_cartesian(xlim=c(-0.09, 1), ylim=c(-0.09, 1))#create room for PC label
p<-p + theme( axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank())
p<-p + geom_rect(data=dat2, mapping=aes(xmin=-0.114, xmax=-0.04, ymin=0, ymax=explPC2), color=NA, fill="blue", alpha=0.2, inherit.aes = FALSE)
p<-p + geom_rect(data=dat2, mapping=aes(xmin=-0.114, xmax=-0.04, ymin=explPC2, ymax=1), color=NA, alpha=0.0, inherit.aes = FALSE)
p<-p + geom_rect(data=dat2, mapping=aes(xmin=-0.114, xmax=-0.04, ymin=0, ymax=1), color="blue", fill=NA, inherit.aes = FALSE)
p<-p + geom_rect(data=dat2, mapping=aes(ymin=-0.12, ymax=-0.04, xmin=0, xmax=explPC1), color=NA, fill="blue", alpha=0.2, inherit.aes = FALSE)
p<-p + geom_rect(data=dat2, mapping=aes(ymin=-0.12, ymax=-0.04, xmin=explPC1, xmax=1), color=NA, alpha=0.0, inherit.aes = FALSE)
p<-p + geom_rect(data=dat2, mapping=aes(ymin=-0.12, ymax=-0.04, xmin=0, xmax=1), color="blue", fill=NA, inherit.aes = FALSE)
p<-p + geom_text(data=dat2, mapping=aes(y=-0.075, x=0.5, label=textPC1), size=3.8, color="black", inherit.aes = FALSE)
p<-p + geom_text(data=dat2, mapping=aes(x=-0.085, y=0.5, label=textPC2), size=3.8, color="black", angle = 90, inherit.aes = FALSE)
p<-p + facet_wrap( ~ nm, nrow=2)
annotate_figure(p,
top = text_grob("Visualization of Batch Effects", face = "bold", size = 14))
annotate_figure(p,
top = text_grob("Visualization of Batch Effects", face = "bold", size = 14))
p<-ggplot(df,aes(x=PC1,y=PC2,color=labs, shape=bulkAndCt))
p<-p + geom_point()
p<-p + scale_color_manual(values=colors)
p$labels$shape = "Sample properties"
p$labels$colour = "Dataset"
p<-p + coord_cartesian(xlim=c(-0.09, 1), ylim=c(-0.09, 1))#create room for PC label
p<-p + theme( axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank())
p<-p + geom_rect(data=dat2, mapping=aes(xmin=-0.114, xmax=-0.04, ymin=0, ymax=explPC2), color=NA, fill="blue", alpha=0.2, inherit.aes = FALSE)
p<-p + geom_rect(data=dat2, mapping=aes(xmin=-0.114, xmax=-0.04, ymin=explPC2, ymax=1), color=NA, alpha=0.0, inherit.aes = FALSE)
p<-p + geom_rect(data=dat2, mapping=aes(xmin=-0.114, xmax=-0.04, ymin=0, ymax=1), color="blue", fill=NA, inherit.aes = FALSE)
p<-p + geom_rect(data=dat2, mapping=aes(ymin=-0.12, ymax=-0.04, xmin=0, xmax=explPC1), color=NA, fill="blue", alpha=0.2, inherit.aes = FALSE)
p<-p + geom_rect(data=dat2, mapping=aes(ymin=-0.12, ymax=-0.04, xmin=explPC1, xmax=1), color=NA, alpha=0.0, inherit.aes = FALSE)
p<-p + geom_rect(data=dat2, mapping=aes(ymin=-0.12, ymax=-0.04, xmin=0, xmax=1), color="blue", fill=NA, inherit.aes = FALSE)
p<-p + geom_text(data=dat2, mapping=aes(y=-0.075, x=0.5, label=textPC1), size=3.8, color="black", inherit.aes = FALSE)
p<-p + geom_text(data=dat2, mapping=aes(x=-0.085, y=0.5, label=textPC2), size=3.8, color="black", angle = 90, inherit.aes = FALSE)
p<-p + facet_wrap( ~ nm, nrow=2)
p
