-
Notifications
You must be signed in to change notification settings - Fork 7
/
linear_mappers_accuracy_kubernetes
executable file
·175 lines (154 loc) · 9.76 KB
/
linear_mappers_accuracy_kubernetes
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
(kubectl delete job xhchang-linear-primary || true) && kubectl apply -f - <<'EOF'
apiVersion: batch/v1
kind: Job
metadata:
name: xhchang-linear-primary
spec:
ttlSecondsAfterFinished: 86400
template:
spec:
containers:
- name: xhchang-map-primary
imagePullPolicy: Always
image: xhchang/vg:giraffe-paper
command:
- /bin/bash
- -c
- |
cd /tmp
export DEBIAN_FRONTEND=noninteractive
# Get the read combiner script
wget https://raw.githubusercontent.com/vgteam/giraffe-sv-paper/master/scripts/mapping/combine_reads.py
printf "graph\talgorithm\treads\tpairing\tcorrect\tmapq60\twrong_mapq60\tidentity\tscore\n" > report_minimap2.tsv
printf "graph\talgorithm\treads\tpairing\tcorrect\tmapq60\twrong_mapq60\tidentity\tscore\n" > report_bwa_mem.tsv
printf "graph\talgorithm\treads\tpairing\tcorrect\tmapq60\twrong_mapq60\tidentity\tscore\n" > report_bowtie2.tsv
printf "correct\tmq\tscore\taligner\n" > roc_stats_minimap2.tsv
printf "correct\tmq\tscore\taligner\n" > roc_stats_bwa_mem.tsv
printf "correct\tmq\tscore\taligner\n" > roc_stats_bowtie2.tsv
THREADS=16
#Get reference genomes
aws s3 cp s3://vg-k8s/profiling/data/hs38d1_renamed/ . --recursive
#Get xgs
aws s3 cp s3://vg-k8s/profiling/graphs/v3/for-NA19239/1000gplo/hs38d1/1000GPlo_hs38d1_filter.xg 1000gp.xg
aws s3 cp s3://vg-k8s/profiling/graphs/v2/for-NA19240/hgsvc/hs38d1/HGSVC_hs38d1.xg ./hgsvc.xg
for SPECIES in human ; do
case "${SPECIES}" in
yeast)
GRAPHS=(S288C)
READSETS=(DBVPG6044 DBVPG6765 N44 UWOPS034614 UWOPS919171 Y12 YPS138)
GRAPH_NAME=yeast
;;
human)
GRAPHS=(1000gp)
READSETS=(novaseq6000 hiseqxten hiseq2500)
GRAPH_NAME=hs38d1
;;
esac
for GRAPH in ${GRAPHS[@]} ; do
for READS in ${READSETS[@]} ; do
case ${GRAPH} in
1kg)
aws s3 cp s3://vg-k8s/profiling/reads/sim/for-NA19239/1kg/hs37d5/${READS}/out_sim_gbwt/sim.gam ./sim.gam
aws s3 cp s3://vg-k8s/profiling/reads/sim/for-NA19239/1kg/hs37d5/${READS}/out_sim_gbwt/sim.fq.gz ./sim.fq.gz
rm -f sim.fq
gunzip sim.fq.gz
;;
1000gp)
aws s3 cp s3://vg-k8s/profiling/reads/sim/for-NA19239/1000gp/hs38d1/liftover_nosegdups/${READS}/out_sim_gbwt/sim.gam ./sim.gam
aws s3 cp s3://vg-k8s/profiling/reads/sim/for-NA19239/1000gp/hs38d1/liftover_nosegdups/${READS}/out_sim_gbwt/sim.fq.gz ./sim.fq.gz
rm -f sim.fq
gunzip sim.fq.gz
;;
hgsvc)
aws s3 cp s3://vg-k8s/profiling/reads/sim/for-NA19240/hgsvc/grch38/${READS}/out_sim_gbwt/sim.gam ./sim.gam
aws s3 cp s3://vg-k8s/profiling/reads/sim/for-NA19240/hgsvc/grch38/${READS}/out_sim_gbwt/sim.fq.gz ./sim.fq.gz
rm -f sim.fq
gunzip sim.fq.gz
;;
S288C)
aws s3 cp s3://vg-k8s/profiling/reads/sim/yeast/sim-${READS}.gam ./sim.gam
aws s3 cp s3://vg-k8s/profiling/reads/sim/yeast/sim-${READS}.fq.gz ./sim.fq.gz
rm -f sim.fq
gunzip sim.fq.gz
;;
esac
sed 's/_1$//g' sim.fq | sed 's/_2$//g' > sim.paired.fq
for ALGORITHM in minimap2 bwa_mem bowtie2 ; do
for PAIRING in single paired ; do
if [[ ${ALGORITHM} == "minimap2" ]] ; then
if [[ ${PAIRING} == "paired" ]] ; then
minimap2 -ax sr --secondary=no -t ${THREADS} ${GRAPH_NAME}.fa.gz sim.paired.fq > mapped.bam
elif [[ ${PAIRING} == "single" ]] ; then
minimap2 -ax sr --secondary=no -t ${THREADS} ${GRAPH_NAME}.fa.gz sim.fq > mapped.bam
fi
elif [[ ${ALGORITHM} == "bwa_mem" ]] ; then
if [[ ${PAIRING} == "paired" ]] ; then
bwa mem -t ${THREADS} -p ${GRAPH_NAME}.fa sim.paired.fq > mapped.bam
elif [[ ${PAIRING} == "single" ]] ; then
bwa mem -t ${THREADS} -p ${GRAPH_NAME}.fa sim.fq > mapped.bam
fi
elif [[ ${ALGORITHM} == "bowtie2" ]] ; then
if [[ ${PAIRING} == "paired" ]] ; then
bowtie2 -t -p ${THREADS} -X 1065 -x ${GRAPH_NAME} --interleaved sim.paired.fq > mapped.bam
elif [[ ${PAIRING} == "single" ]] ; then
bowtie2 -t -p ${THREADS} -x ${GRAPH_NAME} -U sim.fq > mapped.bam
fi
fi
samtools view -F 2048 -b mapped.bam > mapped.primary.bam
samtools view -f 2048 -b mapped.bam > mapped.secondary.bam
vg inject -x ${GRAPH}.xg mapped.primary.bam > mapped.primary.gam
vg inject -x ${GRAPH}.xg mapped.secondary.bam > mapped.secondary.gam
if [[ ${PAIRING} == "paired" ]] ; then
vg view -aj mapped.primary.gam | sed 's/\/1/_1/g' | sed 's/\/2/_2/g' | vg view -aGJ - | vg annotate -m -x ${GRAPH}.xg -a - | vg gamcompare -r 100 -s - sim.gam 2> count | vg view -aj - > compared.primary.json
vg view -aj mapped.secondary.gam | sed 's/\/1/_1/g' | sed 's/\/2/_2/g' | vg view -aGJ - | vg annotate -m -x ${GRAPH}.xg -a - | vg gamcompare -r 100 - sim.gam| vg view -aj - > compared.secondary.json
elif [[ ${PAIRING} == "single" ]] ; then
vg annotate -m -x ${GRAPH}.xg -a mapped.primary.gam | vg gamcompare -s -r 100 - sim.gam 2> count | vg view -aj - > compared.primary.json
vg annotate -m -x ${GRAPH}.xg -a mapped.secondary.gam | vg gamcompare -r 100 - sim.gam | vg view -aj - > compared.secondary.json
fi
python ./combine_reads.py compared.primary.json compared.secondary.json compared.json
sed -i '/^$/d' compared.json
CORRECT_COUNT="$(grep correctly_mapped compared.json | wc -l)"
SCORE="$(sed -n '2p' count | sed 's/[^0-9\.]//g')"
MAPQ="$(grep mapping_quality\":\ 60 compared.json | wc -l)"
MAPQ60="$(grep -v correctly_mapped compared.json | grep mapping_quality\":\ 60 | wc -l)"
IDENTITY="$(jq '.identity' compared.json | awk '{sum+=$1} END {print sum/NR}')"
echo ${GRAPH} ${READS} ${PAIRING} ${SPEED} ${CORRECT_COUNT} ${MAPQ} ${MAPQ60} ${SCORE}
printf "${GRAPH}\t${ALGORITHM}\t${READS}\t${PAIRING}\t-\t${CORRECT_COUNT}\t${MAPQ}\t${MAPQ60}\t${IDENTITY}\t${SCORE}\n" >> report_${ALGORITHM}.tsv
jq -r '(if .correctly_mapped then 1 else 0 end|tostring) + "," + (.mapping_quality|tostring) + "," + (.score|tostring)' compared.json | sed 's/,/\t/g' | sed "s/$/\t${ALGORITHM}_${GRAPH}${READS}${PAIRING}/" | sed 's/single//g ; s/paired/-pe/g ; s/null/0/g' >> roc_stats_${ALGORITHM}.tsv
done
done
done
done
done
aws s3 cp report_minimap2.tsv s3://vg-k8s/users/xhchang/giraffe_experiments_2/liftover/report_minimap2_primary.tsv
aws s3 cp report_bwa_mem.tsv s3://vg-k8s/users/xhchang/giraffe_experiments_2/liftover/report_bwa_primary.tsv
aws s3 cp report_bowtie2.tsv s3://vg-k8s/users/xhchang/giraffe_experiments_2/liftover/report_bowtie2_primary.tsv
gzip roc_stats_minimap2.tsv
gzip roc_statsbwa_mem.tsv
gzip roc_stats_bowtie2.tsv
aws s3 cp roc_stats_minimap2.tsv.gz s3://vg-k8s/users/xhchang/giraffe_experiments_2/liftover/nosegdups/roc_stats_minimap2_primary.tsv.gz
aws s3 cp roc_stats_bwa_mem.tsv.gz s3://vg-k8s/users/xhchang/giraffe_experiments_2/liftover/nosegdups/roc_stats_bwa_primary.tsv.gz
aws s3 cp roc_stats_bowtie2.tsv.gz s3://vg-k8s/users/xhchang/giraffe_experiments_2/liftover/nosegdups/roc_stats_bowtie2_primary.tsv.gz
volumeMounts:
- mountPath: /tmp
name: scratch-volume
- mountPath: /root/.aws
name: s3-credentials
resources:
requests:
cpu: 24
memory: "120Gi"
ephemeral-storage: "150Gi"
limits:
cpu: 24
memory: "120Gi"
ephemeral-storage: "150Gi"
restartPolicy: Never
volumes:
- name: scratch-volume
emptyDir: {}
- name: s3-credentials
secret:
secretName: shared-s3-credentials
backoffLimit: 0
EOF