# 18 Dec 2010, 17h00
# Code from Mark Robinson, 13, 16 Dec 2010
# source("/projects/remc_bigdata/Karsan/R/20101219/edgeR_TSS-rgns_DE-H4ac.20101219.R")

#install.packages("/home/grobertson/linux_x86_64/R/Rsamtools_1.2.2.tar.gz",type="source")
#install.packages("/home/grobertson/linux_x86_64/R/IRanges_1.8.7.tar.gz",type="source")


##############################################################################
options(stringsAsFactors=FALSE)

library(rtracklayer)
session <- browserSession("UCSC")
genome(session) <- "hg18"
#trackNames(session) ## list the track names
q <- ucscTableQuery(session, "knownGene")
cat("getting UCSC knownGene table",date(),"\n")
knownGene <- getTable(q)
cat("done",date(),"\n")


# merge TSSs with identical starts
knownGene$position <- ifelse(knownGene$strand=="+", knownGene$txStart, knownGene$txEnd)
key <- paste(knownGene$chrom, knownGene$position, sep=":")
ukey <- unique(key)
m <- match(ukey, key)
mkg <- knownGene[m,c("chrom","strand","txStart","txEnd","name","position")]  # mkg = merged known genes

# create merged IDs
rownames(mkg) <- ukey
ids <- split(knownGene$name, key)
collapsedids <- sapply(ids, paste, collapse=",")
mkg$newID <- ""
mkg[names(collapsedids),"newID"] <- collapsedids

# merge nearby TSSs 
tssTol <- 500  # if distance b/w TSSs is less than this, they are merged
splitTol <- 250000

lockey <- paste(mkg$chrom, mkg$strand, sep=":")
ind <- seq_len(nrow(mkg))
inds <- split(ind, lockey)

clustno <- lapply(inds, FUN=function(u) {
  w <- diff(mkg$position[u]) > splitTol
  z <- cumsum(c(0,w))
  s <- split(u, z)
  cat(".")
  clusts <- lapply(s, FUN=function(uu) {
     if(length(uu)==1) return(1)
     d <- dist(mkg$position[uu])
     h <- hclust(d,"ave")
     cutree(h,h=tssTol)
  })
  paste( rep(names(s), sapply(s,length)), unlist(clusts, use.names=FALSE), sep=".")
})

clustno <- unsplit( clustno, lockey )
clustkey <- paste( mkg$chrom, mkg$strand, clustno, sep=":" )

uckey <- unique(clustkey)
mc <- match(uckey, clustkey)
mkg1 <- mkg[mc,]

rownames(mkg1) <- uckey
ids <- split(mkg$newID, clustkey)
collapsedids <- sapply(ids, paste, collapse=";")
mkg1$newID <- ""
mkg1[names(collapsedids),"newID"] <- collapsedids

rownames(mkg1) <- NULL

anno <- data.frame(chr=mkg1$chrom, strand=mkg1$strand, start=mkg1$txStart, end=mkg1$txEnd, name=mkg1$name, allIDs=mkg1$newID)

cat("\nTSS's are merged",date(),"\n")


##############################################################################
cat("read in BAM files\n")
library(Rsamtools)
library(Repitools)
library(edgeR)
library(BSgenome.Hsapiens.UCSC.hg18)

cat("set scanBam parameters\n")
# need to add a failed-chastity filter and a MAPQ filter
#p <- ScanBamParam(what=c("rname", "strand", "pos"), flag=scanBamFlag(isUnmappedQuery=FALSE,isDuplicate=FALSE))
# 16 Dec 2010, filter by MAPQ
p <- ScanBamParam(what=c("rname", "strand", "pos", "mapq"), flag=scanBamFlag(isUnmappedQuery=FALSE,isDuplicate=FALSE))

chrNames <- paste("chr", c(1:22, "X", "Y"), sep = "")

fragSize <- 200

##############################################################################
# 16 Dec 2010
#readLen <- c(50,75,50,50)  # vector to match filenames
#filenames <- c("/projects/remc_bigdata/Karsan/HS1238_kd/maq2sam/HS1238.h.sorted.bam",
#               "/projects/remc_bigdata/Karsan/HS1235_ctl/maq2sam/HS1235.h.sorted.bam",
#               "/projects/analysis/analysis5/HS1238/62P7NAAXX_2/bwa/62P7NAAXX_2.bam",
#               "/projects/analysis/analysis5/HS1235/62P7NAAXX_1/bwa/62P7NAAXX_1.bam")
##############################################################################
# 18 Dec 2010 - just 'lane-2' data
readLen <- c(50,50)  # vector to match filenames
filenames <- c("/projects/analysis/analysis5/HS1238/62P7NAAXX_2/bwa/62P7NAAXX_2.bam",
               "/projects/analysis/analysis5/HS1235/62P7NAAXX_1/bwa/62P7NAAXX_1.bam")

cat("now read in .bam files\n")
##############################################################################
# 16 Dec 2010: old way without filtering by MAPQ
#gr <- mapply(FUN=function(u,v) {
#  sb <- scanBam(u, param=p)[[1]]
#  GRanges(seqnames=paste("chr",sb$rname,sep=""), ranges=IRanges(start=sb$pos,width=readLen), strand=sb$strand)
#  },filenames,readLen)
#names(gr) <- c("HS1238_kd","HS1235_ctl","HS1238_kd_2","HS1235_ctl_2")
#names(gr) <- c("HS1238_kd_2","HS1235_ctl_2")

##############################################################################
# 16 Dec, Mark Robinson: add MAPQ filter here
# 18 Dec 2010 - use for lane-2 data
gr <- mapply(FUN=function(u,v) {
  sb <- scanBam(u, param=p)[[1]]
  w <- sb$mapq > 10  # SET CUTOFF HERE
  cat(u, ": Read", length(sb$rname), "mapped non-dup. reads,", round(mean(w)*100,2), "percent with MAPQ>10\n")
  GRanges(seqnames=paste("chr",sb$rname[w],sep=""), ranges=IRanges(start=sb$pos[w],width=v), strand=sb$strand[w])
  },filenames,readLen)
#names(gr) <- c("HS1238_kd","HS1235_ctl","HS1238_kd_2","HS1235_ctl_2")
names(gr) <- c("HS1238_kd_2","HS1235_ctl_2")
##############################################################################
grl <- GRangesList(gr)


cat("TSS regions: get counts\n")
#cat("bpUp=1500, bpDown=1000", "\n")
#counts <- annotationCounts(grl, anno, bpUp=1500, bpDown=1000, seqLen=fragSize, verbose=TRUE)
cat("bpUp=2000, bpDown=1000", "\n")
counts <- annotationCounts(grl, anno, bpUp=2000, bpDown=1000, seqLen=fragSize, verbose=TRUE)
#cat("bpUp=2500, bpDown=1000", "\n")
#counts <- annotationCounts(grl, anno, bpUp=2500, bpDown=1000, seqLen=fragSize, verbose=TRUE)
#cat("bpUp=2500, bpDown=1500", "\n")
#counts <- annotationCounts(grl, anno, bpUp=2500, bpDown=1500, seqLen=fragSize, verbose=TRUE)

colSums(counts)
#HS1238_kd HS1235_ctl 
#  3070716    3288149
k <- rowSums(counts) > 10 
f <- calcNormFactors(counts[k,])


##############################################################################
#print("estimate common dispersion as 'dhack'")
#dhack <- d <- DGEList(counts=counts[k,], 
#	         group=colnames(counts), 
#             lib.size=colSums(counts[k,])*f,
#             genes=anno[k,1:3])
#             #genes=as.data.frame(windows)[k,1:3])
#dhack$samples$group <- factor(c("A","A"))
#dhack <- estimateCommonDisp(dhack)
#d <- estimateCommonDisp(d)  # in my example, no bio replicates, will get warning
#d$common.dispersion <- dhack$common.dispersion
#cat("TSSrgns: dhack$common.dispersion=",dhack$common.dispersion,"\n")


##############################################################################
# Old way, without 'dhack', for a single run pair
# 18 Dec 2010 - resusing this for lane-2 data
d <- DGEList(counts=counts[k,], 
            group=colnames(counts), 
            lib.size=colSums(counts[k,])*f,
            genes=anno[k,1:3])
d <- estimateCommonDisp(d)
cat("d$common.dispersion=",d$common.dispersion,"\n")

#print("exactTest()")
#de <- exactTest(d,pair=c("HS1235_ctl_2","HS1238_kd_2"))
#topTags(de)

#xx <- cbind(anno[k,-5], counts[k,], id=rownames(de$table), de$table, adjp=p.adjust(de$table$p.value, method="BH") )

#print("write.table, all records, regardless of p-val")
#write.table(xx, "/projects/remc_bigdata/Karsan/R/20101218/HS1238_vs_1235.lane-2.mapq10-noDups.knownGene.merge-500.2p5kb-TSS-1p5kb.mapq10_noDups.edgeR.all.20101218.txt", sep="\t", row.names=FALSE, quote=FALSE)
#write.table(xx, "/projects/remc_bigdata/Karsan/R/20101218/HS1238_vs_1235.lane-2.mapq10-noDups.knownGene.merge-500.2p5kb-TSS-1kb.mapq10_noDups.edgeR.all.20101218.txt", sep="\t", row.names=FALSE, quote=FALSE)
#write.table(xx, "/projects/remc_bigdata/Karsan/R/20101218/HS1238_vs_1235.lane-2.mapq10-noDups.knownGene.merge-500.1p5kb-TSS-1kb.mapq10_noDups.edgeR.all.20101218.txt", sep="\t", row.names=FALSE, quote=FALSE)
#write.table(xx, "/projects/remc_bigdata/Karsan/R/20101218/HS1238_vs_1235.lane-2.mapq10-noDups.knownGene.merge-500.2kb-TSS-1kb.mapq10_noDups.edgeR.all.20101218.txt", sep="\t", row.names=FALSE, quote=FALSE)

#print(date())

#cat("save.image()\n")

#save.image("Robjects.2p5kb-TSS-1p5kb.2010-12-18-Vancouver.Rdata") # save all objects in a single file
#save.image("Robjects.2p5kb-TSS-1kb.2010-12-18-Vancouver.Rdata") # save all objects in a single file
#save.image("Robjects.2kb-TSS-1kb.2010-12-18-Vancouver.Rdata") # save all objects in a single file
#save.image("Robjects.1.5kb-TSS-1kb.2010-12-18-Vancouver.Rdata") # save all objects in a single file

##############################################################################
# Mark, 19 Dec 2010
downSamp <- function(d, f=0.5) {
  g <- rownames(d)
  newd <- matrix(0,nr=nrow(d), nc=ncol(d), dimnames=dimnames(d))
  for(j in 1:ncol(d)) {
    cat(".")
    s <- sample(g, f*sum(d[,j]), prob=d[,j], replace=TRUE)
    ts <- table(s)
    newd[names(ts),j] <- ts
  }
  cat("\n")
  newd
}



nSamps <- 10
cat("starting subsampling analysis,",nSamps,"iterations at each stage\n")

subs <- seq(0.3,1,by=0.1)
names(subs) <- subs

nDE <- lapply( subs, function(u) matrix(0,nrow=nrow(d$counts),ncol=nSamps) )
#nDEadjp <- lapply( subs, function(u) matrix(0,nrow=nrow(d$counts),ncol=nSamps) )


cat("starting downsampling at ",date(),"\n")

for(i in 1:length(subs)) {
  for(j in 1:nSamps) {
    DD <- downSamp(d$counts, subs[i])
    dSub <- DGEList( counts=DD, group=colnames(counts), lib.size=colSums(DD)*f )
    dSub <- estimateCommonDisp(dSub)
    deSub <- exactTest( dSub, pair=c("HS1235_ctl_2","HS1238_kd_2") )
    nDE[[i]][,j] <- abs( qnorm(deSub$table$p.value/2) )*sign(deSub$table$logFC)
    #nDEadjp[[i]][,j] <- deSub$table$adjp
  }
  cat("---------- downsampling at", subs[i],"finished ----------\n")
}

cat("finished downsampling at ",date(),"\n")


#write.table(nDE, "/projects/remc_bigdata/Karsan/R/20101219/nDE.lane-2.mapq10-noDups.knownGene.merge-500.2kb-TSS-1kb.mapq10_noDups.edgeR.satn_n10.20101219.txt", sep="\t", row.names=FALSE, quote=FALSE)
cat("generate a saturation boxplot for abs(z)\n")
zCut <- 4
z <- sapply( nDE, function(u) colSums(abs(u) > zCut) )

pdf("/projects/remc_bigdata/Karsan/R/20101219/downSampling.2-TSS-1kb.zGT4.pdf",w=6,h=6);
boxplot(z,at=subs,boxwex=.05,xlab="Percent of dataset",
        ylab="Nbr TSS rgns with abs(Z) > 4",
        ylim=range(unlist(z)),xlim=range(subs),
        main="Merged UCSC TSS rgns (-2000/+1000bp)")
grid()
dev.off()


#write.table(nDEadjp, "/projects/remc_bigdata/Karsan/R/20101219/nDE.lane-2.mapq10-noDups.knownGene.merge-500.2kb-TSS-1kb.mapq10_noDups.edgeR.satn_n10.20101219.txt", sep="\t", row.names=FALSE, quote=FALSE)
#cat("generate a saturation boxplot for adjp\n")
#adjpCut <- 0.0001
#ap <- sapply( nDEadjp, function(u) colSums(abs(u) > adjCut) )

#pdf("/projects/remc_bigdata/Karsan/R/20101219/downSampling.adjp-LE-1e-4.pdf",w=6,h=6);
#boxplot(ap,at=subs,boxwex=.05,xlab="Percent of dataset",
#        ylab="Number of TSS regions with adjp <= 0.0001",
#        ylim=range(unlist(ap)),xlim=range(subs))
#grid()
#dev.off()

cat("done\n")