Delete SpicyBytesStatsCold.R~

SpicyBytesStatsCold.R~

-args <- commandArgs(trailingOnly = TRUE)
-treatment = args[1]
-mock = args[2]
-
-# Adds path to library path on leunissen server
-.libPaths(c('/prog/BIF30806/project/SpicyBytes/R-libraries',.libPaths()))
-# Install DESeq2 and sleuth packages if needed
- if (!requireNamespace("BiocManager", quietly = TRUE))
-   install.packages("BiocManager")
-
-#BiocManager::install("DESeq2")
-if (!requireNamespace("remotes", quietly = TRUE))
-  BiocManager::install("remotes")
-if (!requireNamespace("pachterlab/sleuth", quietly = TRUE))
-  BiocManager::install("pachterlab/sleuth")
-
-## for using sleuth; https://hbctraining.github.io/DGE_workshop_salmon/lessons/09_sleuth.html
-
-#load sleuth
-library(sleuth)
-#set working directory, this is where each kalisto output is stored in its own directory
-setwd("~/prog/BIF30806/project/SpyciBytes/Kallisto_result") #check if this works
-
-
-base_dir = "."
-sample_ids = dir(base_dir,"SRR")
-kal_dirs = sapply(sample_ids, function(id) file.path(base_dir, id))
-
-## design.txt has to be created, should contain 2 columns sample	and condition
-s2c = read.table(file.path(base_dir,"../design1.txt"), header=TRUE, stringsAsFactors=FALSE)
-# condition should be all the conditions, first per sample, then all the names
-#condition = factor(c("ST","ST","ST","HT","HT","HT"),c("ST","HT"))
-condition <- factor(s2c$condition, levels = unique(s2c$condition))
-
-s2c$path = kal_dirs
-
-
-#print(s2c)
-
-so = sleuth_prep(s2c, ~conditionCold,extra_bootstrap_summary=TRUE, transformation_function = function(x) log2(x + 0.5))
-
-# compares two different models
-so = sleuth_fit(so, ~conditionCold, 'full')
-so = sleuth_fit(so, ~1, 'reduced')
-#so = sleuth_lrt(so, 'reduced', 'full')
-so = sleuth_wt(so, which_beta='conditionAbio_mock', which_model='full')
-sleuth_table = sleuth_results(so, 'conditionAbio_mock', 'wt', show_all=FALSE) 
-
-#find all significant differentially expressed genes with a = 0.05
-sleuth_significant <- sleuth_table[sleuth_table$qval<=0.05,]
-write.csv(sleuth_significant, file="../Stat_results/significant_genes_cold.csv", row.names=FALSE)
-
-#head(sleuth_significant)
-
-
-# plot_bootstrap(so, "AT1G56600.1")
-# plot_bootstrap(so,"AT1G01810.1")
-
-
-# so = sleuth_wt(so,which_beta="conditionST", which_model="full")
-
-# ## ------------------------------------------------------------------------
-# sleuth_table <- sleuth_results(so, 'conditionST', 'wt', show_all=FALSE)
-
-
-# ## ------------------------------------------------------------------------
-# sleuth_matrix = sleuth_to_matrix(so, 'obs_norm', 'tpm')
-# transcripts = rownames(sleuth_matrix)
-# t2g = data.frame(target_id=transcripts, gene_id=substring(transcripts,1,9),stringsAsFactors=FALSE)
-
-# ## ------------------------------------------------------------------------
-# so = sleuth_prep(s2c, ~ condition,target_mapping = t2g, aggregation_column = 'gene_id')
-# so = sleuth_fit(so, ~condition, 'full')
-# so = sleuth_fit(so, ~1, 'reduced')
-# so = sleuth_wt(so,which_beta="conditionST", which_model="full")
-# sleuth_table <- sleuth_results(so, 'conditionST', 'wt', show_all=FALSE)
-# sleuth_significant <- sleuth_table[sleuth_table$qval<=0.05,]
-
-
-#nrow(sleuth_significant)
-