Start of codon correction update

bin/config.sh

@@ -252,6 +252,21 @@ KeepPreMappingReads=false
 MaxContigGappiness=0.05
 MinContigHitFrac=0.9
 
+# CodonCorrection.sh variables:
+virulign=~/VIRULIGN/virulign # testing - set to a more 'default' variable before finalising CC.sh
+# Shall we codon-correct the shiver output before remapping?
+CodonCorrect=true
+# The genes that can be analysed by VIRULIGN:
+genes="GAG POL VIF VPR VPU ENV NEF"
+# Which data should be output by VIRULIGN? Separated by ',' and enclosed by double quotes (see CodonCorrection.sh help for more details)
+VirulignOptions="Nucleotides,Mutations,AminoAcids"
+# Which genes should be corrected? Separated by ',' and enclosed by double quotes (see CodonCorrection.sh help for more details)
+GenesToCorrect="all"
+# Shall we delete temporary files within the CodonCorrection.sh output? 
+CC_DeleteTemp=false
+# Should key messages in the CodonCorrection output be printed in colour?
+CC_EnableColours=true
+
 # Suffixes we'll append to the sample ID for output files.
 # If you change the extension (whatever follows the dot) you might break
 # something.

bin/shiver_funcs.sh

@@ -27,6 +27,8 @@ Code_PrintSeqLengths="$ToolsDir/PrintSeqLengths.py"
 Code_AddSNPsToSeqs="$ToolsDir/AddAllPossibleSNPsToSeqs.py"
 Code_KeepBestLinesInDataFile="$ToolsDir/KeepBestLinesInDataFile.py"
 Code_ConvertFastqToFasta="$ToolsDir/ConvertFastqToFasta.py"
+Code_CodonCorrectionInit="$ToolsDir/CodonCorrectionInit.sh"
+Code_CodonCorrection="$ToolsDir/CodonCorrection.sh"
 
 # Only needed if GiveHXB2coords is set to true in the config file
 HXB2file=$(dirname "$ThisDir")'/data/external/B.FR.83.HXB2_LAI_IIIB_BRU.K03455.fasta'

bin/shiver_init.sh

@@ -147,3 +147,18 @@ mkdir -p "$IndividualRefDir" || \
 exit 1; }
 "$python" "$Code_SplitFasta" -G "$RefAlignment" "$IndividualRefDir" || { echo "Problem" \
 "splitting $RefAlignment into one file per sequence. Quitting." >&2 ; exit 1; }
+
+# Make the CodonCorrection init directory
+# NB This assumes that the /CodonCorrectionReferences is part of shiver to avoid changing the init arguments for now.
+# If this is not the case, then the reference file and gene coordinates will need to be supplied as an init argument.
+if [[ "$CodonCorrect" == "true" ]]; then
+  CC_RefFile="$ThisDir"/'CodonCorrectionReferences/HIV1_ALL_2020_genome_DNA_no_ambiguity_only_annotated_aligned_no_gap_columns_aligned.fasta'
+  CC_GeneCoords="$ThisDir"/'CodonCorrectionReferences/AnnotatableGenePositions.txt'
+  CC_Init="$OutDir"/'CodonCorrectionInit'
+  "$Code_CodonCorrectionInit" "$CC_RefFile" "$CC_Init" "$CC_GeneCoords" || \
+  { echo "Problem initialising CodonCorrection. Quitting." >&2 ; exit 1; }
+fi
+
+
+
+

bin/shiver_map_reads.sh

@@ -745,4 +745,12 @@ if [[ "$remap" == "true" ]]; then
     { echo 'Problem remapping to the consensus from the first round of mapping.'\
     'Quitting.' >&2 ; exit 1 ; }
   fi
+
+  # Codon-correct the round 2 consensus
+  CC_Init="$InitDir"/'CodonCorrectionInit'
+  CC_Output=./'CodonCorrection'
+  CC_Consensus="$NewSID"'_consensus_MinCov_'"$MinCov1"'_'"$MinCov2.fasta"
+  "$Code_CodonCorrection" "$CC_Consensus" "$virulign" "$CC_Init" "$VirulignOptions" "$GenesToCorrect" "$CC_Output" || \
+  { echo 'Problem codon-correcting the consensus before remapping. Quitting.' >&2 ; exit 1 ; }
+
 fi