Commit 7700e36a authored by Holger Brandl's avatar Holger Brandl
Browse files

dge workflow improvements

parent 2624510c
......@@ -20,14 +20,15 @@ Options:
opts <- docopt(doc, commandArgs(TRUE)) ## does not work when spining
# opts <- docopt(doc, "--overlay_expr_data ../plot_score_matrix.txt ../degs_by_contrast.txt contrast" )
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.30/R/core_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.30/R/ggplot_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/master/R/bio/cp_utils.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.36/R/core_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.36/R/ggplot_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.36/R/bio/cp_utils.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.36/R/bio/diffex_commons.R")
load_pack(knitr)
load_pack(DT)
#loadpack(clusterProfiler)
#load_pack(clusterProfiler)
#devtools::session_info() # nice!
......@@ -69,7 +70,7 @@ qCutoff <- as.numeric(opts$qcutoff)
#' ## Enrichment Analysis
#' Run configuration was
vec2df(unlist(opts)) %>% filter(!str_detect(name, "^[<-]")) %>% kable()
vec_as_df(unlist(opts)) %>% filter(!str_detect(name, "^[<-]")) %>% kable()
#' This analysis was performed using [clusterProfiler](http://bioconductor.org/packages/release/bioc/html/clusterProfiler.html). The following ontologies were tested: Kegg, Go, Reactome, Dose,
......@@ -92,7 +93,7 @@ annoDb <- guess_anno_db(geneLists$ensembl_gene_id) # e.g. "org.Hs.eg.db"
#' ### Convert to entrez gene ids
install_package(clusterProfiler)
install_package("clusterProfiler")
glMappedRaw <- clusterProfiler::bitr(unique(geneLists$ensembl_gene_id), fromType="ENSEMBL", toType="ENTREZID", OrgDb=annoDb) %>%
set_names("ensembl_gene_id", "entrez_gene_id") %>%
......@@ -128,7 +129,9 @@ glMapped %<>% filter(!is.na(entrez_gene_id))
#}
## prety slow
enrResults <- quote(glMapped %>% do(cp_test(.$entrez_gene_id))) %>% cache_it(paste0("enrdata_", digest(glMapped)))
## user purrr::map here
enrResults <- quote(glMapped %>% do(cp_test(.$entrez_gene_id, annoDb, cpSpecies, q_cutoff=qCutoff))) %>%
cache_it(paste0("enrdata_", digest(glMapped)))
## run the actual enrichment test for all clusters and all ontologies
#library(foreach); library(doMC); registerDoMC(cores=20)
......@@ -151,7 +154,7 @@ write_tsv(enrResults, path=paste0(resultsBaseName, "enrResults.txt"))
# enrResults <- read.delim(paste0(resultsBaseName, "enrResults.txt"))
#' [Enrichment Results](`r paste0(resultsBaseName, "enrResults.txt")`)
loadpack(DT)
load_pack(DT)
datatable(enrResults)
#enrResults %>% ggplot(aes(pvalue)) + geom_histogram() + scale_x_log10()
......@@ -228,6 +231,21 @@ warning("dropping levels")
erPlotData %<>% mutate(ontology=ac(ontology)) ## drop unsused level to get consistent color palette
erPlotData %<>% rename(Term=Description)
## TODO remove with v1.37
create_palette <- function(x, pal = 'Set1'){
load_pack(RColorBrewer)
ux <- sort(unique(x))
n <-length(ux)
if(n==0) return(c())
setNames(brewer.pal(name = pal, n = n)[1:n], ux)
}
term_category_colors <- create_palette(unique(ac(erPlotData$ontology)))
figDir <- "enr_charts"
......@@ -287,7 +305,7 @@ do({
})
#+ results="asis"
l_ply(term_barplot_files$file, function(pngFile){ cat(paste0("<img src='", pngFile, "'><br>"))})
walk(term_barplot_files$file, function(pngFile){ cat(paste0("<img src='", pngFile, "'><br>"))})
########################################################################################################################
#' ## Enriched KEGG Pathways
......@@ -299,12 +317,12 @@ l_ply(term_barplot_files$file, function(pngFile){ cat(paste0("<img src='", pngFi
keggPathways <- enrResults %>% filter(ontology=="kegg") %$% ac(ID) %>% unique() # %>% head(3) ## todo remove debug head
##+ results='asis'
#if(!exists("keggPathways") | nrow(keggPathways)==0){
# cat("No enriched pathways found")
#}
# results='asis'
if(!exists("keggPathways") || nrow(keggPathways)==0){
cat("No enriched pathways found")
}else{
loadpack(naturalsort)
load_pack(naturalsort)
#reload_dplyr()
......@@ -354,12 +372,12 @@ unload_packages()
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.26/R/core_commons.R")
#loadpack(plyr)
#loadpack(dplyr)
#loadpack(dplyr)
#loadpack(stringr)
loadpack(pathview)
loadpack(png)
#load_pack(plyr)
#load_pack(dplyr)
#load_pack(dplyr)
#load_pack(stringr)
load_pack(pathview)
load_pack(png)
......@@ -450,7 +468,7 @@ makeTooltip <- function(entrez_id){
#+ results="asis", echo=FALSE
## #+ results="asis", echo=FALSE
## simple non-clickable plots
## http://stackoverflow.com/questions/12588323/r-how-to-extract-values-for-the-same-named-elements-across-multiple-objects-of
......@@ -507,7 +525,7 @@ createImgMap <- function(plotData){
#keggNodes %>% rowwise() %>% {curNode=.; cat(curNode$name)}
#see http://www.html-world.de/180/image-map/
keggNodes %>% a_ply(1, function(curNode){
keggNodes %>% plyr::a_ply(1, function(curNode){
rectDef=with(curNode, paste(x, y, x+width, y+height, sep=","))
# paste0("<area href='", curNode$link, "' alt='Text' coords='", rectDef , "' shape='rect'>") %>% cat()
paste0("<area href='", curNode$link, "' title='", curNode$tooltip, "' alt='Text' coords='",rectDef , "' shape='rect'>") %>% cat()
......@@ -516,11 +534,12 @@ createImgMap <- function(plotData){
cat("</map></p><br>")
}
pathwayPlots %>% l_ply(createImgMap)
pathwayPlots %>% plyr::l_ply(createImgMap)
## respin it for cild inclusion
# require(knitr); setwd("/Volumes/projects/bioinfo/scripts/ngs_tools/dev/common"); spin("cp_enrichment.R", knit=F)
}
#' ## Notes
#' Other interesting tools to perform enrichment testing
......
......@@ -70,6 +70,10 @@ sampleSheet %>% group_by(bio_sample) %>% summarise(
) %$%
zcat %>%
write_lines("lane_merge.cmd")
sampleSheet %>% distinct(replicate, sample) %>%
arrange(sample, replicate) %>%
write_tsv(path="basic_design.txt")
' | R --vanilla -q
cat lane_merge.cmd | while read line; do
......@@ -100,7 +104,17 @@ mailme "$project: mapping done"
mcdir $baseDir/dge_analysis
rend.R -e ${NGS_TOOLS}/dge_workflow/featcounts_deseq.R ../alignments/star_counts_matrix.txt 2>&1 | tee featcounts_deseq.log
## build custom design matrix
#csvcut -tc replicate,sample $baseDir/lanereps_pooled/renaming_scheme.txt | csvformat -T > design_matrix.txt
#Rscript - <<"EOF"
#devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.36/R/core_commons.R")
#file.path(Sys.getenv("baseDir"), "lanereps_pooled/renaming_scheme.txt") %>%
# read_tsv() %>%
# distinct(replicate, sample) %>%
# arrange(sample, replicate) %>%
# write_tsv("design_matrix.txt")
#EOF
## or use just some contrasts of interest
......@@ -109,12 +123,17 @@ rend.R -e ${NGS_TOOLS}/dge_workflow/featcounts_deseq.R ../alignments/star_counts
#unpolarised,liver_polar_stage3
#" | csvformat -T > contrasts.txt
#
#rend.R -e $NGS_TOOLS/dge_workflow/featcounts_deseq.R --contrasts contrasts.txt ../alignments/star_counts_matrix.txt 2>&1 | tee featcounts_deseq.log
rend.R -e ${NGS_TOOLS}/dge_workflow/featcounts_deseq_mf.R ../alignments/star_counts_matrix.txt 2>&1 | tee featcounts_deseq.log
#rend.R -e ${NGS_TOOLS}/dge_workflow/featcounts_deseq_mf.R --design 'batch+sample' --pcutoff 0.01 --contrasts ../contrasts_wt.txt ../qseq.counts.txt ../qseq.design_matrix.txt 2>&1 | tee featcounts_deseq_mf.log
## Term enrichment
mcdir ${baseDir}/dge_analysis/dge_enrichment_analysis
$NGS_TOOLS/common/cp_enrichment.R --overlay_expr_data ../plot_score_matrix.txt ../degs_by_contrast.txt contrast
#$NGS_TOOLS/common/cp_enrichment.R --overlay_expr_data ../plot_score_matrix.txt ../degs_by_contrast.txt contrast
rend.R -e ${NGS_TOOLS}/common/cp_enrichment.R --overlay_expr_data ../plot_score_matrix.txt ../degs_by_contrast.txt contrast
......@@ -123,7 +142,7 @@ rend.R -e ${NGS_TOOLS}/common/cp_enrichment.R --overlay_expr_data ../plot_score_
## bidirectional sync with project space
## todo define mount path on bioinfo for bidirectional synching
## ~/bin/unison $baseDir ssh://bioinfo///home/brandl/mnt/<<MOUNT_PATH>> -fastcheck true -times -perms 0
# ~/bin/unison ${baseDir:=error} /net/fileserver-nfs/stornext/snfs3/projects/huttnerlab/Paul/fam72d_transfection_seq -fastcheck -times -perms 0 -batch
# or use a uni-directional sync
#rsync -avsn --delete ${baseDir} brandl@fileserver:/projects//file/server/path
......
......@@ -9,12 +9,15 @@ else
fi
export PATH=${BIO_BIN_BASE}/bowtie2-2.2.5:$PATH
export PATH=${BIO_BIN_BASE}/deepTools-2.2.2/bin:$PATH
export PATH=${BIO_BIN_BASE}/FastQC_0.11.2:$PATH
export PATH=${BIO_BIN_BASE}/bedtools2-2.25.0/bin/:$PATH
export PATH=${BIO_BIN_BASE}/samtools-1.3:$PATH
export PATH=${BIO_BIN_BASE}/STAR-2.5.2b/source:$PATH
export PATH=${BIO_BIN_BASE}/ucsc:$PATH # for wigToBigWig
## Fixme use BIO_BIN_BASE for deeptools
#export PATH=${BIO_BIN_BASE}/deepTools-2.2.2/bin:$PATH
export PATH=/home/brandl/.local/bin/:$PATH
#export PATH=/home/brandl/bin/subread-1.4.6-p3-Linux-x86_64/bin:$PATH
......
......@@ -28,9 +28,9 @@ library(DESeq2)
#load_pack(readr)
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.35/R/core_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.35/R/ggplot_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.35/R/bio/diffex_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.36/R/core_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.36/R/ggplot_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.36/R/bio/diffex_commons.R")
## process input arguments
......@@ -54,7 +54,7 @@ assert(designFormula=="sample" || !is.null(design_matrix_file), "custom designs
#' # Differential Expression Analysis
#' Run configuration was
vec2df(unlist(opts)) %>% filter(!str_detect(name, "^[<-]")) %>% kable()
vec_as_df(unlist(opts)) %>% filter(!str_detect(name, "^[<-]")) %>% kable()
########################################################################################################################
#' ## Data Preparation
......@@ -95,6 +95,7 @@ genesAfter = nrow(countMatrix)
#' To understand fold-change shrinking and estimation check out
#' [Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2](http://genomebiology.com/2014/15/12/550/abstract)
## todo this is required now , so why all the checking
if(!is.null(design_matrix_file)){
replicates2samples = read_tsv(design_matrix_file) #%>% select(replicate, sample)
}else{
......@@ -107,7 +108,7 @@ if(!is.null(design_matrix_file)){
if(!is.null(contrasts_file)){
contrasts = read_tsv(contrasts_file)
}else{
stopifnot(is.null(design_matrix_file)) ## cant work yet because sample column with have random name
assert(!is.null(replicates2samples), "Can not auto-create contrasts without design matrix") ## cant work yet because sample column with have random name
contrasts = distinct_all(replicates2samples, sample) %>% select(sample) %>%
merge(.,., suffixes=c("_1", "_2"), by=NULL) %>%
filter(ac(sample_1)>ac(sample_2)) %>%
......@@ -156,7 +157,7 @@ exDesign = data_frame(replicate=colnames(countMatrix)) %>% mutate(col_index=row_
#exDesign = replicates2samples %>% arrange(colnames(countMatrix))
# valiadate that contrasts model is compatible with data
designSamples = as.df(exDesign)[, contrastAttribute] %>% unique %>% sort
designSamples = as_df(exDesign)[, contrastAttribute] %>% unique %>% sort
contrastSamples = contrasts %>% gather %$% value %>% unique %>% sort
if(!all(contrastSamples %in% designSamples)){
warning("design : ", paste(designSamples, collapse=", "))
......
......@@ -20,9 +20,9 @@ python2 setup.py install --user
cd ~/bin
wget -O deepTools-2.3.5.tar.gz wget https://github.com/fidelram/deepTools/archive/2.3.5.tar.gz
tar xvf deepTools-2.3.5.tar.gz
cd deepTools-2.3.5
wget -O deepTools-2.4.2.tar.gz wget https://github.com/fidelram/deepTools/archive/2.4.2.tar.gz
tar xvf deepTools-2.4.2.tar.gz
cd deepTools-2.4.2
#http://effbot.org/pyfaq/when-importing-module-x-why-do-i-get-undefined-symbol-pyunicodeucs2.htm
# https://github.com/galaxyproject/tools-iuc/issues/256
......@@ -33,6 +33,7 @@ python3 setup.py install --user
#which /home/brandl/.local/bin/correctGCBias
#which /home/brandl/.local/bin/bamCoverage
/home/brandl/.local/bin/bamCoverage
/home/brandl/.local/bin/multiBamSummary
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment