--- a
+++ b/a_PreprocessPipeline/WGBSpre.R
@@ -0,0 +1,200 @@
+
+# MESSAGE -----------------------------------------------------------------
+# 
+# author: Yulin Lyu
+# email: lvyulin@pku.edu.cn
+# 
+# require: R whatever
+# 
+# ---
+
+# * 0. Setting ------------------------------------------------------------
+
+settingList <- list(
+  workDir = ".",
+  mapRef = "fasta",
+  lambdaRef = "ref/lambda",
+  BSfilter = "BSfilter.R",
+  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")
+
+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. Map ------------------------------------------------------------
+
+ref <- settingList$mapRef
+
+file_dir <- "0_fastq"
+map_dir <- "2_map"
+map_dir %T>% dir.create() %>% str_c("code/", .) %>% dir.create()
+
+map_cmd <- glue(
+  "bismark --multicore 5 -X 700 --dovetail -p 2 --un --ambiguous --genome {ref} \\
+  {file_dir}/{R1} {file_dir}/{R2} \\
+  {map_dir} > {map_dir}/{sample_name}.log")
+cat(map_cmd[1])
+
+setCMD(map_cmd, str_c("code/", map_dir), 1, T)
+# multiqc -o 2_map -f -n map 2_map/*.txt
+
+# * * 2.3. Dedup ----------------------------------------------------------
+
+dedup_dir <- "3_dedup"
+dedup_dir %T>% dir.create() %>% str_c("code/", .) %>% dir.create()
+
+dedup_cmd <- glue(
+  "deduplicate_bismark -p -o {sample_name} --output_dir {dedup_dir} \\
+  {map_dir}/{sample_name}_1_bismark_bt2_pe.bam")
+cat(dedup_cmd[1])
+
+setCMD(dedup_cmd, str_c("code/", dedup_dir), 1, T)
+
+# * * 2.4. Meth -----------------------------------------------------------
+
+meth_dir <- "4_meth"
+meth_dir %T>% dir.create() %>% str_c("code/", .) %>% dir.create()
+
+meth_cmd <- glue(
+  "bismark_methylation_extractor -p --no_overlap --ignore 5 --ignore_r2 5 \\
+  --multicore 7 --bedGraph --cytosine_report --zero_based --CX --buffer_size 50% \\
+  --genome_folder {ref} -o {meth_dir} {dedup_dir}/{sample_name}.deduplicated.bam")
+cat(meth_cmd[1])
+
+setCMD(meth_cmd, str_c("code/", meth_dir), 1, T)
+
+# * * 2.5. Lambda ---------------------------------------------------------
+
+lambda_ref <- settingList$lambdaRef
+
+lambda_dir <- "5_lambda"
+lambda_dir %T>% dir.create() %>% str_c("code/", .) %>% dir.create()
+
+lambda_cmd <- glue(
+  "bismark --multicore 5 -X 700 --dovetail -p 2 --un --ambiguous --genome {lambda_ref} \\
+  {file_dir}/{R1} {file_dir}/{R2} \\
+  {lambda_dir} > {lambda_dir}/{sample_name}.log
+  
+  deduplicate_bismark -p -o {sample_name} --output_dir {lambda_dir} \\
+  {lambda_dir}/{sample_name}_1_bismark_bt2_pe.bam
+  
+  bismark_methylation_extractor -p --no_overlap --ignore 5 --ignore_r2 5 \\
+  --multicore 7 --bedGraph --cytosine_report --zero_based --CX --buffer_size 50% \\
+  --genome_folder {lambda_ref} -o {lambda_dir} {lambda_dir}/{sample_name}.deduplicated.bam")
+
+setCMD(lambda_cmd, str_c("code/", lambda_dir), 1, T)
+
+# * * 2.6. Extract --------------------------------------------------------
+
+extract_dir <- "6_extract"
+extract_dir %>% dir.create()
+
+extract_cmd <- glue(
+  "{awk_cmd} {meth_dir}/{sample_name}.deduplicated.CX_report.txt > {extract_dir}/{sample_name}.cg.txt &",
+  awk_cmd = "awk -F '\\t' '{if($4 + $5 >= 5 && $6 == \"CG\"){print $0}}'")
+cat(extract_cmd[1])
+
+write.table(c("#!/bin/bash\n", extract_cmd), glue("code/{extract_dir}.sh"), quote = F, row.names = F, col.names = F)
+
+# * * 2.7. Filter ---------------------------------------------------------
+
+p_files <- list.files(lambda_dir, full.names = T)
+p_text <- map(p_files, read_delim, delim = "\n")
+m <- p_text %>% map(~ .x[[1]][13:15]) %>% map(~ str_replace(.x, ".*\\t", "")) %>% map(~ sum(as.numeric(.x)))
+u <- p_text %>% map(~ .x[[1]][16:18]) %>% map(~ str_replace(.x, ".*\\t", "")) %>% map(~ sum(as.numeric(.x)))
+p <- map2_dbl(m, u, ~ {.x / (.x + .y)})
+
+filter_dir <- "7_filter"
+filter_dir %>% dir.create()
+
+BSfilter <- settingList$BSfilter
+
+filter_cmd <- glue(
+  "Rscript {BSfilter} -i {extract_dir}/{sample_name}.cg.txt \\
+  -d {filter_dir}/{sample_name}.dss.txt \\
+  -b {filter_dir}/{sample_name}.bdg \\
+  -p {p} > {filter_dir}/{sample_name}.log 2>&1 &")
+cat(filter_cmd[1])
+
+write.table(c("#!/bin/bash\n", filter_cmd), glue("code/{filter_dir}.sh"), quote = F, row.names = F, col.names = F)
+
+# * * 2.8 Bigwig ----------------------------------------------------------
+
+bw_dir <- "8_bw"
+bw_dir %>% dir.create()
+
+bdg2bw <- settingList$bdg2bwDir
+ref_len <- settingList$refLen
+
+bw_cmd <- glue("{bdg2bw} {filter_dir}/{sample_name}.bdg {ref_len} {bw_dir}/{sample_name}.bw &")
+cat(bw_cmd[1])
+
+write.table(c("#!/bin/bash\n", bw_cmd), glue("code/{bw_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/Bismark-0.22.3:$PATH
+      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)
+}