extensions/net.sf.basedb.reggie/trunk/src/net/sf/basedb/reggie/grid/scripts/bwa-mem2.sh

Code
Comments
Other
Rev Date Author Line
7086 31 Mar 23 nicklas 1 #!/bin/bash
7086 31 Mar 23 nicklas 2 ##
7183 17 May 23 nicklas 3 ## Pipeline script for the bwa-mem alignment step. 
7086 31 Mar 23 nicklas 4 ## 
7086 31 Mar 23 nicklas 5 ##
7086 31 Mar 23 nicklas 6 ## Environment variables that should be defined (in job.sh) before calling this script
7086 31 Mar 23 nicklas 7 ##  -FastqFolder: Path to folder with FASTQ files
7086 31 Mar 23 nicklas 8 ##  -AlignFolder: Folder for output. Will be created or cleared from existing files
7098 06 Apr 23 nicklas 9 ##  -REF_FASTA: Path to the reference FASTA file
7183 17 May 23 nicklas 10 ##  -BWAidx: Path prefix to the index used by bwa-mem
7183 17 May 23 nicklas 11 ##  -BWA_MEM: Path to bwa-mem program
7183 17 May 23 nicklas 12 ##  -BwaMemOptions: Other options for the bwa-mem program
7100 06 Apr 23 nicklas 13 ##  -ReadGroupStatic: Read group tags that are added to all alignments (eg. CN: and PL: tags)
7086 31 Mar 23 nicklas 14 ##  -SAMTOOLS: Path to samtools program
7102 06 Apr 23 nicklas 15 ##  -FixmateOptions: Options for the 'samtools fixmate' step
7102 06 Apr 23 nicklas 16 ##  -SortOptions: Options for the 'samtools sort' step
7102 06 Apr 23 nicklas 17 ##  -CalmdOptions: Options for the 'samtools calmd' step
7102 06 Apr 23 nicklas 18 ##  -MergeOptions: Options for the 'samtools merge' step
7102 06 Apr 23 nicklas 19 ##  -MarkdupOptions: Options for the 'samtools markdup' step
7086 31 Mar 23 nicklas 20 ##  -JAVA: Path to java
7086 31 Mar 23 nicklas 21 ##  -GATK: Path to GATK JAR file
7086 31 Mar 23 nicklas 22 ##  -DBSNP: Path to VCF file with SNP variants that should be genotyped
7086 31 Mar 23 nicklas 23 ##  -HaplotypeCallerOptions: Other options for the the HaplotypeCaller program
7086 31 Mar 23 nicklas 24 ##  -TMPDIR: Path to a local temporary working directory
7086 31 Mar 23 nicklas 25 ##
7086 31 Mar 23 nicklas 26
7088 03 Apr 23 nicklas 27 set -eo pipefail
7086 31 Mar 23 nicklas 28
7086 31 Mar 23 nicklas 29 ## Import utility functions
7086 31 Mar 23 nicklas 30 source ./reggie-utils.sh
7086 31 Mar 23 nicklas 31
7183 17 May 23 nicklas 32 ## If BWA_MEM is not set, we use BWA_MEM2 by default
7183 17 May 23 nicklas 33 if [ -z "${BWA_MEM}" ]; then
7183 17 May 23 nicklas 34   BWA_MEM="${BWA_MEM2}"
7183 17 May 23 nicklas 35 fi
7183 17 May 23 nicklas 36
7086 31 Mar 23 nicklas 37 ## Verify settings in options
7086 31 Mar 23 nicklas 38 rg_var_isdir "FastqFolder" "TMPDIR" "WD"
7183 17 May 23 nicklas 39 rg_var_isfile "BWA_MEM" "SAMTOOLS" "JAVA" "PICARD" "GATK" "REF_FASTA" "DBSNP"
7108 13 Apr 23 nicklas 40 rg_var_isset "BWAidx" "AlignFolder" "ReadLength"
7086 31 Mar 23 nicklas 41
7086 31 Mar 23 nicklas 42 ## Move to the temporary working directory and create subdirectories
7086 31 Mar 23 nicklas 43 cd ${TMPDIR}
7086 31 Mar 23 nicklas 44 mkdir -p fastq
7086 31 Mar 23 nicklas 45 mkdir -p aligned
7088 03 Apr 23 nicklas 46 mkdir -p merged
7105 11 Apr 23 nicklas 47 mkdir -p stats
7091 04 Apr 23 nicklas 48 mkdir -p tmp/sort
7091 04 Apr 23 nicklas 49 mkdir -p tmp/markdup
7091 04 Apr 23 nicklas 50 mkdir -p tmp/gatk
7086 31 Mar 23 nicklas 51 mkdir -p done
7086 31 Mar 23 nicklas 52
7088 03 Apr 23 nicklas 53 ## Read information about paired FASTQ files.
7088 03 Apr 23 nicklas 54 ## Columns in 'fastq_info.txt': 
7088 03 Apr 23 nicklas 55 ##  -FASTQ_R1
7088 03 Apr 23 nicklas 56 ##  -FASTQ_R2
7088 03 Apr 23 nicklas 57 NUM_LINES=`wc -l < ${WD}/fastq_info.txt`
7086 31 Mar 23 nicklas 58
7088 03 Apr 23 nicklas 59 if [ ! -f "done/copy.done" ]; then
7088 03 Apr 23 nicklas 60   CURRENT_LINE=0
7100 06 Apr 23 nicklas 61   while IFS=$'\t' read FASTQ1 FASTQ2 READ_GROUP; do
7088 03 Apr 23 nicklas 62     CURRENT_LINE=$(( ${CURRENT_LINE}+1 ))
7100 06 Apr 23 nicklas 63     rg_var_isset "FASTQ1" "FASTQ2" "READ_GROUP"
7088 03 Apr 23 nicklas 64     
7112 14 Apr 23 nicklas 65     PROGRESS=$(( ${CURRENT_LINE}*8/${NUM_LINES} ))
7112 14 Apr 23 nicklas 66     rg_progress ${PROGRESS} "Copying FASTQ files: ${CURRENT_LINE} of ${NUM_LINES} pairs"
7112 14 Apr 23 nicklas 67     
7088 03 Apr 23 nicklas 68     ## Rename FASTQ files so that we don't expose information in the BAM file
7088 03 Apr 23 nicklas 69     FASTQ1_LOCAL="fastq/${CURRENT_LINE}_R1.fastq.gz"
7088 03 Apr 23 nicklas 70     FASTQ2_LOCAL="fastq/${CURRENT_LINE}_R2.fastq.gz"
7088 03 Apr 23 nicklas 71     if [ "${JOB_DEBUG}" ]; then
7088 03 Apr 23 nicklas 72       # In debug mode we copy only parts of the FASTQ files
7088 03 Apr 23 nicklas 73       # The ' || true' is critical since 'zcat' will exit with 
7088 03 Apr 23 nicklas 74       # error when 'head' has closed the pipe
7112 14 Apr 23 nicklas 75       rg_time "STARTING Copy ${FASTQ1_LOCAL}"
7088 03 Apr 23 nicklas 76       (zcat ${FastqFolder}/${FASTQ1} || true) | head -n ${MaxFASTQLines} | bgzip -c > ${FASTQ1_LOCAL}
7112 14 Apr 23 nicklas 77       rg_time "STARTING Copy ${FASTQ2_LOCAL}"
7088 03 Apr 23 nicklas 78       (zcat ${FastqFolder}/${FASTQ2} || true) | head -n ${MaxFASTQLines} | bgzip -c > ${FASTQ2_LOCAL}
7088 03 Apr 23 nicklas 79     else
7112 14 Apr 23 nicklas 80       rg_time "STARTING Copy ${FASTQ1_LOCAL}"
7112 14 Apr 23 nicklas 81       cp ${FastqFolder}/${FASTQ1} ${FASTQ1_LOCAL}
7112 14 Apr 23 nicklas 82       rg_time "STARTING Copy ${FASTQ2_LOCAL}"
7112 14 Apr 23 nicklas 83       cp ${FastqFolder}/${FASTQ2} ${FASTQ2_LOCAL}
7088 03 Apr 23 nicklas 84     fi
7088 03 Apr 23 nicklas 85   done < ${WD}/fastq_info.txt
7112 14 Apr 23 nicklas 86
7112 14 Apr 23 nicklas 87   rg_time "ENDED Copy"
7088 03 Apr 23 nicklas 88   touch "done/copy.done"
7088 03 Apr 23 nicklas 89 fi
7088 03 Apr 23 nicklas 90
7088 03 Apr 23 nicklas 91 ## Align and samtools fixmate, sort and calmd 
7088 03 Apr 23 nicklas 92 if [ ! -f "done/align.done" ]; then
7088 03 Apr 23 nicklas 93   CURRENT_LINE=0
7100 06 Apr 23 nicklas 94   while IFS=$'\t' read FASTQ1 FASTQ2 READ_GROUP; do
7088 03 Apr 23 nicklas 95     CURRENT_LINE=$(( ${CURRENT_LINE}+1 ))
7088 03 Apr 23 nicklas 96
7112 14 Apr 23 nicklas 97     PROGRESS=$(( 10+${CURRENT_LINE}*30/${NUM_LINES} ))
7183 17 May 23 nicklas 98     rg_progress ${PROGRESS} "Running ${BWA_MEM##*/}: ${CURRENT_LINE} of ${NUM_LINES} FASTQ pairs (${NumThreads} threads)"
7088 03 Apr 23 nicklas 99
7088 03 Apr 23 nicklas 100     FASTQ1_LOCAL="fastq/${CURRENT_LINE}_R1.fastq.gz"
7088 03 Apr 23 nicklas 101     FASTQ2_LOCAL="fastq/${CURRENT_LINE}_R2.fastq.gz"
7088 03 Apr 23 nicklas 102     BAM_OUT="aligned/${CURRENT_LINE}.bam"
7088 03 Apr 23 nicklas 103   
7183 17 May 23 nicklas 104     rg_time "STARTING ${BWA_MEM##*/} ${BAM_OUT}"
7183 17 May 23 nicklas 105     ${BWA_MEM} mem ${BwaMemOptions} -R "@RG\\tID:${CURRENT_LINE}\\t${READ_GROUP}\\t${ReadGroupStatic}" \
7100 06 Apr 23 nicklas 106       -t ${NumThreads} ${BWAidx} ${FASTQ1_LOCAL} ${FASTQ2_LOCAL} \
7088 03 Apr 23 nicklas 107       | ${SAMTOOLS} fixmate ${FixmateOptions} --output-fmt bam,level=0 - - \
7091 04 Apr 23 nicklas 108       | ${SAMTOOLS} sort ${SortOptions} --output-fmt bam,level=0 -T tmp/sort - \
7088 03 Apr 23 nicklas 109       | ${SAMTOOLS} calmd ${CalmdOptions} --output-fmt bam,level=1 - ${REF_FASTA} \
7088 03 Apr 23 nicklas 110       > ${BAM_OUT}
7183 17 May 23 nicklas 111     rg_time "ENDED ${BWA_MEM##*/} ${BAM_OUT}"
7088 03 Apr 23 nicklas 112   
7088 03 Apr 23 nicklas 113   done < ${WD}/fastq_info.txt
7088 03 Apr 23 nicklas 114   touch "done/align.done"
7088 03 Apr 23 nicklas 115 fi
7088 03 Apr 23 nicklas 116
7088 03 Apr 23 nicklas 117 ## Merge and mark duplicates
7088 03 Apr 23 nicklas 118 if [ ! -f "done/merge.done" ]; then
7112 14 Apr 23 nicklas 119   rg_progress 50 "Merging ${NUM_LINES} BAM files and checking duplicates"
7088 03 Apr 23 nicklas 120
7109 13 Apr 23 nicklas 121   rg_time "STARTING merge and index"
7088 03 Apr 23 nicklas 122   ALL_BAM=`ls -w 0 aligned/*.bam`
7088 03 Apr 23 nicklas 123   ${SAMTOOLS} merge ${MergeOptions} --output-fmt bam,level=0 - ${ALL_BAM} \
7105 11 Apr 23 nicklas 124     | ${SAMTOOLS} markdup ${MarkdupOptions} \
7109 13 Apr 23 nicklas 125     -f stats/samtools_markdup.txt \
7105 11 Apr 23 nicklas 126     --output-fmt bam -T tmp/markdup - - \
7088 03 Apr 23 nicklas 127     > merged/alignment.bam
7106 12 Apr 23 nicklas 128   rg_time "ENDED merge"
7088 03 Apr 23 nicklas 129
7112 14 Apr 23 nicklas 130   rg_progress 60 "Creating BAM index"
7109 13 Apr 23 nicklas 131   ${SAMTOOLS} index -b -@ 8 merged/alignment.bam merged/alignment.bai
7109 13 Apr 23 nicklas 132   rg_time "ENDED index"
7105 11 Apr 23 nicklas 133
7088 03 Apr 23 nicklas 134   touch "done/merge.done"
7088 03 Apr 23 nicklas 135 fi
7088 03 Apr 23 nicklas 136
7105 11 Apr 23 nicklas 137 if [ ! -f "done/statistics.done" ]; then
7112 14 Apr 23 nicklas 138   rg_progress 65 "Calculating alignment statistics"
7106 12 Apr 23 nicklas 139   rg_time "STARTING statistics"
7110 13 Apr 23 nicklas 140   
7110 13 Apr 23 nicklas 141   ## Start all steps in the background...
7110 13 Apr 23 nicklas 142   ## Typical time ~15min
7110 13 Apr 23 nicklas 143   ${WD}/alignment_statistics.sh merged/alignment.bam > stats/alignment_statistics.txt &
7110 13 Apr 23 nicklas 144   ALIGNMENT_STATISTICS_PID=$!
7091 04 Apr 23 nicklas 145
7110 13 Apr 23 nicklas 146   ## Picard: CollectInsertSizeMetrics -- typical time ~25min
7105 11 Apr 23 nicklas 147   ${WD}/stdwrap.sh ${JAVA} ${JavaOptions} \
7105 11 Apr 23 nicklas 148     -jar ${PICARD} CollectInsertSizeMetrics \
7105 11 Apr 23 nicklas 149     -INPUT merged/alignment.bam \
7105 11 Apr 23 nicklas 150     -OUTPUT stats/picard_CollectInsertSizeMetrics.txt \
7105 11 Apr 23 nicklas 151     -Histogram_FILE stats/picard_CollectInsertSizeMetrics.pdf \
7106 12 Apr 23 nicklas 152     -QUIET true -VERBOSITY WARNING \
7110 13 Apr 23 nicklas 153     > picard_CollectInsertSizeMetrics.out &
7110 13 Apr 23 nicklas 154   INSERT_SIZE_PID=$!
7105 11 Apr 23 nicklas 155
7110 13 Apr 23 nicklas 156   ## Picard: CollectQualityYieldMetrics -- typical time ~30min
7106 12 Apr 23 nicklas 157   ${WD}/stdwrap.sh ${JAVA} ${JavaOptions} \
7110 13 Apr 23 nicklas 158     -jar ${PICARD} CollectQualityYieldMetrics \
7106 12 Apr 23 nicklas 159     -INPUT merged/alignment.bam \
7110 13 Apr 23 nicklas 160     -OUTPUT stats/picard_CollectQualityYieldMetrics.txt \
7106 12 Apr 23 nicklas 161     -QUIET true -VERBOSITY WARNING \
7110 13 Apr 23 nicklas 162     > picard_CollectQualityYieldMetrics.out &
7110 13 Apr 23 nicklas 163   QUALITY_YIELD_PID=$!
7106 12 Apr 23 nicklas 164
7110 13 Apr 23 nicklas 165   ## Picard: CollectAlignmentSummaryMetrics -- typical time ~1h
7106 12 Apr 23 nicklas 166   ${WD}/stdwrap.sh ${JAVA} ${JavaOptions} \
7106 12 Apr 23 nicklas 167     -jar ${PICARD} CollectAlignmentSummaryMetrics \
7106 12 Apr 23 nicklas 168     -INPUT merged/alignment.bam \
7106 12 Apr 23 nicklas 169     -OUTPUT stats/picard_CollectAlignmentSummaryMetrics.txt \
7106 12 Apr 23 nicklas 170     -QUIET true -VERBOSITY WARNING \
7110 13 Apr 23 nicklas 171     > picard_CollectAlignmentSummaryMetrics.out &
7110 13 Apr 23 nicklas 172   ALIGNMENT_SUMMARY_PID=$!
7106 12 Apr 23 nicklas 173
7111 13 Apr 23 nicklas 174   ## Picard: CollectWgsMetrics -- typical time ~1h15min
7106 12 Apr 23 nicklas 175   ${WD}/stdwrap.sh ${JAVA} ${JavaOptions} \
7110 13 Apr 23 nicklas 176     -jar ${PICARD} CollectWgsMetrics \
7106 12 Apr 23 nicklas 177     -INPUT merged/alignment.bam \
7110 13 Apr 23 nicklas 178     -OUTPUT stats/picard_CollectWgsMetrics.txt \
7110 13 Apr 23 nicklas 179     -R ${REF_FASTA} \
7110 13 Apr 23 nicklas 180     -READ_LENGTH ${ReadLength} \
7111 13 Apr 23 nicklas 181     -USE_FAST_ALGORITHM true \
7106 12 Apr 23 nicklas 182     -QUIET true -VERBOSITY WARNING \
7110 13 Apr 23 nicklas 183     > picard_CollectWgsMetrics.out &
7110 13 Apr 23 nicklas 184   WGS_PID=$!
7110 13 Apr 23 nicklas 185   
7110 13 Apr 23 nicklas 186   ## Wait for the steps to finish in the order they are expected to be completed
7110 13 Apr 23 nicklas 187   wait ${ALIGNMENT_STATISTICS_PID}
7110 13 Apr 23 nicklas 188   rg_time "ENDED alignment_statistics.sh"
7110 13 Apr 23 nicklas 189
7112 14 Apr 23 nicklas 190   rg_progress 70 "Calculating statistics: picard CollectInsertSizeMetrics"
7110 13 Apr 23 nicklas 191   wait ${INSERT_SIZE_PID}
7110 13 Apr 23 nicklas 192   rg_time "ENDED picard CollectInsertSizeMetrics"
7110 13 Apr 23 nicklas 193
7112 14 Apr 23 nicklas 194   rg_progress 75 "Calculating statistics: picard CollectQualityYieldMetrics"
7110 13 Apr 23 nicklas 195   wait ${QUALITY_YIELD_PID}
7110 13 Apr 23 nicklas 196   rg_time "ENDED picard CollectQualityYieldMetrics"
7110 13 Apr 23 nicklas 197
7112 14 Apr 23 nicklas 198   rg_progress 80 "Calculating statistics: picard CollectAlignmentSummaryMetrics"
7110 13 Apr 23 nicklas 199   wait ${ALIGNMENT_SUMMARY_PID}
7110 13 Apr 23 nicklas 200   rg_time "ENDED picard CollectAlignmentSummaryMetrics"
7110 13 Apr 23 nicklas 201   
7112 14 Apr 23 nicklas 202   rg_progress 85 "Calculating statistics: picard CollectWgsMetrics"
7110 13 Apr 23 nicklas 203   wait ${WGS_PID}
7110 13 Apr 23 nicklas 204   rg_time "ENDED picard CollectWgsMetrics"
7106 12 Apr 23 nicklas 205   rg_time "ENDED statistics"
7110 13 Apr 23 nicklas 206   
7089 04 Apr 23 nicklas 207   touch "done/statistics.done"
7089 04 Apr 23 nicklas 208 fi
7089 04 Apr 23 nicklas 209
7091 04 Apr 23 nicklas 210 ## Genotype QC
7091 04 Apr 23 nicklas 211 if [ ! -f "done/genotypeqc.done" ]; then
7091 04 Apr 23 nicklas 212   rg_progress 90 "Running HaplotypeCaller"
7106 12 Apr 23 nicklas 213   rg_time "STARTING Genotype QC"
7091 04 Apr 23 nicklas 214   ${WD}/stdwrap.sh ${JAVA} ${JavaOptions} \
7091 04 Apr 23 nicklas 215     -jar ${GATK} HaplotypeCaller \
7091 04 Apr 23 nicklas 216     -R ${REF_FASTA} \
7091 04 Apr 23 nicklas 217     -L ${DBSNP} --dbsnp ${DBSNP} --alleles ${DBSNP} \
7091 04 Apr 23 nicklas 218     --tmp-dir tmp/gatk \
7091 04 Apr 23 nicklas 219     ${HaplotypeCallerOptions} \
7091 04 Apr 23 nicklas 220     -I merged/alignment.bam \
7091 04 Apr 23 nicklas 221     -O merged/qc_genotype.vcf \
7091 04 Apr 23 nicklas 222     > gatk_HaplotypeCaller.out
7106 12 Apr 23 nicklas 223   rg_time "ENDED Genotype QC"
7102 06 Apr 23 nicklas 224     
7102 06 Apr 23 nicklas 225   # Remove '##contig' entries from the qc_genotype.vcf since they are not really needed
7102 06 Apr 23 nicklas 226   sed -i '/##contig/d' merged/qc_genotype.vcf
7091 04 Apr 23 nicklas 227   touch "done/genotypeqc.done"
7091 04 Apr 23 nicklas 228 fi
7089 04 Apr 23 nicklas 229
7091 04 Apr 23 nicklas 230
7089 04 Apr 23 nicklas 231 rg_progress 95 "Copying result files to project archive"
7112 14 Apr 23 nicklas 232 rg_time "STARTING Copy"
7112 14 Apr 23 nicklas 233
7105 11 Apr 23 nicklas 234 cp stats/*.txt ${WD}
7089 04 Apr 23 nicklas 235 mkdir -p ${AlignFolder}
7089 04 Apr 23 nicklas 236 rm -rf ${AlignFolder}/*
7089 04 Apr 23 nicklas 237 cp merged/* ${AlignFolder}
7105 11 Apr 23 nicklas 238 cp stats/* ${AlignFolder}
7105 11 Apr 23 nicklas 239 # cp *.out ${AlignFolder}
7112 14 Apr 23 nicklas 240 ls -1 ${AlignFolder}/* > ${WD}/files.out
7089 04 Apr 23 nicklas 241
7112 14 Apr 23 nicklas 242 rg_time "ENDED Copy"
7086 31 Mar 23 nicklas 243 rg_progress 99 "Analysis completed, cleaning up..."

Generated on Tue Nov 21 03:33:37 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)