LoadExonIntronCounts <- function(inf){
# inf <- "/home/shared/4c_seq/bowtie2_WT_Cry1Intron/Liver_Liv06_S6_R1_filtered.counts_mat"
bname <- basename(inf)
jtiss <- strsplit(bname, split = "_")[[1]][[1]]
if (jtiss == "Kidney"){
# Kidney_ZT8_k_42_filtered.counts_mat
samp.numb <- as.numeric(gsub("[^\\d]+", "", strsplit(bname, split = "_")[[1]][[4]], perl=TRUE))
jtime <- as.numeric(gsub("ZT", "", strsplit(bname, split = "_")[[1]][[2]]))
jgeno.sym <- strsplit(bname, split = "_")[[1]][[3]] # k or w
if (jgeno.sym == "k"){
jgeno <- "Cry1IntronKO"
} else if (jgeno.sym == "w"){
jgeno <- "WT"
} else {
stop(paste("jgeno must be k or w", jgeno.sym))
}
} else if (jtiss == "Liver") {
# Liver_Liv09_S9_R1_filtered.counts_mat
samp.numb <- as.numeric(gsub("[^\\d]+", "", strsplit(bname, split = "_")[[1]][[2]], perl=TRUE))
jtime <- SeriesToTime(samp.numb)
jgeno <- SeriesToGeno(samp.numb)
} else {
stop(paste("jtiss must be Liver or Kidney", jtiss))
}
dat <- read.table(inf, header = TRUE, sep = "\t", stringsAsFactors=FALSE)
dat$samp.numb <- samp.numb
dat$tissue <- jtiss
dat$geno <- jgeno
dat$time <- jtime
return(dat)
}
PhaseDelta <- function(p.vec){
if (abs(diff(p.vec)) > 12){
# take smallest p and add 24
p.vec[which.min(p.vec)] <- p.vec[which.min(p.vec)] + 24
}
# take 2nd minus 1st, as in diff function
return(diff(p.vec))
}
SeriesToTime <- function(samp.numb){
samp.numb.by.12 <- (samp.numb - 1) %% 12
# odds are WT, evens are KO, timepoints 0, 4, 8 ...
jtimes <- ifelse(samp.numb %% 2, (floor(samp.numb.by.12/2)) * 4, (floor(samp.numb.by.12/2)) * 4) # 0,0,4,4,8,8 ... 68,68
}
SeriesToGeno <- function(samp.numb){
jgenos <- ifelse(samp.numb %% 2 == 1, "WT", "Cry1IntronKO")
}
LoadDataset <- function(jtissue, jeps = 1e-3){
inf <- paste0("/home/shared/4c_seq/kallisto_WT_Cry1Intron/Mermet_", jtissue, "_Series_1_WT_Cry1IntronKO.txt")
dat <- read.table(inf, header = TRUE, sep = "\t")
# Ensembl to gene ID ------------------------------------------------------
dat$gene <- Transcript2Gene(dat$target_id, return.original = TRUE)
# If gene name all numbers, then use transcript name
dat$gene[grepl("^[1-9]*$", dat$gene)] <- dat$target_id[grepl("^[1-9]*$", dat$gene)]
colnames(dat) <- gsub("Liv", "Liver", colnames(dat))
dat.exprs <- dat[, grepl(paste0("^", jtissue), colnames(dat))]
sampnames <- colnames(dat)[grepl(paste0("^", jtissue), colnames(dat))]
if (jtissue == "Liver"){
# colnames need to be readjusted to be "Tissue_Time_Geno"
samp.numb <- as.numeric(sapply(sampnames, function(s) substr(s, start=6, stop = 7), USE.NAMES = TRUE))
# odds are WT, evens are KO, timepoints 0, 4, 7 ...
# times
samp.numb.by.12 <- (samp.numb - 1) %% 12
jtimes <- ifelse(samp.numb %% 2, (floor(samp.numb.by.12/2)) * 4, (floor(samp.numb.by.12/2)) * 4) # 0,0,4,4,8,8 ... 68,68
jgenos <- ifelse(samp.numb %% 2 == 1, "WT", "Cry1IntronKO")
jtissues <- rep("Liver", length(jtimes))
sampnames <- make.names(paste(jtissues, jtimes, jgenos, sep = "_"), unique=TRUE)
}
# fix sampnames if necessary
if (any(grepl("\\.", sampnames))){
psplit <- sapply(sampnames, function(s) strsplit(s, "_")[[1]][[3]], simplify = TRUE, USE.NAMES = FALSE)
# set biorep to 1 if no period, otherwise take number after "." and add 1
biorep <- as.numeric(sapply(psplit, function(p){
tryCatch(as.numeric(strsplit(p, "\\.")[[1]][[2]]) + 1, error = function(e) 1)
}, simplify = TRUE, USE.NAMES = FALSE))
# fix sampnames
samp.fix.i <- grepl("\\.[1-9]$", sampnames)
sampnames[samp.fix.i] <- substr(sampnames[samp.fix.i], start = 1, stop = nchar(sampnames[samp.fix.i]) - 2)
} else {
biorep <- rep(1, length(sampnames))
}
tissues <- sapply(sampnames, function(s) strsplit(s, "_")[[1]][[1]], USE.NAMES = FALSE)
times <- sapply(sampnames, function(s) as.numeric(gsub("ZT", "", strsplit(s, "_")[[1]][[2]])), USE.NAMES = FALSE)
genos <- sapply(sampnames, function(s) strsplit(s, "_")[[1]][[3]], USE.NAMES = FALSE)
dat.long <- data.frame(gene = dat$gene,
tissue = rep(tissues, each = nrow(dat.exprs)),
time = rep(times, each = nrow(dat.exprs)),
geno = rep(genos, each = nrow(dat.exprs)),
biorep = rep(biorep, each = nrow(dat.exprs)),
sampnames = rep(sampnames, each = nrow(dat.exprs)),
exprs = unlist(dat.exprs))
dat.long$gene <- as.character(dat.long$gene)
dat.long$geno <- factor(dat.long$geno, levels = c("WT", "Cry1IntronKO"))
dat.long$tissgeno <- paste(dat.long$tissue, dat.long$geno, sep = "_")
dat.long$tissue <- dat.long$tissgeno
# fix time to plot over entire series
dat.long$time <- mapply(function(jtime, jbiorep) jtime + (jbiorep - 1) * 24, dat.long$time, dat.long$biorep)
dat.long <- dat.long %>%
group_by(gene, tissue, time, geno) %>%
summarise(exprs = sum(exprs))
# jeps <- 1e-3
jcutoff <- -2
ggplot(dat.long, aes(x = log2(exprs + jeps))) + geom_density() + geom_vline(xintercept = jcutoff)
# log2 and tpm separate
dat.long$tpm <- dat.long$exprs
dat.long$exprs <- log2(dat.long$exprs + jeps)
dat.long$tissue.orig <- dat.long$tissue
dat.long$tissue <- paste(dat.long$tissue.orig, dat.long$geno, sep = "_")
return(dat.long)
}
PlotRnaseq <- function(dat.long.merged.sub, single.day=FALSE, time.space = 8, jtitle="", linear.eps = FALSE){
if (linear.eps != FALSE){
if (is.na(as.numeric(linear.eps))) stop("linear.eps must be numeric if not FALSE")
dat.long.merged.sub$exprs <- dat.long.merged.sub$tpm
jylab <- "TPM"
} else {
jylab <- "Log2(TPM)"
}
if (jtitle == ""){
jtitle <- dat.long.merged.sub$gene[[1]]
}
if (single.day){
source("/home/yeung/projects/tissue-specificity/scripts/functions/PlotGeneAcrossTissues.R")
dat.long.merged.sub <- dat.long.merged.sub %>%
ungroup() %>%
mutate(experiment = "RNASeq")
dat.long.merged.sub <- ConvertToSingleDay(dat.long.merged.sub)
}
m <- ggplot(dat.long.merged.sub, aes(x = time, y = exprs, linetype = geno, colour = tissue)) + geom_point() + geom_line() +
theme_bw(24) + xlab("ZT") + ylab(jylab) +
theme(aspect.ratio=1, panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
ggtitle(jtitle) +
scale_x_continuous(breaks = seq(0, max(dat.long.merged.sub$time), time.space))
if (single.day){
m <- m + geom_errorbar(data = dat.long.merged.sub, aes(ymin=exprs-sem, ymax=exprs+sem, linetype = geno, colour = tissue), linetype = "solid", size=0.5, width=0.5)
}
return(m)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.