# MESSAGE -----------------------------------------------------------------

#

# author: Yulin Lyu

# email: lvyulin@pku.edu.cn

#

# require: R whatever

#

# ---

# * 0. Setting ------------------------------------------------------------

settingList <- list(

  workDir = ".",

  trimmomaticDir = "app/Trimmomatic-0.39",

  trimAdapter = "TruSeq3-PE-2.fa",

  mapRef = "fasta/genome",

  picardDir = "app/picard.jar",

  bdg2bwDir = "app/bedGraphToBigWig",

  refLen = "fasta/genome.fa.fai"

)

# * 1. Load packages ------------------------------------------------------

setwd(settingList$workDir)

library(tidyverse)

library(magrittr)

library(glue)

# * 2. Preprocess ---------------------------------------------------------

setwd("0_fastq")

list.files(".")

from_file <- list.files(".", "")

to_file <- gsub("", "", from_file)

file.rename(from_file, to_file)

setwd("..")

# * * 2.0. Load data ------------------------------------------------------

file_name <- list.files("0_fastq", ".gz")

file_name <- list.files("2_trim", ".gz") # after qc

R1 <- grep("_1\\.", file_name, value = T)

R2 <- grep("_2\\.", file_name, value = T)

sample_name <- gsub("_1\\..*", "", R1)

# * * 2.1. QC for raw data ------------------------------------------------

dir.create("code")

qc_dir <- "1_qc"

qc_dir %>% dir.create()

qc_cmd <- glue("fastqc -o {qc_dir} -t 8 0_fastq/{file_name} &")

cat(qc_cmd[1])

write.table(c("#!/bin/bash\n", qc_cmd), glue("code/{qc_dir}.sh"), quote = F, row.names = F, col.names = F)

# multiqc -o 1_qc -f -n qc.raw 1_qc/*.zip

# * * 2.2. Trim -----------------------------------------------------------

dir.create("2_trim")

dir.create(".2_untrim")

trim <- settingList$trimmomaticDir

adapter <- settingList$trimAdapter

trim_cmd <- glue(

  "java -jar {trim}/trimmomatic-0.39.jar PE -threads 10 \\

  0_fastq/{R1} 0_fastq/{R2} \\

  2_trim/{sample_name}_1.trim.fastq.gz \\

  .2_untrim/{sample_name}_1.untrim.fastq.gz \\

  2_trim/{sample_name}_2.trim.fastq.gz \\

  .2_untrim/{sample_name}_2.untrim.fastq.gz \\

  ILLUMINACLIP:{trim}/adapters/{adapter}:2:30:7:1:true \\

  LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 > 2_trim/{sample_name}.trim.log 2>&1")

cat(trim_cmd[1])

trim_dir <- "2_trim"

trim_dir %>% str_c("code/", .) %>% dir.create()

setCMD(trim_cmd, str_c("code/", trim_dir), 1, F)

# * * 2.3. QC for trimmed data --------------------------------------------

# reload trimmed data first

qc_dir <- "3_qc"

qc_dir %>% dir.create()

qc_cmd <- glue("fastqc -o {qc_dir} -t 8 {trim_dir}/{file_name} &")

cat(qc_cmd[1])

write.table(c("#!/bin/bash\n", qc_cmd), glue("code/{qc_dir}.sh"), quote = F, row.names = F, col.names = F)

# multiqc -o 3_qc -f -n qc.trim 3_qc/*.zip

# * * 2.4. Map ------------------------------------------------------------

ref <- settingList$mapRef

map_dir <- "4_map"

map_dir %T>% dir.create() %>% str_c("code/", .) %>% dir.create()

map_cmd <- glue(

  "bowtie2 -p 20 --very-sensitive -X 2000 -x {ref} \\

  -1 {trim_dir}/{R1} -2 {trim_dir}/{R2} \\

  -S {map_dir}/{sample_name}.sam > {map_dir}/{sample_name}.log 2>&1")

cat(map_cmd[1])

setCMD(map_cmd, str_c("code/", map_dir), 1, T)

# * * 2.5. Sort -----------------------------------------------------------

sort_dir <- "5_sort"

sort_dir %T>% dir.create() %>% str_c("code/", .) %>% dir.create()

sort_cmd <- glue(

  "samtools view -@ 10 -bS {map_dir}/{sample_name}.sam | \\

  samtools sort -@ 10 > {sort_dir}/{sample_name}.bam")

cat(sort_cmd[1])

setCMD(sort_cmd, str_c("code/", sort_dir), 1, F)

# * * 2.6. Dedup ----------------------------------------------------------

picard <- settingList$picardDir

dedup_dir <- "6_dedup"

dedup_dir %T>% dir.create() %>% str_c("code/", .) %>% dir.create()

dedup_cmd <- glue(

  "java -Xms2g -Xmx8g -XX:ParallelGCThreads=8 -jar {picard} MarkDuplicates \\

  I={sort_dir}/{sample_name}.bam \\

  O={dedup_dir}/{sample_name}.dedup.bam \\

  M={dedup_dir}/{sample_name}.dedup.txt \\

  REMOVE_DUPLICATES=true")

cat(dedup_cmd[1])

setCMD(dedup_cmd, str_c("code/", dedup_dir), 1, F)

# * * 2.7. Filter ---------------------------------------------------------

fil_dir <- "7_filter"

fil_dir %>% dir.create()

filter_cmd <- glue(

  "samtools view -h -f 2 -q 30 {dedup_dir}/{sample_name}.dedup.bam | \\

  egrep -v 'MT|GL|JH' | samtools sort -@ 4 -O bam > {fil_dir}/{sample_name}.filter.bam &")

cat(filter_cmd[1])

write.table(c("#!/bin/bash\n", filter_cmd), glue("code/{fil_dir}.sh"), quote = F, row.names = F, col.names = F)

# * * 2.8. Finger ---------------------------------------------------------

fin_dir <- "8_finger"

fin_dir %>% dir.create()

index_cmd <- glue("samtools index {fil_dir}/{sample_name}.filter.bam {fil_dir}/{sample_name}.filter.bam.bai &")

cat(index_cmd)

write.table(c("#!/bin/bash\n", index_cmd), "code/8_index.sh", quote = F, row.names = F, col.names = F)

finger_group <- tapply(sample_name, gsub(".*-", "", sample_name), function(x) glue("{fil_dir}/{x}.filter.bam"))

finger_group %<>% map(~ paste(.x, collapse = " "))

finger_cmd <- glue("

  plotFingerprint -b {finger_group} -o {fin_dir}/{names(finger_group)}.finger.pdf &")

cat(finger_cmd[1])

write.table(c("#!/bin/bash\n", finger_cmd), glue("code/{fin_dir}.sh"), quote = F, row.names = F, col.names = F)

# * * 2.9. Peak -----------------------------------------------------------

peak_dir <- "9_peak"

peak_dir %>% dir.create()

input_sample <- grep("in", sample_name, value = T)

chip_sample <- grep("in", sample_name, value = T, invert = T)

input_sample %<>% rep(each = 3)

org <- "hs"

peak_cmd <- glue(

  "macs2 callpeak -t {fil_dir}/{chip_sample}.filter.bam -c {fil_dir}/{input_sample}.filter.bam \\

  -f BAMPE -g {org} --broad --keep-dup all -n {chip_sample} \\

  --outdir {peak_dir} > {peak_dir}/{chip_sample}.log 2>&1 &")

cat(peak_cmd[1])

write.table(c("#!/bin/bash\n", peak_cmd), glue("code/{peak_dir}.sh"), quote = F, row.names = F, col.names = F)

# * * 2.10. Bigwig --------------------------------------------------------

bw_dir <- "10_bw"

bw_dir %>% dir.create()

bw_cmd <- glue("

  bamCoverage --normalizeUsing RPKM -of bigwig \\

  -b {fil_dir}/{sample_name}.filter.bam \\

  -o {bw_dir}/{sample_name}.bw > {bw_dir}/{sample_name}.log 2>&1 &")

cat(bw_cmd[1])

bw_comp_cmd <- glue(

  "bamCompare -p 2 -of bigwig --operation ratio \\

  -b1 {fil_dir}/{chip_sample}.filter.bam \\

  -b2 {fil_dir}/{input_sample}.filter.bam \\

  -o {bw_dir}/{chip_sample}.comp.bw > {bw_dir}/{chip_sample}.comp.log 2>&1 &")

cat(bw_comp_cmd[1])

write.table(c("#!/bin/bash\n", bw_cmd), glue("code/{bw_dir}.sh"), quote = F, row.names = F, col.names = F)

write.table(c("#!/bin/bash\n", bw_comp_cmd), glue("code/{bw_dir}_comp.sh"), quote = F, row.names = F, col.names = F)

# * * 2.11. Heat ----------------------------------------------------------

chip_type <- chip_sample %>% str_replace(".*-", "")

bw_group <- tapply(chip_sample, chip_type, c) %>%

  map(~ glue("{bw_dir}/{.x}.comp.bw")) %>%

  map(~ str_c(.x, collapse = " "))

heat_dir <- "11_heat"

heat_dir %>% dir.create()

bed_file <- "all.bed"

out_prefix <- str_c("all_", names(bw_group))

mtx_cmd <- glue(

  "computeMatrix scale-regions -p 6 -R {heat_dir}/{bed_file} \\

  -S {bw_group} -o {heat_dir}/{out_prefix}_mtx.gz &"

)

cat(mtx_cmd[1])

heat_cmd <- glue(

  "plotHeatmap -m {heat_dir}/{out_prefix}_mtx.gz --colorMap RdBu_r \\

  -o {heat_dir}/{out_prefix}_heatmap.pdf \\

  --outFileNameMatrix {heat_dir}/{out_prefix}_value.txt \\

  --outFileSortedRegions {heat_dir}/{out_prefix}_region.bed &"

)

cat(heat_cmd[1])

write.table(c("#!/bin/bash\n", mtx_cmd), glue("code/{heat_dir}_mtx.sh"), quote = F, row.names = F, col.names = F)

write.table(c("#!/bin/bash\n", heat_cmd), glue("code/{heat_dir}.sh"), quote = F, row.names = F, col.names = F)

# * 3. Function -----------------------------------------------------------

setCMD <- function(cmd, dir = ".", sepN = 1, clu = F) {

  cmd %>% tapply(seq_along(.) %% sepN, c) %>% imap(~ {

    ifelse(clu, glue(

      "#!/bin/bash

      #SBATCH -J batch{.y}

      #SBATCH -o batch{.y}.%j.out

      #SBATCH -e batch{.y}.%j.err

      #SBATCH -p cn-long

      #SBATCH -N 1

      #SBATCH --ntasks-per-node=20

      #SBATCH --no-requeue

      #SBATCH -A hkdeng_g1

      #SBATCH --qos=hkdengcnl

      export PATH=/gpfs1/hkdeng_pkuhpc/lvyl/app/anaconda3/envs/lvyl/bin:$PATH"),

      "#!/bin/bash") %>%

      c(.x)}) %T>%

    iwalk(~ write.table(.x, glue("{dir}/batch{.y}.sh"), quote = F, row.names = F, col.names = F)) %>%

    names() %>% map_chr(~ glue("{head} {dir}/batch{.x}.sh {tail}",

                               head = ifelse(clu, "pkubatch", "sh"),

                               tail = ifelse(clu, "; sleep 1", "&"))) %>%

    c("#!/bin/bash", .) %>% as_tibble() %>%

    write_delim(glue("{dir}/submit.sh"), "\n", col_names = F)

}