-
Notifications
You must be signed in to change notification settings - Fork 0
/
09_BQSR.sh
66 lines (59 loc) · 1.84 KB
/
09_BQSR.sh
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
#!/bin/bash
#Base quality score recalibration
#load modules
ml arch/haswell
ml arch/arch/haswell
ml modulefiles/arch/haswell
ml GATK/3.5-Java-1.8.0_71
#set variables
ref=/hpcfs/groups/acad_users/Refs/Homo_sapiens/GATK/GRCh38/GRCh38_full_analysis_set_plus_decoy_hla.fa
knownsites=/hpcfs/users/a1717363/mapping_resources/00-All.vcf.gz
indata=./05-markdup
outdata=./06-bqsr
recal=./06-bqsr/recal
tmp=./06-bsqr/tmp
samples=(Ka24 K24 K29 PUN67 PUN68 PUN76)
parallel=6
#recalibrate base quality scores
for (( i=0 ; i<${#samples[@]} ; i += $parallel ))
do
for j in $(seq $parallel)
do
individ=${samples[$i + $j - 1]}
if [[ -z $individ ]] ; then break 2; fi
java -Xmx16G -Djava.io.tmpdir=$tmp/ -jar $EBROOTGATK/GenomeAnalysisTK.jar -T BaseRecalibrator \
-nct 8 -R $ref -I $indata/${individ}_markdup.bam --knownSites $knownsites -o $recal/${individ}.table &
pids[$j]=$!
done
for pid in ${pids[@]}; do wait $pid; done
done
#write post BQSR bam
samples=(Ka24 K24 K29 PUN67 PUN68 PUN76)
parallel=6
for (( i=0 ; i<${#samples[@]} ; i += $parallel ))
do
for j in $(seq $parallel)
do
individ=${samples[$i + $j - 1]}
if [[ -z $individ ]] ; then break 2; fi
java -Xmx16G -Djava.io.tmpdir=$tmp/ -jar $EBROOTGATK/GenomeAnalysisTK.jar -T PrintReads \
-nct 16 -R $ref -I $indata/${individ}_markdup.bam -BQSR $recal/${individ}.table -o $outdata/${individ}_bqsr.bam &
pids[$j]=$!
done
for pid in ${pids[@]}; do wait $pid; done
done
#index bqsr bam
module load SAMtools/1.9-foss-2016b
samples=(Ka24 K24 K29 PUN67 PUN68 PUN76)
parallel=6
for (( i=0 ; i<${#samples[@]} ; i += $parallel ))
do
for j in $(seq $parallel)
do
individ=${samples[$i + $j - 1]}
if [[ -z $individ ]] ; then break 2; fi
samtools index $outdata/${individ}_bqsr.bam &
pids[$j]=$!
done
for pid in ${pids[@]}; do wait $pid; done
done