Quality_Assessment now generates a combined png image of all FASTQC plots

Clarify error messages for Adapter_Trimming and Read_Mapping
RTC and Indel_Realigner now support fixing quality scores that are too high
Automatically check for .bai indexing in Haplotype_Caller, RTC, and Indel_Realigner

Author: Aerin13

Config

@@ -1,5 +1,6 @@
 #!/bin/bash
 
+#   Config version 2.00
 #   More complete information on how to fill out
 #       this Config file can be found at:
 #       https://github.com/MorrellLAB/sequence_handling/wiki/Configuration
@@ -37,6 +38,11 @@ REF_GEN=
 #       Choose from: "true" or "false"
 BARLEY=
 
+#   Do the quality scores need to be adjusted for GATK? Default: false
+#       Change to true if you get errors from GATK like:
+#       "<Sample> appears to be using the wrong encoding for quality scores: we encountered an extremely high quality score"
+FIX_QUALITY_SCORES=false
+
 ############################################
 ##########     GBS_Demultiplex    ##########
 ############################################
@@ -368,11 +374,6 @@ FINISHED_BAM_LIST=
 #       For soybean: 0.001
 THETA=
 
-#   Do the quality scores need to be adjusted? Default: false
-#       Change to true only if you get errors from GATK like:
-#       "<Sample> appears to be using the wrong encoding for quality scores: we encountered an extremely high quality score"
-FIX_QUALITY_SCORES=false
-
 ############################################
 ##########     Genotype_GVCFs     ##########
 ############################################
@@ -401,9 +402,9 @@ REF_DICT=
 #       For soybean: 20
 NUM_CHR=
 
-#   Optional: Include the full file path to a list of the names of any and all scaffolds or parts of the reference not covered by the chromosomes above
-#       One scaffold name per line, must end in .intervals, leave blank if you wish to not call SNPs on non-chromosomal sequence
-#       Samtools style intervals are also acceptable, one per line (ex: chr1:100-200)
+#   Optional: Leave blank if you do not wish to call SNPs on non-chromosomal sequence
+#       Include the full file path to a list of the names of any and all scaffolds or parts of the reference not covered by the chromosomes above
+#       One scaffold name per line, must end in .intervals, SAMtools style intervals are also acceptable, one per line (ex: chr1:100-200)
 CUSTOM_INTERVALS=
 
 #   What is the sample ploidy? (integer value)
@@ -428,26 +429,26 @@ CHS_VCF_LIST=
 #   If not exome capture, put "NA"
 CAPTURE_REGIONS=
 
-#   What is the depth per sample (DP) cutoff?
+#   What is the depth per sample (DP) cutoff? (integer)
 #       If a sample's DP is below this threshold, it will count as a "bad" sample for that site.
 #       Recommended value: 5
 CHS_DP_PER_SAMPLE_CUTOFF=5
 
-#   What is the genotyping quality (GQ) cutoff?
+#   What is the genotyping quality (GQ) cutoff? (integer)
 #       If a sample's GQ is below this threshold, it will count as a "bad" sample for that site.
 #       Recommended value: 10th percentile of the raw GQ percentile table 
 CHS_GQ_CUTOFF=
 
-#   What is the maximum number of "bad" (low GQ, low DP, or missing genotype data) samples allowed at a site? 
+#   What is the maximum number of "bad" (low GQ, low DP, or missing genotype data) samples allowed at a site? (integer)
 #       Sites with more "bad" samples than this threshold will be filtered out.
 #       Recommended value: total number of samples * 0.2 (rounded to the nearest whole number)
 CHS_MAX_BAD=
 
-#   What is the maximum number of samples at a site that can be heterozygous?
+#   What is the maximum number of samples at a site that can be heterozygous? (integer)
 #       Recommended value: total number of samples * 0.9 (rounded to the nearest whole number)
 CHS_MAX_HET=
 
-#   What is the QUAL score cutoff?
+#   What is the QUAL score cutoff? (integer)
 #       Sites with a QUAL below this cutoff will be excluded.
 #       Recommended value: 40
 CHS_QUAL_CUTOFF=40
@@ -513,39 +514,39 @@ VF_VCF=${OUT_DIR}/Variant_Recalibrator/${PROJECT}_recalibrated.vcf
 #   If not exome capture, put "NA"
 VF_CAPTURE_REGIONS=
 
-#   What is the minimum number of reads needed to support a genotype?
+#   What is the minimum number of reads needed to support a genotype? (integer)
 #       Recommended value: 5
 MIN_DP=5
 
-#   What is the maximum number of reads allowed to support a genotype? (too many = gene duplication problems)
+#   What is the maximum number of reads allowed to support a genotype? (too many = gene duplication problems) (integer)
 #       Recommended value: 95th percentile of the raw DP per sample percentile table
 MAX_DP=
 
-#   What is the maximum percent deviation from 50/50 reference/alternative reads allowed in heterozygotes?
-#       For example, MAX_DEV=0.1 allows 60/40 ref/alt and also 40/60 ref/alt but not 70/30 or 30/70 ref/alt reads
+#   What is the maximum percent deviation from 50/50 reference/alternative reads allowed in heterozygotes? (decimal)
+#       For example, MAX_DEV=0.1 allows 60/40 ref/alt and also 40/60 ref/alt but not anything more skewed
 #       Recommended value: 0.1
 MAX_DEV=0.1
 
-#   What is the depth per sample (DP) cutoff?
+#   What is the depth per sample (DP) cutoff? (integer)
 #       If a sample's DP is below this threshold, it will count as a "bad" sample for that site.
 #       Recommended value: 5
 DP_PER_SAMPLE_CUTOFF=5
 
-#   What is the genotyping quality (GQ) cutoff?
+#   What is the genotyping quality (GQ) cutoff? (integer)
 #       If a sample's GQ is below this threshold, it will count as a "bad" sample for that site.
 #       Recommended value: 10th percentile of the raw GQ percentile table 
 GQ_CUTOFF=
 
-#   What is the maximum number of "bad" (low GQ, low DP, or missing genotype data) samples allowed at a site? 
+#   What is the maximum number of "bad" (low GQ, low DP, or missing genotype data) samples allowed at a site? (integer)
 #       Sites with more "bad" samples than this threshold will be filtered out.
 #       Recommended value: total number of samples * 0.2 (rounded to the nearest whole number)
 MAX_BAD=
 
-#   What is the maximum number of samples at a site that can be heterozygous?
+#   What is the maximum number of samples at a site that can be heterozygous? (integer)
 #       Recommended value: total number of samples * 0.9 (rounded to the nearest whole number)
 MAX_HET=
 
-#   What is the QUAL score cutoff?
+#   What is the QUAL score cutoff? (integer)
 #       Sites with a QUAL below this cutoff will be excluded.
 #       Recommended value: 40
 QUAL_CUTOFF=40
@@ -597,6 +598,11 @@ module load fastqc/0.11.5
 #FASTQC=
 #export PATH=${FASTQC}:${PATH}
 
+#   Do we have Riss-util installed?
+module load riss_util/1.0
+#RISS_UTIL=
+#export PATH=${RISS_UTIL}:${PATH}
+
 #   Do we have Seqqs installed?
 module load seqqs_ML/3d05375
 #SEQQS=

Handlers/Adapter_Trimming.sh

@@ -51,7 +51,7 @@ function trimAdapters() {
     #   Make a named for the trimmed sample
     local trimmed="${out}/${name}_ScytheTrimmed.fastq.gz"
     #   Trim the sample
-    scythe -a "${adapters}" -p "${prior}" -q "${platform}" --quiet "${toTrim}" | gzip -c > "${trimmed}"
+    (set -x; scythe -a "${adapters}" -p "${prior}" -q "${platform}" --quiet "${toTrim}" | gzip -c > "${trimmed}")
 }
 
 #   Export the function

Handlers/Indel_Realigner.sh

@@ -21,6 +21,7 @@ function Indel_Realigner() {
 	local intervals_list="$6" # Where is the list of intervals files?
     local lod="$7" # What is the LOD threshold?
     local entropy="$8" # What is the entropy threshold?
+    local qscores="$9" # Do we fix quality scores?
     declare -a intervals_array=($(grep -E ".intervals" "${intervals_list}")) # Put the intervals list into array format
     local current_intervals="${intervals_array[${PBS_ARRAYID}]}" # Get the intervals file for the current sample
     declare -a sample_array=($(grep -E ".bam" "${sample_list}")) # Put the sample list into array format
@@ -30,6 +31,8 @@ function Indel_Realigner() {
     mkdir -p "${out}"
     #   Change into the output directory
     cd "${out}"
+    if [[ "${qscores}" == true ]]
+    then
     #   Run GATK using the parameters given
     (set -x; java -Xmx"${memory}" -jar "${gatk}" \
 	    -T IndelRealigner \
@@ -38,7 +41,19 @@ function Indel_Realigner() {
 	    --LODThresholdForCleaning "${lod}" \
 	    --targetIntervals "${current_intervals}" \
 	    -I "${current}" \
+        --fix_misencoded_quality_scores \
 	    -o "${out}/${name}_realigned.bam")
+    else
+    #   Run GATK using the parameters given
+    (set -x; java -Xmx"${memory}" -jar "${gatk}" \
+	    -T IndelRealigner \
+	    -R "${reference}" \
+        --entropyThreshold "${entropy}" \
+	    --LODThresholdForCleaning "${lod}" \
+	    --targetIntervals "${current_intervals}" \
+	    -I "${current}" \
+	    -o "${out}/${name}_realigned.bam")
+    fi
 }
 
 #   Export the function

Handlers/Quality_Assessment.sh

@@ -73,6 +73,12 @@ function Quality_Assessment() {
     tail -n +2 "${out}/${project}_quality_summary_unfinished.txt" | sort >> "${out}/${project}_quality_summary.txt"
     # Remove the unsorted file
     rm "${out}/${project}_quality_summary_unfinished.txt"
+    # Change into the Zip_Files directory
+    cd "${out}/Zip_Files"
+    # Combine all plots into one file
+    fasterqc.pl -o "${out}/${project}_quality_plots.png"
+    # Remove the other file generated by fasterqc.pl
+    rm "${out}/Zip_Files/fasterqc.overrep.txt"
 }
 
 export -f Quality_Assessment

Handlers/Read_Mapping.sh

@@ -78,7 +78,7 @@ function Read_Mapping_Paired() {
     mkdir -p "${outDirectory}" # Make our outdirectory
     local memSettings=$(ParseBWASettings) # Assemble our settings for BWA mem
     local readGroupID=$(createReadGroupID "${sampleName}" "${project}" "${platform}") # Assemble our read group ID
-    bwa mem "${memSettings}" -R "${readGroupID}" "${reference}" "${forwardSample}" "${reverseSample}" > "${outDirectory}"/"${sampleName}".sam # Read map our sample
+    (set -x; bwa mem "${memSettings}" -v 2 -R "${readGroupID}" "${reference}" "${forwardSample}" "${reverseSample}" > "${outDirectory}"/"${sampleName}".sam)
 }
 
 #   Export the function
@@ -95,7 +95,7 @@ function Read_Mapping_Singles() {
     mkdir -p "${outDirectory}" # Make our outdirectory
     local memSettings=$(ParseBWASettings) # Assemble our settings for BWA mem
     local readGroupID=$(createReadGroupID "${sampleName}" "${project}" "${platform}") # Assemble our read group ID
-    bwa mem "${memSettings}" -R "${readGroupID}" "${reference}" "${sampleFile}" > "${outDirectory}"/"${sampleName}".sam # Read map our sample
+    (set -x; bwa mem "${memSettings}" -v 2 -R "${readGroupID}" "${reference}" "${sampleFile}" > "${outDirectory}"/"${sampleName}".sam)
 }
 
 #   Export the function

Handlers/Realigner_Target_Creator.sh

@@ -18,18 +18,31 @@ function Realigner_Target_Creator() {
     local gatk="$3" # Where is the GATK jar?
     local reference="$4" # Where is the reference sequence?
     local memory="$5" # How much memory can Java use?
+    local qscores="$6" # Do we fix quality scores?
     declare -a sample_array=($(grep -E ".bam" "${sample_list}")) # Put the sample list into array format
     local current="${sample_array[${PBS_ARRAYID}]}" # Pull out one sample to work on
     local name=$(basename ${current} .bam) # Get the name of the sample without the extension
     #	Make sure the out directory exists
     mkdir -p "${out}"
+    if [[ "${qscores}" == true ]]
+    then
     #   Run GATK using the parameters given
     (set -x; java -Xmx"${memory}" -jar "${gatk}" \
 	    -T RealignerTargetCreator \
 	    -R "${reference}" \
         -nt 1 \
 	    -I "${current}" \
+        --fix_misencoded_quality_scores \
 	    -o "${out}/${name}.intervals")
+    else
+    #   Run GATK using the parameters given
+    (set -x; java -Xmx"${memory}" -jar "${gatk}" \
+	    -T RealignerTargetCreator \
+	    -R "${reference}" \
+        -nt 1 \
+	    -I "${current}" \
+	    -o "${out}/${name}.intervals")
+    fi
 }
 
 #   Export the function

Handlers/SAM_Processing_Picard.sh

@@ -36,59 +36,65 @@ function SAM_Processing(){
     if [[ -z ${tmp} ]] # If tmp is left blank
     then
         #   Sort the SAM file and convert to BAM
-        (set -x; java -Xmx"${maxMem}" -jar ${picardJar} SortSam \
+        java -Xmx"${maxMem}" -jar ${picardJar} SortSam \
             INPUT="${SAMFile}" \
             OUTPUT="${outDirectory}/Intermediates/Sorted/${sampleName}_sorted.bam" \
             SO="coordinate" \
-            VALIDATION_STRINGENCY="SILENT")
+            VERBOSITY="WARNING" \
+            VALIDATION_STRINGENCY="SILENT"
         #   Deduplicate the BAM file
-        (set -x; java -Xmx"${maxMem}" -jar ${picardJar} MarkDuplicates \
+        java -Xmx"${maxMem}" -jar ${picardJar} MarkDuplicates \
             INPUT="${outDirectory}/Intermediates/Sorted/${sampleName}_sorted.bam" \
             OUTPUT="${outDirectory}/Intermediates/Deduplicated/${sampleName}_deduped.bam" \
             METRICS_FILE="${outDirectory}/Statistics/Deduplicated_BAM_Stats/${sampleName}_deduplicated.txt" \
             REMOVE_DUPLICATES="true" \
             ASSUME_SORTED="true" \
-            MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=${maxFiles})
+            VERBOSITY="WARNING" \
+            MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=${maxFiles}
         #   Add read group information to the BAM file
-        (set -x; java -Xmx"${maxMem}" -jar ${picardJar} AddOrReplaceReadGroups \
+        java -Xmx"${maxMem}" -jar ${picardJar} AddOrReplaceReadGroups \
             INPUT="${outDirectory}/Intermediates/Deduplicated/${sampleName}_deduped.bam" \
             OUTPUT="${outDirectory}/${sampleName}.bam" \
             RGID="${sampleName}" \
             RGLB="${sampleName}" \
             RGPL="${platform}" \
             RGPU="${sampleName}" \
-            RGSM="${sampleName}")
+            VERBOSITY="WARNING" \
+            RGSM="${sampleName}"
     else    # If a tmp is provided
         #   Make sure tmp exists
         mkdir -p ${tmp}
         #   Make sure we have write permissions for tmp
         if ! [[ -w "${tmp}" ]]; then echo "You do not have write permissions for ${tmp}, exiting..." >&2; exit 1; fi
         #   Sort the SAM file and convert to BAM
-        (set -x; java -Xmx"${maxMem}" -jar ${picardJar} SortSam \
+        java -Xmx"${maxMem}" -jar ${picardJar} SortSam \
             INPUT="${SAMFile}" \
             OUTPUT="${outDirectory}/Intermediates/Sorted/${sampleName}_sorted.bam" \
             SO="coordinate" \
             VALIDATION_STRINGENCY="SILENT" \
-            TMP_DIR="${tmp}")
+            VERBOSITY="WARNING" \
+            TMP_DIR="${tmp}"
         #   Deduplicate the BAM file
-        (set -x; java -Xmx"${maxMem}" -jar ${picardJar} MarkDuplicates \
+        java -Xmx"${maxMem}" -jar ${picardJar} MarkDuplicates \
             INPUT="${outDirectory}/Intermediates/Sorted/${sampleName}_sorted.bam" \
             OUTPUT="${outDirectory}/Intermediates/Deduplicated/${sampleName}_deduped.bam" \
             METRICS_FILE="${outDirectory}/Statistics/Deduplicated_BAM_Stats/${sampleName}_deduplicated.txt" \
             REMOVE_DUPLICATES="true" \
             ASSUME_SORTED="true" \
             MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=${maxFiles} \
-            TMP_DIR="${tmp}")
+            VERBOSITY="WARNING" \
+            TMP_DIR="${tmp}"
         #   Add read group information to the BAM file
-        (set -x; java -Xmx"${maxMem}" -jar ${picardJar} AddOrReplaceReadGroups \
+        java -Xmx"${maxMem}" -jar ${picardJar} AddOrReplaceReadGroups \
             INPUT="${outDirectory}/Intermediates/Deduplicated/${sampleName}_deduped.bam" \
             OUTPUT="${outDirectory}/${sampleName}.bam" \
             RGID="${sampleName}" \
             RGLB="${sampleName}" \
             RGPL="${platform}" \
             RGPU="${sampleName}" \
             RGSM="${sampleName}" \
-            TMP_DIR="${tmp}")
+            VERBOSITY="WARNING" \
+            TMP_DIR="${tmp}"
     fi
     #   Generate metrics on the finished BAM file
     samtools flagstat "${outDirectory}/${sampleName}.bam" > "${outDirectory}/Statistics/Finished_BAM_Stats/${sampleName}_finished.txt"

HelperScripts/utils.sh

@@ -2,7 +2,7 @@
 
 #   Check to make sure our samples exist
 function checkSamples() {
-    local sample_list="$1" # Sample ist
+    local sample_list="$1" # Sample list
     if [[ -f "${sample_list}" ]] # If the sample list exists
     then
         if [[ -r "${sample_list}" ]] # If the sample list is readable
@@ -11,22 +11,22 @@ function checkSamples() {
             do
                 if ! [[ -f "${sample}" ]] # If the sample doesn't exist
                 then
-                    echo "${sample} does not exist, exiting..." >&2 # Exit out with error
+                    echo "The sample ${sample} does not exist, exiting..." >&2 # Exit out with error
                     return 1
                 else
                     if ! [[ -r "${sample}" ]] # If the sample isn't readable
                     then 
-                        echo "${sample} does not have read permissions, exiting..." >&2
+                        echo "The sample ${sample} does not have read permissions, exiting..." >&2
                         return 1
                     fi
                 fi
             done
         else # If the sample isn't readable
-            echo "${sample_list} does not have read permissions, exiting..." >&2
+            echo "The sample list ${sample_list} does not have read permissions, exiting..." >&2
             return 1
         fi
     else # If the sample list doesn't exist
-        echo "$sample_list does not exist, exiting..." >&2 # Exit out with error
+        echo "The sample list $sample_list does not exist, exiting..." >&2 # Exit out with error
         return 1
     fi
 }
@@ -149,4 +149,22 @@ function createDict() {
 }
 
 #   Export the function
-export -f createDict
+export -f createDict
+
+#   Check to make sure our BAM files are indexed
+function checkBaiIndex() {
+    local sample_list="$1" # Sample list of BAM files, already checked by checkSamples (above)
+    for sample in $(cat "${sample_list}") # For each sample in the sample list
+    do
+        local basename=$(basename "${sample}" .bam)
+        local dirname=$(dirname "${sample}")
+        if [[ ! -f "${dirname}/${basename}.bai" && ! -f "${dirname}/${basename}.bam.bai" ]] # If the sample doesn't have a .bai index file of either naming convention
+        then
+            echo "The sample ${sample} does not have a .bai index, exiting..." >&2
+            exit 32
+        fi
+    done
+}
+
+#   Export the function
+export -f checkBaiIndex