Commit c074b546 authored by Holger Brandl's avatar Holger Brandl

cont. igv session creator

parent d297c25f
......@@ -22,7 +22,7 @@ Bioinfo
To use the structure from above when working on bioinformatics-srv1 just use
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/home/brandl/mnt/mack/bioinfo/scripts/ngs\_tools/v1.0
/home/brandl/mnt/mack/bioinfo/scripts/ngs_tools/v1.0
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
ToDo
......
#!/usr/bin/env Rscript
#' # Create an IGV Session File
#' ## Create an IGV Session File
#' hallo
suppressMessages(library(docopt))
doc <- 'Usage: igvtrack.R --genome [options] <genome> <bam_files>...
Options:
--cov_quantils <cov_quantils> Use a q-value cutoff of 0.01 instead of a q-value cutoff [default: 0.8]
'
opts <- docopt(doc)
#opts <- docopt(doc, "--cov_quantils 0.7 --genome musmus test.bam naha.bam", quoted_args=T)
print(opts)
'
# https://github.com/holgerbrandl/datautils/blob/master/R/spinr/spin.R
export PATH=/projects/bioinfo/scripts/ngs_tools/dev/chipseq_workflow:$PATH
create_igv_session.R --genome asdf lala.bam
spin.R $(which create_igv_session.R) \""--genome asdf lala.bam"\"
'
require(knitr)
require(ggplot2)
require(dplyr)
iris %>% head(5) %>% kable()
iris %>% ggplot(aes(Sepal.Width)) + geom_histogram()
#require(whisker)
#bamFile="surprise.bam"
#whisker.render("<igg file={{bamFile}}</igv>") %>% writeLines()
#paste0("<igg file=", bamFile, "</igv>")
\ No newline at end of file
##todo merge with create_igv_session.R
docopts...
opts$bamfile
......@@ -9,12 +11,20 @@ bamBaseName=$(basename $bamFile .bam)
genomeCoverageBed -d -ibam $bamFile | cut -f3 | sort | uniq -c | trim > ${bamBaseName}.coveragedist.txt
echo '
setwd("/Volumes/albert_chipseq/Naha_analysis/BigWig_Files/Combined_BAMS_20141024")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.9/R/core_commons.R")
require(ggplot2)
covDist = read.delim("K4M_Oblong_mmf.coveragedist.txt", header=F, sep=" ") %>% rename(num_bases=V1, coverage=V2)
covDist = read.delim("NEC_K27_Combined.remdup.bam.sorted.coveragedist.txt", header=F, sep="")
names(covDist)=c("num_bases","coverage")
covDist %<>% arrange(coverage) %>% mutate(total_cov=sum(coverage), cum_coverage=cumsum(coverage), cum_cov_prop=cum_coverage/total_cov)
covDist %>% ggplot(aes(coverage, cum_cov_prop)) + geom_area()
' | R --vanilla
\ No newline at end of file
covDist %>% ggplot(aes(coverage, cum_cov_prop)) + geom_area()+geom_hline(y=0.80, color="red")
cutoff=0.8
igv_track_range<- subset(covDist,cum_cov_prop >= cutoff)[1,]$coverage
print(igv_track_range)
' | R --vanilla
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