From 5ca6ef8f01ed5253f506c2c5f94e255866278e23 Mon Sep 17 00:00:00 2001 From: Bryan Roessler Date: Fri, 2 Aug 2024 17:52:44 -0400 Subject: [PATCH] Modularize everything and refactor gtf --- workflow/qhtcp-workflow | 270 +++++++++++++++++++++++----------------- 1 file changed, 158 insertions(+), 112 deletions(-) diff --git a/workflow/qhtcp-workflow b/workflow/qhtcp-workflow index b3c37097..b8068300 100755 --- a/workflow/qhtcp-workflow +++ b/workflow/qhtcp-workflow @@ -57,16 +57,9 @@ # @option -y | --yes | --auto Assume yes answer to all questions (non-interactive mode) # @option -d | --debug Turn on extra debugging output # @option -h | --help Print help message and exit (overrides other options) - -DEBUG=1 # Turn debugging ON by default during development shopt -s extglob -# @section Libraries -# @description Change these variables to use different libraries -JAVA="${JAVA:-java}" -PYTHON="${PYTHON:-python3}" -PERL="${PERL:-perl}" -RSCRIPT="${RSCRIPT:-"R_LIBS_USER=~/R/qhtcp-workflow Rscript"}" +DEBUG=1 # Turn debugging ON by default during development # @section Help # @description Print a helpful message @@ -138,7 +131,7 @@ print_help() { # @description Creates array and switches from user input # parse_input() takes all of the arguments passed to the script parse_input() { - debug "Running: ${FUNCNAME[0]}" "$@" + debug "Running: ${FUNCNAME[0]} $*" long_opts="project:,module:,nomodule:,markdown,yes,auto,debug,help" short_opts="+p:m:n:ydh" @@ -267,7 +260,7 @@ random_words() { # @description Backup one or more files to an incremented .bk file # @exitcode backup iterator max 255 backup() { - debug "Running: ${FUNCNAME[0]}" "$@" + debug "Running: ${FUNCNAME[0]} $*" for f in "$@"; do [[ -e $f ]] || continue count=1 @@ -428,7 +421,7 @@ module install_dependencies # * install.packages(c('ontologyIndex', 'ggrepel', 'tidyverse', 'sos', 'openxlsx'), dep=TRUE) # @noargs install_dependencies() { - debug "Running: ${FUNCNAME[0]}" "$@" + debug "Running: ${FUNCNAME[0]} $*" # Dependency arrays depends_rpm=(graphviz pandoc pdftk-java gd-devel perl-CPAN shdoc nano rsync coreutils) @@ -487,7 +480,7 @@ install_dependencies() { "$RSCRIPT" -e "install.packages(c(\"$depends_r_str), dep=TRUE, repos=\"https://cloud.r-project.org\")" "$RSCRIPT" -e "BiocManager::install(\"${depends_bioc[0]}\")" - hash matlab &>/dev/null || echo "You will also need MATLAB installed for GUI modules" + hash "$MATLAB" &>/dev/null || echo "You will also need MATLAB installed for GUI modules" } @@ -735,11 +728,10 @@ easy() { # Add EASY directory to the Matlab path # If this does not work we can try changing the -sd argument and if that fails then pushing/popping debug "Adding EASY directory to the Matlab path" - hash matlab &>/dev/null && - matlab -nodisplay -nosplash -nodesktop -nojvm -batch "addpath('$EASY_DIR')" + "$MATLAB" -nodisplay -nosplash -nodesktop -nojvm -batch "addpath('$EASY_DIR')" # Launch matlab # matlab -nosplash -sd "$PROJECT_SCANS_DIR" -r "run $script" - matlab -nosplash -r "run $script" + "$MATLAB" -nosplash -r "run $script" fi } @@ -757,7 +749,7 @@ ezview() { # Make EZview dirs # Start EZview - matlab -nosplash -r "run $script" + "$MATLAB" -nosplash -r "run $script" fi } @@ -962,7 +954,7 @@ module remc # * This allows us to abstract the program away in script-run-workflow and treat it like a module # @arg $1 string studyInfo file remc() { - debug "Running: ${FUNCNAME[0]}" "$@" + debug "Running: ${FUNCNAME[0]} $*" # If any submodules fail the rest will not run, this is fundamental to module design # Remove leading && to run regardless @@ -993,37 +985,58 @@ remc() { module gtf # @section GTF # @description GTF module for QHTCP +# @arg $1 string output directory +# @arg $2 string gene_association.sgd +# @arg $3 string gene_ontology_edit.obo +# @arg $4 string ORF_List_Without_DAmPs.txt gtf() { debug "Running: ${FUNCNAME[0]}" - process_dir="$QHTCP_PROJECT_DIR/out/gtf/process" - function_dir="$QHTCP_PROJECT_DIR/out/gtf/function" - component_dir="$QHTCP_PROJECT_DIR/out/gtf/component" - out_dir="$QHTCP_PROJECT_DIR/out/gtf" + gtf_out_dir="${1:-$QHTCP_PROJECT_DIR/out/gtf}" + process_dir="$gtf_out_dir/process" + function_dir="$gtf_out_dir/function" + component_dir="$gtf_out_dir/component" + + gene_association_sgd="${2:-"$APPS_DIR/perl/gene_association.sgd"}" + gene_ontology_obo="${3:-"$APPS_DIR/perl/gene_ontology_edit.obo"}" + orf_list="${4:-"$APPS_DIR/perl/ORF_List_Without_DAmPs.txt"}" py_gtf_dcon \ "$process_dir" \ - "$out_dir" + "$gtf_out_dir" - # Perform operations on each directory in parallel + # Reproduce the function and components dirs from the process dir for d in "$function_dir" "$component_dir"; do debug "rsync -a $process_dir/ $d/" rsync -a "$process_dir/" "$d/" done + for d in "$process_dir" "$function_dir" "$component_dir"; do out_file="${d##*/}Results.txt" # Use the dirname to create each Results filename - pl_gtf "$d" "$out_dir" & # parallelize - py_gtf_concat "$d" "$out_dir" "$out_file" + shopt -s nullglob + txts=("$d"/*.txt) # glob all txt files from each dir + shopt -u nullglob + for txt in "${txts[@]}"; do + debug "pl_gtf_analyze -an $gene_association_sgd -as P -o $gene_ontology_obo -b $orf_list $txt" + pl_gtf_analyze \ + '-an' "$gene_association_sgd" \ + '-as' 'P' \ + '-o' "$gene_ontology_obo" \ + '-b' "$orf_list" \ + "$txt" + debug "pl_terms2tsv $txt" + pl_gtf_terms2tsv "$txt" + done + debug "py_gtf_concat $gtf_out_dir $out_file" + py_gtf_concat "$gtf_out_dir" "$out_file" done - r_compile_gtf "$out_dir" + r_compile_gtf "$gtf_out_dir" } module gta # @section GTA # @description GTA module for QHTCP -# NOTES -# * Heavily modified GTAtemplate.R # TODO # * # * @@ -1076,7 +1089,12 @@ gta() { zscores_file="$QHTCP_PROJECT_DIR/Exp$s/$zscores_file" if [[ -f $zscores_file ]]; then mkdir "$GTA_OUT_DIR/Exp$s" - r_gta "Exp$s" "$zscores_file" "$sgd_terms_tfile" "$sgd_features_file" "$GTA_OUT_DIR" + r_gta \ + "Exp$s" \ + "$zscores_file" \ + "$sgd_terms_tfile" \ + "$sgd_features_file" \ + "$GTA_OUT_DIR" fi done @@ -1117,6 +1135,7 @@ submodule r_gta # @description GTAtemplate R script # TODO: # * Is GTAtemplate.R actually a template? +# * Do we need to allow user customization? # # Files: # * gene_association.sgd: https://downloads.yeastgenome.org/curation/chromosomal_feature/gene_association.sgd @@ -1136,17 +1155,20 @@ submodule r_gta # @arg $5 string output directory # r_gta() { - debug "Running: ${FUNCNAME[0]}" "$@" + debug "Running: ${FUNCNAME[0]} $*" cat <<-EOF EOF - script="$APPS_DIR/r/GTAtemplate.R" - [[ -d $5 ]] || mkdir -p "$5" - debug "$RSCRIPT $script $*" - "$RSCRIPT" "$script" "$@" + "$RSCRIPT" "$script" \ + "$1" \ + "$2" \ + "$3" \ + "$4" \ + "$5" \ + "${@:6}" } @@ -1216,15 +1238,13 @@ submodule r_gta_heatmaps # @arg $7 string output directory # r_gta_heatmaps() { - debug "Running: ${FUNCNAME[0]}" "$@" + debug "Running: ${FUNCNAME[0]} $*" cat <<-EOF EOF - script="$APPS_DIR/r/TSHeatmaps5dev2.R" - [[ -d $5 ]] || mkdir -p "$5" - + [[ -d $7 ]] || mkdir -p "$7" debug "$RSCRIPT $script $*" "$RSCRIPT" "$script" "$@" } @@ -1265,7 +1285,7 @@ r_gta_heatmaps() { # script="ExpFrontend.m" # if ! ((YES)) && # ask "Start MATLAB to run $script? This requires a GUI."; then -# matlab -nosplash -r "$script" +# $MATLAB -nosplash -r "$script" # fi # } @@ -1292,7 +1312,9 @@ r_interactions() { script="$APPS_DIR/r/interactions.R" debug "$RSCRIPT $script" "$@" - "$RSCRIPT" "$script" "$@" + "$RSCRIPT" "$script" \ + "$1" \ + "${@:2}" # optional arguments } @@ -1302,10 +1324,14 @@ submodule r_join_interactions # @arg $2 string The sd value # @arg $3 string The studyInfo file r_join_interactions() { - debug "Running: ${FUNCNAME[0]}" "$@" + debug "Running: ${FUNCNAME[0]} $*" script="$APPS_DIR/r/joinInteractExps.R" debug "$RSCRIPT $script $*" - "$RSCRIPT" "$script" "$@" + "$RSCRIPT" "$script" \ + "$1" \ + "$2" \ + "$3" \ + "${@:4}" # optional arguments local out_files=("$1/REMcRdy_lm_only.csv" "$1/Shift_only.csv" "$1/parameters.csv") for f in "${out_files[@]}"; do [[ -f $f ]] || (echo "$f does not exist"; return 1) @@ -1314,15 +1340,18 @@ r_join_interactions() { submodule java_extract -# @description Jingyu's REMc java utility using file input file REMcRdy_lm_only.csv -# and output REMcRdy_lm_only.csv-finalTable.csv -# I'm not sure if the output dir is configurable so we can copy data around or push/pop +# @description Jingyu's REMc java utility +# Input file: +# * REMcRdy_lm_only.csv +# Output file: +# * REMcRdy_lm_only.csv-finalTable.csv +# NOTE: +# * Closed-source w/ hardcoded output directory, so have to pushd/popd to run (not ideal) # @arg $1 string The output directory java_extract() { debug "Running: ${FUNCNAME[0]}" classpath="$APPS_DIR/java/javaExtract.jar" - out_file="$1/REMcRdy_lm_only.csv-finalTable.csv" - + # backup REMcRdy_lm_only.csv-finalTable.csv if ! backup "$out_file"; then ask "Backup of $out_file failed, continue?" || return 1 @@ -1337,7 +1366,7 @@ java_extract() { debug "pushd && ${java_cmd[*]} && popd" pushd "$1" && "${java_cmd[@]}" && popd || return 1 - + out_file="$1/REMcRdy_lm_only.csv-finalTable.csv" [[ -f $out_file ]] || (echo "$out_file does not exist"; return 1) } @@ -1345,17 +1374,22 @@ java_extract() { submodule r_add_shift_values # @description Add shift values back to REMcRdy_lm_only.csv-finalTable.csv # and output "REMcWithShift.csv" for use with the REMc heat maps -# @arg $1 string The output csv file REMcRdy_lm_only.csv-finalTable.csv +# @arg $1 string REMcRdy_lm_only.csv-finalTable.csv # @arg $2 string Shift_only.csv -# @arg $3 string REMcWithShift.csv +# @arg $3 string StudyInfo.csv file # @arg $4 string The sd value r_add_shift_values() { - debug "Running: ${FUNCNAME[0]}" "$@" + debug "Running: ${FUNCNAME[0]} $*" script="$APPS_DIR/r/addShiftVals.R" - out_file="$QHTCP_PROJECT_DIR/REMcWithShift.csv" debug "$RSCRIPT $script $*" - "$RSCRIPT" "$script" "$@" - rm -f "REMcHeatmaps/"*.pdf + "$RSCRIPT" "$script" \ + "$1" \ + "$2" \ + "$3" \ + "$4" \ + "${@:5}" # optional arguments + rm -f "$QHTCP_PROJECT_DIR/REMcHeatmaps/"*.pdf + out_file="$QHTCP_PROJECT_DIR/REMcWithShift.csv" [[ -f $out_file ]] || (echo "$out_file does not exist"; return 1) } @@ -1365,92 +1399,87 @@ submodule r_create_heat_maps # @arg $1 string The final shift table (REMcWithShift.csv) # @arg $2 string The output directory r_create_heat_maps() { - debug "Running: ${FUNCNAME[0]}" "$@" + debug "Running: ${FUNCNAME[0]} $*" script="$APPS_DIR/r/createHeatMaps.R" - out_file="$QHTCP_PROJECT_DIR/compiledREMcHeatmaps.pdf" debug "$RSCRIPT $script $*" - "$RSCRIPT" "$script" "$@" + "$RSCRIPT" "$script" \ + "$1" \ + "$2" \ + "${@:3}" # optional arguments pdfs=(REMcHeatmaps/*.pdf) debug "pdftk ${pdfs[*]} output $out_file" pdftk "${pdfs[@]}" output "$out_file" + out_file="$2/compiledREMcHeatmaps.pdf" [[ -f $out_file ]] || (echo "$out_file does not exist"; return 1) } submodule r_heat_maps_homology # @description Execute createHeatMapsAll.R -# @arg $1 string The final shift table (REMcRdy_lm_only.csv-finalTable.csv) -# @arg $2 string The (Shift_only.csv) +# @arg $1 string REMcRdy_lm_only.csv-finalTable.csv +# @arg $2 string Shift_only.csv # @arg $3 string The (Yeast_Human_Homology_Mapping_biomaRt_18_0920.csv) # @arg $4 string The output directory r_heat_maps_homology() { - debug "Running: ${FUNCNAME[0]}" "$@" + debug "Running: ${FUNCNAME[0]} $*" script="$APPS_DIR/r/createHeatMapsHomology.R" - out_file="$4/compiledREMcHomologyHeatmaps.pdf" - # Clean old output + # Remove old output + debug "Removing old pdfs and csvs from $4" rm "$4/"*.{pdf,csv} "$RSCRIPT" "$script" \ - REMcWithShift.csv \ - Homology \ - "$APPS_DIR/r/170503_DAmPs_Only.txt" \ - Yeast_Human_Homology_Mapping_biomaRt_18_0920.csv + "$1" \ + "$2" \ + "$3" \ + "$4" \ + "${@:5}" # optional arguments pdfs=("$work_dir"/homology/*.pdf) pdftk "${pdfs[@]}" output "$out_file" - + out_file="$4/compiledREMcHomologyHeatmaps.pdf" [[ -f $out_file ]] || (echo "$out_file does not exist"; return 1) } submodule py_gtf_dcon # @description Perform python dcon portion of GTF +# Output file: +# * 1-0-0-finaltable.csv # @arg $1 string Directory to process # @arg $2 string Output directory name py_gtf_dcon() { - debug "Running: ${FUNCNAME[0]}" "$@" + debug "Running: ${FUNCNAME[0]} $*" script="$APPS_DIR/python/DconJG2.py" - debug "$PYTHON $SCRIPT $1 $2/" - "$PYTHON" "$SCRIPT" "$1" "$2/" + debug "$PYTHON $script $1 $2/" + "$PYTHON" "$script" \ + "$1" \ + "$2/" \ + "${@:3}" # optional arguments out_file="$2/1-0-0-finaltable.csv" [[ -f $out_file ]] || (echo "$out_file does not exist"; return 1) } -submodule pl_gtf -# @description Perl modules for GTF -# @arg $1 string Working directory -# @arg $2 string Output directory name to look for txt files -pl_gtf() { - debug "Running: ${FUNCNAME[0]}" "$@" - set1="$APPS_DIR/perl/ORF_List_Without_DAmPs.txt" - shopt -s nullglob - set2=("$2"/*.txt) # glob them all - shopt -u nullglob - for s2 in "${set2[@]}"; do - debug "pl_gtf_analyze $set1 $s2" - pl_gtf_analyze "$set1" "$s2" - debug "pl_terms2tsv $s2" - pl_gtf_terms2tsv "$s2" - done -} - - submodule pl_gtf_analyze # @description Perl analyze submodule # This seems weird to me because we're just overwriting the same data for all set2 members # https://metacpan.org/dist/GO-TermFinder/view/examples/analyze.pl # Is there a reason you need a custom version and not the original from cpan? -# @arg $1 string Set 1 TODO naming -# @arg $2 string Set 2 TODO naming +# @arg $1 string gene_association.sgd +# @arg $2 string gene_ontology_edit.obo +# @arg $3 string ORF_List_Without_DAmPs.txt +# @arg $4 string TODO txt to anaylze? I'm not sure what this is called pl_gtf_analyze() { - debug "Running: ${FUNCNAME[0]}" + debug "Running: ${FUNCNAME[0]} $*" script="$APPS_DIR/perl/analyze_v2.pl" - an="$APPS_DIR/perl/gene_association.sgd" - obo="$APPS_DIR/perl/gene_ontology_edit.obo" - debug "$PERL $script -an $an -as P -o $obo -b $1 $2" - "$PERL" "$script" -an "$an" -as P -o "$obo" -b "$1" "$2" + debug "$PERL $script $*" + "$PERL" "$script" \ + "$1" \ + "$2" \ + "$3" \ + "$4" \ + "${@:5}" # optional arguments } @@ -1460,7 +1489,7 @@ submodule pl_gtf_terms2tsv # # @arg $1 string Terms file TODO naming pl_gtf_terms2tsv() { - debug "Running: ${FUNCNAME[0]}" "$@" + debug "Running: ${FUNCNAME[0]} $*" script="$APPS_DIR/perl/terms2tsv.pl" debug "$PERL $script $1.terms > $1.tsv" "$PERL" "$script" "$1.terms" > "$1.tsv" @@ -1471,15 +1500,14 @@ submodule py_gtf_concat # @description Python concat submodule for GTF # Concat the process ontology outputs from the /REMcReady_lm_only folder # Probably should be translated to bash -# @arg $1 string working directory -# @arg $2 string output directory name to look for txt files -# @arg $3 string output file +# @arg $1 string output directory name to look for txt files +# @arg $2 string output file py_gtf_concat() { - debug "Running: ${FUNCNAME[0]}" + debug "Running: ${FUNCNAME[0]} $*" script="$APPS_DIR/python/concatGTFResults.py" - debug "$PYTHON $script $2/ $3" - "$PYTHON" "$script" "$2/" "$3" - [[ -f $3 ]] || (echo "$3 does not exist"; return 1) + debug "$PYTHON $script $1/ $2" + "$PYTHON" "$script" "$1/" "$2" + [[ -f $2 ]] || (echo "$2 does not exist"; return 1) } @@ -1487,7 +1515,7 @@ submodule r_compile_gtf # @description Compile GTF in R # @arg $1 string gtf output directory r_compile_gtf() { - debug "Running: ${FUNCNAME[0]}" + debug "Running: ${FUNCNAME[0]} $*" script="$APPS_DIR/r/CompileGTF.R" debug "$RSCRIPT $script $1" "$RSCRIPT" "$script" "$1" @@ -1599,9 +1627,21 @@ documentation() { # Most variables in main() are user configurable or can be overriden by env # @internal main() { - debug "Running: ${FUNCNAME[0]}" "$@" - - # Some global vars + debug "Running: ${FUNCNAME[0]} $*" + + # Libraries + declare -g JAVA="${JAVA:-$(which java 2>/dev/null)}" + declare -g PYTHON="${PYTHON:-$(which python3 2>/dev/null)}" + declare -g PERL="${PERL:-$(which perl 2>/dev/null)}" + declare -g RSCRIPT="${RSCRIPT:-$(which Rscript 2>/dev/null)}" + declare -g MATLAB="${MATLAB:-$(which matlab 2>/dev/null)}" + + # Use a custom R library + declare -gx R_LIBS_USER="$HOME/R/$SCRIPT_NAME" + [[ -d "$R_LIBS_USER" ]] || mkdir -p "$R_LIBS_USER" + + # Global vars + SCRIPT_NAME="${BASH_SOURCE[0]}" SCRIPT=$(realpath -s "${BASH_SOURCE[0]}") SCRIPT_DIR=$(dirname "$SCRIPT") APPS_DIR="$SCRIPT_DIR/apps" @@ -1610,13 +1650,18 @@ main() { DATE="$(date +%Y%m%d)" # change in EASYconsole.m to match 'hardcode' # Find a scans directory - # TODO change back for production, avoid actual scan dirs during testing # local scans_heirarchy=("./scans" "/mnt/data/scans" "/mnt/data/ExpJobs" "./scans") - local scans_heirarchy=( "$SCANS_DIR" "$SCRIPT_DIR/scans" "/mnt/data/scans" "templates/scans-demo" "./scans") + local scans_heirarchy=( + "$SCANS_DIR" + "$SCRIPT_DIR/scans" + "/mnt/data/scans" + "$SCRIPT_DIR/templates/scans-demo" + "$SCRIPT_DIR/scans" # fallback and create if others not found + ) [[ -z $SCANS_DIR ]] && for d in "${scans_heirarchy[@]}"; do if [[ -d $d ]]; then - declare -g SCANS_DIR="$d" + declare -gx SCANS_DIR="$d" fi done if ! [[ -d $SCANS_DIR ]]; then @@ -1644,7 +1689,8 @@ main() { if [[ -z $OUT_DIR ]]; then echo "No output directory found" - declare -g OUT_DIR="$SCRIPT_DIR/out" + declare -gx OUT_DIR="$SCRIPT_DIR/out" + # This is not something we do often, so ask if ask "Create $SCRIPT_DIR/out?"; then debug "mkdir $SCRIPT_DIR/out" mkdir "$SCRIPT_DIR/out"