.Rhistory

metadf
bakRDataobj <- bakRData(cB, metadf)
# Assess raw mutation rates
cB %>% group_by(sample) %>%
summarise(mutrate = sum(TC*n)/sum(nT*n),
reads = sum(n))
Fit <- bakRFit(bakRDataobj, NSS = TRUE)
# Assess raw mutation rates
cB %>% group_by(sample) %>%
summarise(mutrate = sum(TC*n)/sum(nT*n),
reads = sum(n))
plotPCA(Fit$Fast_Fit)
bakR::FnPCA(Fit$Fast_Fit)
bakR::FnPCA(Fit$Fast_Fit, log_kdeg = TRUE)
Fns <- Fit$Fast_Fit$Fn_Estimates
Fns[order(Fns$logit_fn, decreasing = TRUE),]
plot(density(Fns$logit_fn))
plotVolcano(Fit$Fast_Fit)
plot(density(Fit$Fast_Fit$Effects_df$L2FC_kdeg))
Fit_rerun <- bakRFit(Fit, FastRerun = TRUE, StanRateEst = TRUE, NSS = TRUE, RateEst_size = 50, high_reads = 100000000, low_reads = 10)
plotVolcano(Fit_rerun$Fast_Fit)
plot(density(Fit_rerun$Fast_Fit$Fn_Estimates$logit_fn))
plot(density(Fit_rerun$Fast_Fit$Effects_df$L2FC_kdeg))
plotVolcano(Fit_rerun$Fast_Fit)
plot(density(Fit_rerun$Fast_Fit$Fn_Estimates$logit_fn))
plot(density(Fit_rerun$Fast_Fit$Effects_df$L2FC_kdeg))
plot(Fit$Fast_Fit$Fn_Estimates$logit_fn[Fit$Fast_Fit$Fn_Estimates$Exp_ID == 1],
Fit_rerun$Fast_Fit$Fn_Estimates$logit_fn[Fit_rerun$Fast_Fit$Fn_Estimates$Exp_ID == 1])
abline(0,1)
plot(Fit$Fast_Fit$Effects_df$L2FC_kdeg, Fit_rerun$Fast_Fit$Effects_df$L2FC_kdeg)
abline(0,1)
setwd("C:/Users/isaac/Documents/Simon_Lab/Revisions_bakR/Fits/4reps/")
setwd("C:/Users/isaac/Documents/Simon_Lab/Revisions_bakR/Fits/SPT6")
saveRDS(Fit_rerun, file = "MLE_fit_StanRateEst_2022_10_17.rds")
nrow(Fit_rerun$Fast_Fit$Effects_df)
sim <- Simulate_bakRData(5000, fn_mean = -3.5, fn_sd = 0.25)
fit_sim <- bakRFit(sim$bakRData)
# What are raw mutation rates
sim$bakRData$cB %>%
group_by(sample) %>%
summarise(mutrate = sum(TC*n)/sum(nT*n))
sim <- Simulate_bakRData(5000, fn_mean = -4, fn_sd = 0.25)
fit_sim <- bakRFit(sim$bakRData)
# What are raw mutation rates
sim$bakRData$cB %>%
group_by(sample) %>%
summarise(mutrate = sum(TC*n)/sum(nT*n))
fit_sim <- bakRFit(sim$bakRData, StanRateEst = TRUE)
fit_rerun <- bakRFit(fit_sim, FastRerun = TRUE, StanRateEst = TRUE, RateEst_size = 50)
fit_rerun <- bakRFit(fit_sim, FastRerun = TRUE, StanRateEst = TRUE, RateEst_size = 100)
sim <- Simulate_bakRData(5000, fn_mean = -4, fn_sd = 0.25, p_new = 0.1)
fit_sim <- bakRFit(sim$bakRData, StanRateEst = TRUE)
# What are raw mutation rates
sim$bakRData$cB %>%
group_by(sample) %>%
summarise(mutrate = sum(TC*n)/sum(nT*n))
fit_rerun <- bakRFit(fit_sim)
fit_rerun <- bakRFit(fit_sim, FastRerun = TRUE)
Fit <- readRDS("MLE_fit_StanRateEst_2022_10_17.rds")
plotMA(Fit, Model = "MLE")
Counts <- Fit$Data_lists$Count_Matrix
colnames(Counts)
Counts
colnames(Counts)
s4U <- as.factor(rep(rep(c("s4U", "nos4U"), each = 3), times = 2))
s4U
conditions <- data.frame(conditions = conditions,
s4U = s4U)
conditions <- as.factor(rep(c("Auxin", "Ctl"), each = 6))
s4U <- as.factor(rep(rep(c("s4U", "nos4U"), each = 3), times = 2))
conditions <- data.frame(conditions = conditions,
s4U = s4U)
row.names(conditions) <- colnames(Counts)
library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = Counts,
colData = conditions,
design = ~conditions+s4U)
ddso <- DESeq(dds)
plotDispEsts(ddso)
ddso
plotMA(ddso)
results(ddso, contrast = c("conditions", "Auxin", "Ctl"))
res <- results(ddso, contrast = c("conditions", "Auxin", "Ctl"))
plot(density(res$log2FoldChange))
Fit$Fast_Fit$Fn_Estimates
plot(Fit$Fast_Fit$Effects_df$L2FC_kdeg, res$log2FoldChange)
abline(0,1)
plot(Fit$Fast_Fit$Effects_df$effect, res$log2FoldChange)
abline(0,1)
res <- as.data.frame(reso)
reso <- as.data.frame(res)
DE_df <- data.frame(XF = row.names(reso),
log2FoldChange = reso$log2FoldChange,
stat = reso$stat,
padj = reso$padj)
# NSS Heatmap
NSSheatmap <- NSSHeat(Fit, DE_df)
library(pheatmap)
pheatmap(NSSheatmap, cluster_cols = FALSE, show_rownames = FALSE)
sizeFactors(ddso)
# ID some apparent stabilization for further analysis
NSSheatmap
plotVolcano(Fit$Fast_Fit)
view(NSSheatmap)
pheatmap(NSSheatmap, cluster_cols = FALSE, show_rownames = FALSE)
pheatmap(NSSheatmap, cluster_cols = FALSE, show_rownames = FALSE, cluster_rows = FALSE)
?SummarizedExperiment::SummarizedExperiment
library(data.table)
library(tidyverse)
library(bakR)
library(rstan)
library(DEseq2)
library(pheatmap)
theme_mds <-    ggplot2::theme(panel.grid.major = ggplot2::element_blank(),
panel.grid.minor = ggplot2::element_blank(),
panel.background = ggplot2::element_blank(),
axis.line.x = ggplot2::element_line(colour = "black"),
axis.line.y = ggplot2::element_line(colour = "black"),
axis.ticks = ggplot2::element_line(colour = "black"),
axis.text = ggplot2::element_text(color="black", size = 12),
axis.title = ggplot2::element_text(color="black", size = 18),
strip.background = ggplot2::element_blank())
theme_set(theme_mds)
library(DESeq2)
pheatmap(NSSheatmap, cluster_cols = FALSE, show_rownames = FALSE)
Feature_keep <- unique(Fit$Fast_Fit$Fn_Estimates$XF[Fit$Fast_Fit$Fn_Estimates$Exp_ID == 1 & Fit$Fast_Fit$Fn_Estimates$logit_fn > -2,])
Feature_keep <- unique(Fit$Fast_Fit$Fn_Estimates$XF[Fit$Fast_Fit$Fn_Estimates$Exp_ID == 1 & Fit$Fast_Fit$Fn_Estimates$logit_fn > -2])
NSSheatmap_filter <- NSSheatmap[row.names(NSSheatmap) %in% Feature_keep,]
pheatmap(NSSheatmap_filter, cluster_cols = FALSE, show_rownames = FALSE)
view(NSSheatmap_filter)
Feature_keep <- unique(Fit$Fast_Fit$Fn_Estimates$XF[Fit$Fast_Fit$Fn_Estimates$Exp_ID == 1 & Fit$Fast_Fit$Fn_Estimates$logit_fn > -1])
NSSheatmap_filter <- NSSheatmap[row.names(NSSheatmap) %in% Feature_keep,]
pheatmap(NSSheatmap_filter, cluster_cols = FALSE, show_rownames = FALSE)
view(NSSheatmap_filter)
GO_check <- row.names(NSSheatmap_filter[NSSheatmap_filter$Mech_score < -1.2 ])
GO_check <- row.names(NSSheatmap_filter[NSSheatmap_filter$Mech_score < -1.2 ,])
GO_check
library(topGO)
g_list <- goana(GO_check)
library(limma)
g_list <- goana(GO_check)
g_list <- goana(GO_check)
GO_result <- topGO(g_list)
GO_result
g_list
?goana
library(biomaRt)
mart <- biomaRt::useMart(biomart = "ENSEMBL_MART_ENSEMBL",
dataset = "hsapiens_gene_ensembl",
host = "http://www.ensembl.org")
genes <- getBM(filters = "ensembl_gene_id",
attributes = c("ensembl_gene_id","entrezgene_id"),
values = GO_check,
mart = mart)
entrez_list <- genes$entrezgene_id
g_list <- goana(entrez_list,
species = "Hs")
g_list
GO_result <- topGO(g_list)
GO_result
GO_result
# GSEA?
library(fgsea)
GSEA = function(gene_list, GO_file, pval) {
set.seed(54321)
library(dplyr)
library(fgsea)
if ( any( duplicated(names(gene_list)) )  ) {
warning("Duplicates in gene names")
gene_list = gene_list[!duplicated(names(gene_list))]
}
if  ( !all( order(gene_list, decreasing = TRUE) == 1:length(gene_list)) ){
warning("Gene list not sorted")
gene_list = sort(gene_list, decreasing = TRUE)
}
myGO = fgsea::gmtPathways(GO_file)
fgRes <- fgsea::fgsea(pathways = myGO,
stats = gene_list,
minSize=15, ## minimum gene set size
maxSize=400, ## maximum gene set size
nperm=10000) %>%
as.data.frame() %>%
dplyr::filter(padj < !!pval) %>%
arrange(desc(NES))
message(paste("Number of signficant gene sets =", nrow(fgRes)))
message("Collapsing Pathways -----")
concise_pathways = collapsePathways(data.table::as.data.table(fgRes),
pathways = myGO,
stats = gene_list)
fgRes = fgRes[fgRes$pathway %in% concise_pathways$mainPathways, ]
message(paste("Number of gene sets after collapsing =", nrow(fgRes)))
fgRes$Enrichment = ifelse(fgRes$NES > 0, "Up-regulated", "Down-regulated")
filtRes = rbind(head(fgRes, n = 10),
tail(fgRes, n = 10 ))
total_up = sum(fgRes$Enrichment == "Up-regulated")
total_down = sum(fgRes$Enrichment == "Down-regulated")
header = paste0("Top 10 (Total pathways: Up=", total_up,", Down=",    total_down, ")")
colos = setNames(c("firebrick2", "dodgerblue2"),
c("Up-regulated", "Down-regulated"))
g1= ggplot(filtRes, aes(reorder(pathway, NES), NES)) +
geom_point( aes(fill = Enrichment, size = size), shape=21) +
scale_fill_manual(values = colos ) +
scale_size_continuous(range = c(2,10)) +
geom_hline(yintercept = 0) +
coord_flip() +
labs(x="Pathway", y="Normalized Enrichment Score",
title=header) +
th
output = list("Results" = fgRes, "Plot" = g1)
return(output)
}
genes <- row.names(NSSheatmap)
genes <- getBM(filters = "ensembl_gene_id",
attributes = c("ensembl_gene_id","entrezgene_id"),
values = GO_check,
mart = mart)
entrez_all <- genes$entrezgene_id
scores <- NSSheatmap$Mech_score
names(scores) <- entrez_all
GO_file <- "C:/Users/isaac/Documents/Simon_Lab/Revisions_bakR/Data/Gene_sets/c2.all.v2022.1.Hs.entrez.gmt"
GSEA_res <- GSEA(gene_list = scores, GO_file, pval = 0.05)
GSEA_res$Plot
scores
entrez_all
scores
scores <- NSSheatmap$Mech_score
scores
names(scores) <- as.character(entrez_all)
scores
length(entrez_all)
genes <- row.names(NSSheatmap)
genes <- getBM(filters = "ensembl_gene_id",
attributes = c("ensembl_gene_id","entrezgene_id"),
values = genes,
mart = mart)
entrez_all <- genes$entrezgene_id
scores <- NSSheatmap$Mech_score
names(scores) <- as.character(entrez_all)
genes
sum(is.na(genes$entrezgene_id))
genes[is.na(genes$entrezgene_id),]
filter_out <- genes$ensembl_gene_id[is.na(genes$entrezgene_id),]
filter_out <- genes$ensembl_gene_id[is.na(genes$entrezgene_id)]
NSSheatmap_entrez <- NSSheatmap[!(row.names(NSSheatmap) %in% filter_out),]
scores <- NSSheatmap_entrez$Mech_score
genes_filter <_ genes[!(genes$ensembl_gene_id %in% filter_out),]
genes_filter <- genes[!(genes$ensembl_gene_id %in% filter_out),]
names(scores) <- as.character(genes_filter$entrezgene_id)
sum(is.duplicated(genes_filter$entrezgene_id))
sum(duplicated(genes_filter$entrezgene_id))
# Add Entrez ID
NSSheatmap_entrez
# Add Entrez ID
NSSheatmap_entrez$ENSEMBL <- row.names(NSSheatmap_entrez)
NSSheatmap_entrez
NSSheatmap_entrez <- as_tibble(NSSheatmap_entrez)
NSSheatmap_entrez
genes_filter
colnames(genes_filter)
colnames(genes_filter) <- c("ENSEMBL", "entrez")
NSSheatmap_entrez <- left_join(NSSheatmap_entrez, genes_filter, by = "ENSEMBL")
NSSheatmap_entrez
scores <- NSSheatmap_entrez$Mech_score
names(scores) <- NSSheatmap_entrez$entrez
# Perform GSEA
GO_file <- "C:/Users/isaac/Documents/Simon_Lab/Revisions_bakR/Data/Gene_sets/c2.all.v2022.1.Hs.entrez.gmt"
GSEA_res <- GSEA(gene_list = scores, GO_file, pval = 0.05)
scores
sum(is.infinite(scores))
scores
myGO <- fgsea::gmtPathways(GO_file)
fgsea::fgsea(pathways = myGO,
stats = scores,
minSize=15, ## minimum gene set size
maxSize=400, ## maximum gene set size
nperm=10000)
fgRes <- fgsea::fgsea(pathways = myGO,
stats = scores,
minSize=15, ## minimum gene set size
maxSize=400, ## maximum gene set size
nperm=10000)
fgRes <- fgRes %>%
as.data.frame() %>%
dplyr::filter(padj < !!pval) %>%
arrange(desc(NES))
fgRes <- fgRes %>%
as.data.frame() %>%
dplyr::filter(padj < 0.05) %>%
arrange(desc(NES))
nrow(fgRes)
concise_pathways = collapsePathways(data.table::as.data.table(fgRes),
pathways = myGO,
stats = scores)
scores
scores <- scores[order(scores)]
scores
scores <- scores[order(scores, decreasing = TRUE)]
scores
concise_pathways = collapsePathways(data.table::as.data.table(fgRes),
pathways = myGO,
stats = scores)
fgsea:::preparePathwaysAndStats
sum(!is.finite(scores))
any(!is.finite(scores))
all(is.finite(scores))
?any
any(!is.finite(scores))
fgsea::collapsePathways
fgRes
fgRes$pathway
# GSEA?
library("org.Hs.eg.db")
scores
scores <- scores[order(scores, decreasing = TRUE)]
scores
# Perform GSEA
library(clutserProfiler)
# Perform GSEA
library(clusterProfiler)
library(enrichplot)
keytypes(org.Hs.eg.db)
scores <- NSSheatmap$Mech_score
names(scores) <- genes
scores <- scores[order(scores, decreasing = TRUE)]
?gseGO
gse <- gseGO(geneList = scores,
ont = "ALL",
keyType = "ENSEMBL",
nPerm = 1000,
minGSSize = 3,
minGSSize = 800,
pvalueCutoff = 0.05,
verbose = TRUE,
OrgDb = organism)
gse <- gseGO(geneList = scores,
ont = "ALL",
keyType = "ENSEMBL",
nPerm = 1000,
minGSSize = 3,
maxGSSize = 800,
pvalueCutoff = 0.05,
verbose = TRUE,
OrgDb = organism)
scores
scores <- NSSheatmap$Mech_score
genes
genes <- row.names(NSSheatmap)
scores <- NSSheatmap$Mech_score
names(scores) <- genes
scores <- scores[order(scores, decreasing = TRUE)]
scores
gse <- gseGO(geneList = scores,
ont = "ALL",
keyType = "ENSEMBL",
nPerm = 1000,
minGSSize = 3,
maxGSSize = 800,
pvalueCutoff = 0.05,
verbose = TRUE,
OrgDb = organism)
organism
keytypes(organism)
gse <- gseGO(geneList = scores,
ont = "ALL",
keyType = "ENSEMBL",
nPerm = 1000,
minGSSize = 3,
maxGSSize = 800,
pvalueCutoff = 0.05,
verbose = TRUE,
OrgDb = org.Hs.eg.db)
gse <- gseGO(geneList = scores,
ont = "ALL",
keyType = "ENSEMBL",
minGSSize = 3,
maxGSSize = 800,
pvalueCutoff = 0.05,
verbose = TRUE,
OrgDb = org.Hs.eg.db)
library(DOSE)
dotplot(gse, showCategory=10, split=".sign") + facet_grid(.~.sign)
dotplot(gse, showCategory=7, split=".sign") + facet_grid(.~.sign)
NSSheatmap_rescore <- NSSheatmap
NSSheatmap_rescore$XF <- row.names(NSSheatmap_rescore)
NSSheatmap_rescore <- as_tibble(NSSheatmap_rescore)
NSSheatmap_rescore <- NSSheatmap_rescore %>%
mutate(Mech_score = ifelse(XF %in% Feature_keep, Mech_score, ifelse(Mech_score < 0, Mech_score*-1, Mech_score)) )
# Make named list
genes <- row.names(NSSheatmap_rescore)
scores <- NSSheatmap_rescore$Mech_score
names(scores) <- genes
scores <- scores[order(scores, decreasing = TRUE)]
# Perform GSEA
gse <- gseGO(geneList = scores,
ont = "ALL",
keyType = "ENSEMBL",
minGSSize = 3,
maxGSSize = 800,
pvalueCutoff = 0.05,
verbose = TRUE,
OrgDb = org.Hs.eg.db)
scores
genes
# Make named list
genes <- NSSheatmap_rescore$XF
scores <- NSSheatmap_rescore$Mech_score
names(scores) <- genes
scores <- scores[order(scores, decreasing = TRUE)]
# Perform GSEA
gse <- gseGO(geneList = scores,
ont = "ALL",
keyType = "ENSEMBL",
minGSSize = 3,
maxGSSize = 800,
pvalueCutoff = 0.05,
verbose = TRUE,
OrgDb = org.Hs.eg.db)
dotplot(gse, showCategory=7, split=".sign") + facet_grid(.~.sign)
scores
scores[order(scores)]
library(data.table)
library(tidyverse)
library(bakR)
library(rstan)
library(DESeq2)
library(pheatmap)
library(topGO)
library(limma)
library(biomaRt)
theme_mds <-    ggplot2::theme(panel.grid.major = ggplot2::element_blank(),
panel.grid.minor = ggplot2::element_blank(),
panel.background = ggplot2::element_blank(),
axis.line.x = ggplot2::element_line(colour = "black"),
axis.line.y = ggplot2::element_line(colour = "black"),
axis.ticks = ggplot2::element_line(colour = "black"),
axis.text = ggplot2::element_text(color="black", size = 12),
axis.title = ggplot2::element_text(color="black", size = 18),
strip.background = ggplot2::element_blank())
theme_set(theme_mds)
setwd("G:/Shared drives/Matthew_Simon/IWV/bakR_Revisions/Covid_SLAM/")
cB <- as_tibble(fread("cB-221018.csv.gz"))
unique(cB$sample)
cB <- cB %>% filter(!grepl("__no_feature", XF) & !grepl("__ambiguous", XF))
cB <- cB %>% filter(!grepl("__no_feature", XF) | !grepl("__ambiguous", XF))
cB
cB <- cB[, c("sample", "XF", "nT", "TC", "n")]
metadf <- data.frame(tl = c(2, 2, 0, 2, 2, 0),
Exp_ID = c(1, 1, 1, 2, 2, 2))
row.names(metadf) <- unique(cB$sample)
metadf
bakRDataobj <- bakRData(cB, metadf)
# Assess raw mutation rates
cB %>% group_by(sample) %>%
summarise(mutrate = sum(TC*n)/sum(nT*n),
reads = sum(n))
library(data.table)
library(tidyverse)
library(bakR)
library(DESeq2)
library(pheatmap)
?load_all()
?load_all()
getwd()
setwd("C:/Users/isaac/Documents/Simon_Lab/bakR/")
pkgload::load_all()
## Test updates
metadf
getwd()
setwd("C:/Users/isaac/")
pkgload::load_all()
setwd("C:/Users/isaac/Documents/Simon_Lab/HDProbe/")
devtools::load_all()
setwd("C:/Users/isaac/Documents/Simon_Lab/HDProbe/")
devtools::load_all()
setwd("C:/Users/isaac/Documents/Simon_Lab/bakR/")
devtools::install()
# Or test from installation
library(bakR)
# Simulate data
sim <- Simulate_bakRData(100, nreps = 2)
devtools::install()
pkgbuild::compile_dll(force = TRUE)
.Last.error
getwd()
pkgbuild::compile_dll(force = TRUE)
devtools::load_all()
# Or test from installation
library(bakR)
# Simulate data
sim <- Simulate_bakRData(100, nreps = 2)
bakRData <- sim$bakRData
metadf <- bakRData$metadf
cB <- bakRData$cB
metadf$Exp_ID <- ifelse(metadf$Exp_ID == 2, 3, 1)
metadf <- bakRData$metadf
metadf$tl[1] <- 0
metadf <- bakRData$metadf
metadf$Exp_ID <- ifelse(metadf$Exp_ID == 1, 3, 2)
bakRData(cB, metadf)
#pkgbuild::clean_dll()
devtools::load_all()
# Analyze with MLE implementation
Fit <- bakRFit(sim$bakRData, no_ctl = TRUE)
library(rstan)
options(mc.cores = parallel::detectCores())
Fit <- bakRFit(Fit, HybridFit = TRUE, chains = 2, iter = 3000, warmup = 1000)
ls
setwd("C:/Users/isaac/Documents/Simon_Lab/rstanlm/")
pkgbuild::clean_dll()
pkgbuild::compile_dll()
devtools::load_all()
devtools::load_all()