/
WGBSpre.R
200 lines (142 loc) · 6.34 KB
/
WGBSpre.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
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)
}