igv_track_range.sh

##todo merge with create_igv_session.R

docopts...

opts$bamfile
...


bamFile=/home/brandl/mnt/chip-seq_study/ChIPSeq_March_2015/data/alignments_mmfilt/K4M_Oblong_mmf.bam
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("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()+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