forked from shiquan/PISA
/
sciATAC_Cusanovich_et_al._2018.wdl
121 lines (119 loc) · 3.18 KB
/
sciATAC_Cusanovich_et_al._2018.wdl
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
workflow main {
String root
String fastq1
String fastq2
String outdir
String refdir
String bwapath
String sambambapath
String ID
String ?runID
String config
String macspath
call makedir {
input:
Dir=outdir
}
call parseFastq {
input:
config=config,
fastq1=fastq1,
fastq2=fastq2,
outdir=makedir.Outdir,
runID=runID,
root=root
}
call fastq2bam {
input:
fastq=parseFastq.fastq,
outdir=outdir,
bwapath=bwapath,
refdir=refdir,
root=root
}
call sortBam {
input:
bam=fastq2bam.bam,
sambambapath=sambambapath,
root=root,
outdir=outdir
}
call callPeak {
input:
macspath=macspath,
bam=sortBam.sorted,
outdir=outdir,
root=root,
ID=ID
}
}
task makedir {
String Dir
command {
echo "[`date +%F` `date +%T`] workflow start" > ${Dir}/workflowtime.log
mkdir -p ${Dir}
mkdir -p ${Dir}/outs
mkdir -p ${Dir}/temp
}
output {
String Outdir="${Dir}"
}
}
task parseFastq {
String config
String fastq1
String fastq2
String outdir
String ?runID
String root
command {
${root}/SingleCellTools parse -t 15 -config ${config} -cbdis ${outdir}/temp/barcode_counts_raw.txt -run ${default="1" runID} -report ${outdir}/temp/sequencing_report.json ${fastq1} ${fastq2} | ${root}/SingleCellTools trim -report ${outdir}/temp/trim_report.json -mode Tn5 -p -t 3 > ${outdir}/temp/reads.fq
}
output {
String rawtable="${outdir}/temp/barcode_counts_raw.txt"
String fastq="${outdir}/temp/reads.fq"
String sequencingReport="${outdir}/temp/sequencing_report.json"
}
}
task fastq2bam {
String fastq
String outdir
String bwapath
String refdir
String root
command {
${bwapath} mem -t 20 -p ${refdir} ${fastq} | ${root}/SingleCellTools sam2bam -p -o ${outdir}/temp/aln.bam -report ${outdir}/temp/alignment_report.json -maln ${outdir}/temp/mito.bam /dev/stdin
}
output {
String bam="${outdir}/temp/aln.bam"
String alnReport="{outdir}/temp/alignment_report.json"
}
}
task sortBam {
String bam
String sambambapath
String root
String outdir
command {
${sambambapath} sort -t 20 -o ${outdir}/temp/sorted.bam ${outdir}/temp/aln.bam
${root}/SingleCellTools rmdup -tag CB -t 20 -o ${outdir}/temp/rmdup.bam ${outdir}/temp/sorted.bam
}
output {
String sorted="${outdir}/temp/rmdup.bam"
}
}
task callPeak {
String bam
String root
String outdir
String macspath
String ID
command <<<
${macspath} callpeak -t ${bam} -f BAM --keep-dup all --nomodel --shift -100 --extsize 200 -g mm -n ${ID} --outdir ${outdir}/outs
cut -f1,2,3 ${outdir}/outs/${ID}_peaks.narrowPeak > ${outdir}/temp/peak.bed
${root}/SingleCellTools anno -bed ${outdir}/temp/peak.bed -tag PK -o ${outdir}/outs/processed.bam ${bam}
${root}/SingleCellTools attrcnt -cb CB -tag PK -o ${outdir}/temp/readcount.report.txt ${outdir}/outs/processed.bam
awk '{if($1 !~ /CELL_BARCODE/ && $2>1000 && $3/$2>0.1){print $1;}}' ${outdir}/temp/readcount.report.txt > ${outdir}/temp/barcodes_called.txt
${root}/SingleCellTools count -tag CB -anno_tag PK -list ${outdir}/temp/barcodes_called.txt -o ${outdir}/outs/count_matrix.txt ${outdir}/outs/processed.bam
>>>
}