Commit f3d79628 authored by Holger Brandl's avatar Holger Brandl

cosmetics

parent 7973fe51
......@@ -479,23 +479,13 @@ deResults %>%
facet_grid(sample_1 ~ sample_2)
# Define absolute binding categories
#rawCounts = counts(dds,normalized=F) %>%
# set_names(exDesign(dds)$condition) %>% rownames2column("ensembl_gene_id")
#ggplot(rawCounts, aes(H3HA_Sphere)) + geom_histogram() + scale_x_log10() + ggtitle("raw H3HA_Sphere counts distribution")
#require_auto(Hmisc)
#
#bindCats = rawCounts %>% transmute(ensembl_gene_id, bind_category=cut2(H3HA_Sphere, cuts=c(10, 100)))
#with(bindCats, as.data.frame(table(bind_category))) %>% kable()
########################################################################################################################
#' ## Count Table Exports
## optionally install bioconducor package
install_package("biomaRt")
#if(!any(.packages(all.available=TRUE)=="biomaRt")){
# source("http://bioconductor.org/biocLite.R")
# biocLite("biomaRt", ask=FALSE)
#}
# Load transcriptome annotations needed for results annotation
geneInfo = quote({
......@@ -519,7 +509,7 @@ write_tsv(geneInfo, path = paste0(resultsBase, "geneInfo.txt"))
#' Export counts per replicate as double-normalized FPKMs=(read_count * 10^9)/(gene_length * library_size)
fpkms = countMatrix %>%
as.data.frame %>%
add_rownames("ensembl_gene_id") %>%
rownames_to_column("ensembl_gene_id") %>%
gather(replicate, num_tags, - ensembl_gene_id) %>%
tbl_df %>%
left_join(transmute(geneInfo, ensembl_gene_id, gene_length = end_position - start_position), by = "ensembl_gene_id") %>%
......@@ -542,7 +532,7 @@ fpkms %>%
sampleFpkms = countMatrix %>%
as_df() %>%
add_rownames("ensembl_gene_id") %>%
rownames_to_column("ensembl_gene_id") %>%
gather(replicate, num_tags, - ensembl_gene_id) %>%
tbl_df %>%
## grouping if sample column is known to be called sample
......@@ -571,12 +561,12 @@ sampleFpkms %>%
# Export the complete dataset for later analysis
deAnnot = deResults %>%
deAnnot = deResults %>% left_join(geneInfo)
#inner_join(normCounts) %>%
#merge(., contrasts) %>%
#merge(.,., suffixes=c("_1", "_2"), by=NULL) %>%
#inner_join(baseMeanPerLvl, id='ensembl_gene_id') %>%
left_join(geneInfo)
write_tsv(deAnnot, path = paste0(resultsBase, "de_results.txt"))
#deAnnot = read_tsv(paste0(resultsBase, "de_results.txt"))
......@@ -719,8 +709,9 @@ kableDegs %>%
########################################################################################################################
#' ## MA Plots
#' Redo MA plots but now including hit overlays
## todo round of MA should be enough.
#' Redo MA plots but now including hit overlays
maPlots = deResults %>%
group_by(sample_1, sample_2) %>%
do(gg = { maData <- .
......@@ -733,46 +724,3 @@ maPlots = deResults %>%
#+ fig.height=6*ceiling(length(maPlots)/2)
multiplot(plotlist = maPlots, cols = min(2, length(maPlots)))
#'
##--------------------------------------------------------
## Useful link: https://gist.github.com/stephenturner/f60c1934405c127f09a6
## TODO later, reenable child-inclusion of enrichment analysis
########################################################################################################################
## ## Term enrichment Data Preparation
##+ echo=FALSE
#
## This analysis was performed using [David](http://david.abcc.ncifcrf.gov/). The following ontologies were tested: `r paste(ontologies, collapse=', ')`
#
#geneLists = degs %>%
# transmute(ensembl_gene_id, list_id=paste(sample_1, "vs", sample_2))
#
#if(nrow(contrasts)<4){
# geneLists = bind_rows(
# geneLists,
# degs %>% filter(s1_overex) %>% transmute(ensembl_gene_id, list_id=paste(sample_1, ">", sample_2)),
# degs %>% filter(!s1_overex) %>% transmute(ensembl_gene_id, list_id=paste(sample_1, "<", sample_2))
# )
#}
#
### additional overlaps before doing the intersection analysis
#geneLists %>% count(list_id) %>% kable()
#
#intersectLists = function(geneLists, listIdA, listIdB, intersectListId) {
# commonGenes = setdiff(filter(geneLists, list_id==listIdA)$ensembl_gene_id, filter(geneLists, list_id==listIdB)$ensembl_gene_id)
# data.frame(list_id=intersectListId, ensembl_gene_id=commonGenes)
#}
#
#geneLists %<>% group_by(list_id)
#save(geneLists, file=".enrGeneLists.RData")
## geneLists = local(get(load("enrGeneLists.RData")))
#
#
### redefine opts arguments and tweak knitr options
#opts_knit$set(root.dir = getwd())
#commandArgs = function(x) c(paste("--overlay_expr_data ", count_matrix_file, " enrGeneLists.RData"))
##source("/home/brandl/mnt/mack/bioinfo/scripts/ngs_tools/dev/common/david_enrichment.R")
## child='/home/brandl/mnt/mack/bioinfo/scripts/ngs_tools/dev/common/david_enrichment.Rmd'
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