library(millefy)

# Path to bigWig files
bwfiles = Sys.glob(file.path(system.file("extdata", package="millefy"), "*.bw"))
bamfiles = Sys.glob(file.path(system.file("extdata/bam", package="millefy"), "*.bam"))
# FROM CODE at https://github.com/yuifu/millefy/blob/cf6cf0c8494df394b71702bfb928ef6bce2c7273/R/millefy-bamImport.R#L10 : bam_files = Sys.glob(file.path(system.file("extdata/bam", package="millefy"), "*.bam"))
print(bwfiles)
print(bamfiles)

# Group labels for bam files (same length as bamfiles)
groups = c("00h", "00h", "00h", "00h", "00h", "72h", "72h", "72h", "72h", "72h")

# Color labels for bigWig files (A named vector with the same length as the number of kinds of \\code{groups})
color_labels <- colorRampPalette(c("yellow", "red"))(length(unique(groups))+1)[1:length(unique(groups))]
names(color_labels)  <- unique(groups)
print(color_labels)

# Load gene models (It takes a little time)
path_gtf = system.file("extdata", "example.gtf", package="millefy")
dt_gtf_exon <- gtfToDtExon(path_gtf)

# Set tracks
## Single-cell track
max_value = 500
scTrackBw <- list(path_bam_files = bamfiles, groups = groups, group_colors = color_labels, max_value = max_value, isBw=FALSE, normFactors = calcBamNormFactors(bamfiles))

## Gene annotation track
geneTrack1 <- list(path_gtf = path_gtf, dt_gtf = dt_gtf_exon, label = "GENCODE")


# Prepare arguments for millefyPlot()
## List of tracks
tdlist <- list(scTrackBw, geneTrack1)

## List of track types
tt <- c("sc", "gene")

## List of track hights
heights = c(12, 2)

# Location to visualize
chr =  "chr19" # character
start = 5824708 # integer
end = 5845478 # integer

text_main = "mESC 00h, 72h (FROM Bam files!!!)"

l <- millefyPlot(track_data=tdlist, track_type=tt, heights=heights,
          sc_type = "heatmap",
          chr = chr, start = start, end = end,
          sc_avg = TRUE, sc_avg_height = 1,
          title = text_main)

system("pip install deeptools")

system("deeptools > verify_installed.txt")

cat(paste0(readLines("verify_installed.txt"), collapse="\n")) # based on [R - how can I dump contents of file to console output? `cat` equivalent in R](https://stackoverflow.com/a/59799268/8508004)

getwd()

system("cp -r /srv/conda/envs/notebook/lib/R/library/millefy/extdata/bam .")

system("mv bam bam2bw_extdata")

# based on https://gensoft.pasteur.fr/docs/deepTools/3.4.1/content/tools/bamCoverage.html?highlight=bamcoverage#usage-example-for-chip-seq
# and https://github.com/yuifu/millefy/blob/cf6cf0c8494df394b71702bfb928ef6bce2c7273/tutorial/Tutorial.md under [Converting BAM to BigWig files using deepTools](https://github.com/yuifu/millefy/blob/cf6cf0c8494df394b71702bfb928ef6bce2c7273/tutorial/Tutorial.md#1-1-converting-bam-to-bigwig-files-using-deeptools)
# and size of mouse chromsome 19 (learned it was mouse in Millefy publication) as `effectiveGenomeSize` because bamCoverage documentation says 'The effective genome size is the portion of the genome that is mappable'
# Note because stdout isn't displayed in R kernel in Jupyter, command worked out over in terminal in the same session
system("bamCoverage --bam bam2bw_extdata/RamDA_00h_A04.uniq.q40.chr19.bam -o bam2bw_extdata/RamDA_00h_A04.uniq.q40.chr19.binsize200.cbw --binSize 200 --normalizeUsing RPGC  --effectiveGenomeSize 61420004")

system("curl -OL https://gist.githubusercontent.com/fomightez/bab4bb92880b9545d20e2d0efc523a24/raw/b35a7e932da23be305b1952611d4ffe114c307ac/ZeroAnd72h_big_wigs.tar.gz")
system("tar xzf ZeroAnd72h_big_wigs.tar.gz")
system("mv *.bw bam2bw_extdata")

# Group labels for bam files (same length as bamfiles)
groups = c("00h", "00h", "00h", "00h", "00h", "72h", "72h", "72h", "72h", "72h")
converted_bamfiles = Sys.glob(file.path("bam2bw_extdata/", "*.bw"))

# Color labels for bigWig files (A named vector with the same length as the number of kinds of \\code{groups})
color_labels <- colorRampPalette(c("yellow", "red"))(length(unique(groups))+1)[1:length(unique(groups))]
names(color_labels)  <- unique(groups)
print(color_labels)

# Load gene models (It takes a little time)
path_gtf = system.file("extdata", "example.gtf", package="millefy")
dt_gtf_exon <- gtfToDtExon(path_gtf)

# Set tracks
## Single-cell track
max_value = 500 # NEEDS TO BE LOWER FOR THE CONVERTED BAM FILES
scTrackBw <- list(path_bam_files = converted_bamfiles, groups = groups, group_colors = color_labels, max_value = max_value, isBw=TRUE)

## Gene annotation track
geneTrack1 <- list(path_gtf = path_gtf, dt_gtf = dt_gtf_exon, label = "GENCODE")


# Prepare arguments for millefyPlot()
## List of tracks
tdlist <- list(scTrackBw, geneTrack1)

## List of track types
tt <- c("sc", "gene")

## List of track hights
heights = c(12, 2)

# Location to visualize
chr =  "chr19" # character
start = 5824708 # integer
end = 5845478 # integer

text_main = "mESC 00h, 72h (FROM Bam files!!!)"

converted_bamfiles

l <- millefyPlot(track_data=tdlist, track_type=tt, heights=heights,
          sc_type = "heatmap",
          chr = chr, start = start, end = end,
          sc_avg = TRUE, sc_avg_height = 1,
          title = text_main)