Merge branch 'master' of https://github.com/friendsofstrandseq/pipeline

analysis/CN_BFG.R

+brs <- c(33998805, 34185186, 34211372, 35497754, 37452507, 39578189)
+breakpoints <- data.table(k = length(unique(round(brs/100000))), chrom = "chr10", bps = unique(round(brs/100000)))
+brfile <- "/home/maryam/research/hackathons/C7-7Jun/segmentation2/C7/manual-regions.txt"
+write.table(breakpoints, file = brfile, quote = F, sep = "\t", row.names = F)
+
+dir <- "/home/maryam/research/hackathons/C7-7Jun"
+countsFile <- file.path(dir, "100000_fixed.txt.gz")
+brFile <- file.path(dir, "segmentation2/C7/manual-regions.txt")
+infoFile <- file.path(dir, "100000_fixed.info")
+stateFile <- file.path(dir, "final.txt")
+
+# counts[chrom=="chr10" & start > 33998805 & end < 34185186]
+counts <- fread(paste("zcat", countsFile))
+segs <- fread(brFile)
+info <- fread(infoFile)
+strand <- fread(stateFile)
+
+counts[class=="None", class:="WW"]
+probs <- mosaiClassifierPrepare(counts, info, strand, segs)
+probs[, nb_r:=expected*nb_p/(1-nb_p)]
+
+maximumCN <- 200
+
+probs <- probs[, cbind(.SD, CN=1:maximumCN), 
+               by=.(sample, cell, chrom, start, end)]
+
+
+CNprobs <- probs[, .(cov=sum(W+C), 
+                     expected=sum(expected), 
+                     nb_p=nb_p[1], 
+                     nb_r=sum(nb_r)), 
+                 by=.(sample, chrom, start, end)]
+# first observation
+(CNprobs$cov / CNprobs$expected) * 2
+
+CNprobs <- CNprobs[, cbind(.SD, CN=1:maximumCN), 
+               by=.(sample, chrom, start, end)]
+
+
+# computing NB probs and calling most probable CN in probs table per segment per cell
+probs[, CN_ll:=dnbinom(x=C+W, size=CN*(nb_r/2), prob=nb_p)]
+cellCNcalls <- probs[, .SD[which.max(CN_ll)], 
+        by=.(sample, cell, chrom, start, end)]
+
+# computing NB probs and calling CN in CNprobs table per segment
+CNprobs[, CN_ll:=dnbinom(x=cov, size=CN*(nb_r/2), prob=nb_p)]
+CNcalls <- CNprobs[, .SD[which.max(CN_ll)], 
+        by=.(sample, chrom, start, end)]
+
+
+# removing extra columns for cleaning the data for plotting phylogenetic tree
+cellCNcalls <- cellCNcalls[start!=0]
+cellCNcalls[, `:=`(nb_p=NULL, class=NULL, C=NULL, W=NULL,expected=NULL, scalar=NULL, nb_r=NULL, CN_ll=NULL)]
+write.table(cellCNcalls, file.path(dir, "BFB_cell_CNs.table"))

analysis/plot_phylo_tree.R

+library(ggplot2)
+source("https://bioconductor.org/biocLite.R")
+biocLite("ggtree")
+library('ggtree')
+library("ape")
+
+dir <- "/home/maryam/research/hackathons/C7-7Jun"
+
+tab <- data.table::fread(file.path(dir, "BFB_cell_CNs.table"))
+CN.mat <- split(tab$CN, tab$cell) 
+CN.mat <- do.call(rbind, CN.mat)
+euc.dist <- dist(CN.mat, method = "euclidean")
+hc <- hclust(euc.dist, method = "ward.D2")
+plot(hc)
+
+plot(hclust(CN.mat, ))
+ggtree(as.phylo(hc)) + facet_grid()
+
+ggplot(data=cellCNcalls)
\ No newline at end of file