bryan/hartman-server
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Generalize REMc analysis

Browse Source
This commit is contained in:
Bryan Roessler
2024-07-28 05:19:50 -04:00
parent c8eba4efd4
commit 8de5005055
8 changed files with 179 additions and 11994 deletions
Download Patch File Download Diff File Expand all files Collapse all files

115
workflow/apps/r/joinInteractExps.R
View File

@@ -1,3 +1,4 @@
#!/usr/bin/env Rscript
# JoinInteractExps.R
library(plyr)
@@ -46,40 +47,40 @@ print(length(inputFiles))
# Join the two files at a time as a function of how many inputFile, list the larger file first ? in this example X2 has the larger number of genes.
# If X1 has a larger number of genes, switch the order of X1 and X2
if(length(inputFiles)==2) {
X1 <- read.csv(file=inputFiles[1],stringsAsFactors = FALSE)
X2 <- read.csv(file=inputFiles[2],stringsAsFactors = FALSE)
X1 <- read.csv(file=inputFiles[1],stringsAsFactors=FALSE)
X2 <- read.csv(file=inputFiles[2],stringsAsFactors=FALSE)
X <- join(X1,X2,by="OrfRep")
OBH=X[,order(colnames(X))] #OrderByHeader
headSel= select(OBH, contains('OrfRep'), matches('Gene'), contains('Z_lm_K'), contains('Z_Shift_K'),contains('Z_lm_L'), contains('Z_Shift_L'))
headSel= select(headSel, -'Gene.1') #remove 'Gene.1 column
headSel2 = select(OBH, contains('OrfRep'), matches('Gene')) #Frame for interleaving Z_lm with Shift colums
headSel2 = select(headSel2, -'Gene.1') #remove 'Gene.1 column #Frame for interleaving Z_lm with Shift colums
headSel=select(OBH, contains('OrfRep'), matches('Gene'), contains('Z_lm_K'), contains('Z_Shift_K'),contains('Z_lm_L'), contains('Z_Shift_L'))
headSel=select(headSel, -'Gene.1') #remove 'Gene.1 column
headSel2=select(OBH, contains('OrfRep'), matches('Gene')) #Frame for interleaving Z_lm with Shift colums
headSel2=select(headSel2, -'Gene.1') #remove 'Gene.1 column #Frame for interleaving Z_lm with Shift colums
}else if(length(inputFiles)==3){
X1 <- read.csv(file=inputFiles[1],stringsAsFactors = FALSE) #exp1File,stringsAsFactors = FALSE)
X2 <- read.csv(file=inputFiles[2],stringsAsFactors = FALSE) #exp2File,stringsAsFactors = FALSE)
X3 <- read.csv(file=inputFiles[3],stringsAsFactors = FALSE) #exp3File,stringsAsFactors = FALSE)
X1 <- read.csv(file=inputFiles[1],stringsAsFactors=FALSE) #exp1File,stringsAsFactors=FALSE)
X2 <- read.csv(file=inputFiles[2],stringsAsFactors=FALSE) #exp2File,stringsAsFactors=FALSE)
X3 <- read.csv(file=inputFiles[3],stringsAsFactors=FALSE) #exp3File,stringsAsFactors=FALSE)
X <- join(X1,X2,by="OrfRep")
X <- join(X,X3,by="OrfRep")
OBH=X[,order(colnames(X))] #OrderByHeader
headSel= select(OBH, contains('OrfRep'), matches('Gene'), contains('Z_lm_K'), contains('Z_Shift_K'),contains('Z_lm_L'), contains('Z_Shift_L'))
headSel= select(headSel, -'Gene.1',-'Gene.2')
headSel2 = select(OBH, contains('OrfRep'), matches('Gene'))
headSel2 = select(headSel2, -'Gene.1',-'Gene.2')
headSel=select(OBH, contains('OrfRep'), matches('Gene'), contains('Z_lm_K'), contains('Z_Shift_K'),contains('Z_lm_L'), contains('Z_Shift_L'))
headSel=select(headSel, -'Gene.1',-'Gene.2')
headSel2=select(OBH, contains('OrfRep'), matches('Gene'))
headSel2=select(headSel2, -'Gene.1',-'Gene.2')
}else if(length(inputFiles)==4){
X1 <- read.csv(file=inputFiles[1],stringsAsFactors = FALSE) #exp1File,stringsAsFactors = FALSE)
X2 <- read.csv(file=inputFiles[2],stringsAsFactors = FALSE) #exp2File,stringsAsFactors = FALSE)
X3 <- read.csv(file=inputFiles[3],stringsAsFactors = FALSE) #exp3File,stringsAsFactors = FALSE)
X4 <- read.csv(file=inputFiles[4],stringsAsFactors = FALSE) #exp4File,stringsAsFactors = FALSE)
X1 <- read.csv(file=inputFiles[1],stringsAsFactors=FALSE) #exp1File,stringsAsFactors=FALSE)
X2 <- read.csv(file=inputFiles[2],stringsAsFactors=FALSE) #exp2File,stringsAsFactors=FALSE)
X3 <- read.csv(file=inputFiles[3],stringsAsFactors=FALSE) #exp3File,stringsAsFactors=FALSE)
X4 <- read.csv(file=inputFiles[4],stringsAsFactors=FALSE) #exp4File,stringsAsFactors=FALSE)
X <- join(X1,X2,by="OrfRep")
X <- join(X,X3,by="OrfRep")
X <- join(X,X4,by="OrfRep")
OBH=X[,order(colnames(X))] #OrderByHeader
headSel= select(OBH, contains('OrfRep'), matches('Gene'), contains('Z_lm_K'), contains('Z_Shift_K'),contains('Z_lm_L'), contains('Z_Shift_L'))
headSel= select(headSel, -'Gene.1',-'Gene.2',-'Gene.3')
headSel2 = select(OBH, contains('OrfRep'), matches('Gene'))
headSel2 = select(headSel2, -'Gene.1',-'Gene.2',-'Gene.3')
headSel=select(OBH, contains('OrfRep'), matches('Gene'), contains('Z_lm_K'), contains('Z_Shift_K'),contains('Z_lm_L'), contains('Z_Shift_L'))
headSel=select(headSel, -'Gene.1',-'Gene.2',-'Gene.3')
headSel2=select(OBH, contains('OrfRep'), matches('Gene'))
headSel2=select(headSel2, -'Gene.1',-'Gene.2',-'Gene.3')
}
#headSel$contains('Z_Shift') %>% replace_na(0.001)
@@ -88,10 +89,10 @@ i=0
for(i in 1:length(headers)){
if(grepl("Shift",headers[i])) {
headSel[headers[i]][is.na(headSel[headers[i]])] = 0.001
headSel[headers[i]][is.na(headSel[headers[i]])]=0.001
}
if(grepl("Z_lm_",headers[i])) {
headSel[headers[i]][is.na(headSel[headers[i]])] = 0.0001
headSel[headers[i]][is.na(headSel[headers[i]])]=0.0001
}
}
@@ -101,7 +102,7 @@ REMcRdy=select(headSel, contains('OrfRep'), matches('Gene'), contains('Z_lm_'))
shiftOnly=select(headSel, contains('OrfRep'), matches('Gene'), contains('Z_Shift'))
# Code to replace the numeric (.1 .2 .3) headers with experiment names from StudyInfo.txt
Labels <- read.csv(file= "../Code/StudyInfo.csv",stringsAsFactors = FALSE,sep= ",")
Labels <- read.csv(file="../Code/StudyInfo.csv",stringsAsFactors=FALSE,sep=",")
# Using Text search grepl to relabel headers
REMcRdyHdr=colnames(REMcRdy)
@@ -152,25 +153,25 @@ for(i in 3:length(colnames(REMcRdy))){
combI=cbind.data.frame(combI, REMcRdy[i])
}
Vec1= NA
Vec2= NA
Vec3= NA
Vec4= NA
Vec5= NA
Vec6= NA
Vec7= NA
Vec8= NA
Vec1=NA
Vec2=NA
Vec3=NA
Vec4=NA
Vec5=NA
Vec6=NA
Vec7=NA
Vec8=NA
if(length(REMcRdy)==6){
Vec1=abs(REMcRdy[,3])>=std
Vec2=abs(REMcRdy[,4])>=std
Vec3=abs(REMcRdy[,5])>=std
Vec4=abs(REMcRdy[,6])>=std
bolVec= Vec1 | Vec2 |Vec3 |Vec4
bolVec=Vec1 | Vec2 |Vec3 |Vec4
REMcRdyGT2=REMcRdy[bolVec,1:2]
REMcRdyGT2[ ,3:6]= REMcRdy[bolVec,3:6]
REMcRdyGT2[ ,3:6]=REMcRdy[bolVec,3:6]
shiftOnlyGT2=shiftOnly[bolVec,1:2]
shiftOnlyGT2[ ,3:6]= shiftOnly[bolVec,3:6]
shiftOnlyGT2[ ,3:6]=shiftOnly[bolVec,3:6]
}
if(length(REMcRdy)==8){
@@ -180,11 +181,11 @@ if(length(REMcRdy)==8){
Vec4=abs(REMcRdy[,6])>=std
Vec5=abs(REMcRdy[,7])>=std
Vec6=abs(REMcRdy[,8])>=std
bolVec= Vec1 | Vec2 |Vec3 | Vec4 |Vec5 |Vec6
REMcRdyGT2= REMcRdy[bolVec,1:2]
REMcRdyGT2[ ,3:8]= REMcRdy[bolVec,3:8]
shiftOnlyGT2= shiftOnly[bolVec,1:2]
shiftOnlyGT2[ ,3:8]= shiftOnly[bolVec,3:8]
bolVec=Vec1 | Vec2 |Vec3 | Vec4 |Vec5 |Vec6
REMcRdyGT2=REMcRdy[bolVec,1:2]
REMcRdyGT2[ ,3:8]=REMcRdy[bolVec,3:8]
shiftOnlyGT2=shiftOnly[bolVec,1:2]
shiftOnlyGT2[ ,3:8]=shiftOnly[bolVec,3:8]
}
if(length(REMcRdy)==10){
@@ -196,33 +197,33 @@ if(length(REMcRdy)==10){
Vec6=abs(REMcRdy[,8])>=std
Vec7=abs(REMcRdy[,9])>=std
Vec8=abs(REMcRdy[,10])>=std
bolVec= Vec1 | Vec2 |Vec3 |Vec4|Vec5|Vec6|Vec7|Vec8
REMcRdyGT2= REMcRdy[bolVec,1:2]
REMcRdyGT2[ ,3:10]= REMcRdy[bolVec,3:10]
shiftOnlyGT2= shiftOnly[bolVec,1:2]
shiftOnlyGT2[ ,3:10]= shiftOnly[bolVec,3:10]
bolVec=Vec1 | Vec2 |Vec3 |Vec4|Vec5|Vec6|Vec7|Vec8
REMcRdyGT2=REMcRdy[bolVec,1:2]
REMcRdyGT2[ ,3:10]=REMcRdy[bolVec,3:10]
shiftOnlyGT2=shiftOnly[bolVec,1:2]
shiftOnlyGT2[ ,3:10]=shiftOnly[bolVec,3:10]
}
if(std!=0){
REMcRdy= REMcRdyGT2 # [,2:length(REMcRdyGT2)]
shiftOnly= shiftOnlyGT2 # [,2:length(shiftOnlyGT2)]
REMcRdy=REMcRdyGT2 # [,2:length(REMcRdyGT2)]
shiftOnly=shiftOnlyGT2 # [,2:length(shiftOnlyGT2)]
}
if(std==0){
REMcRdy= REMcRdy # [,2:length(REMcRdy)]
shiftOnly= shiftOnly # [,2:length(shiftOnly)]
REMcRdy=REMcRdy # [,2:length(REMcRdy)]
shiftOnly=shiftOnly # [,2:length(shiftOnly)]
}
# R places hidden "" around the header names. The following
# is intended to remove those quote so that the "" do not blow up the Java REMc.
# Use ,quote=F in the write.csv statement to fix R output file.
#write.csv(combI,file = file.path(outDir,"CombinedKLzscores.csv"),row.names = FALSE)
write.csv(REMcRdy,file=file.path(outDir,"REMcRdy_lm_only.csv"),row.names = FALSE, quote=F)
write.csv(shiftOnly,file=file.path(outDir,"Shift_only.csv"),row.names = FALSE, quote=F)
#LabelStd <- read.table(file= "./parameters.csv",stringsAsFactors = FALSE,sep= ",")
#write.csv(combI,file=file.path(outDir,"CombinedKLzscores.csv"),row.names=FALSE)
write.csv(REMcRdy,file=file.path(outDir,"REMcRdy_lm_only.csv"),row.names=FALSE, quote=F)
write.csv(shiftOnly,file=file.path(outDir,"Shift_only.csv"),row.names=FALSE, quote=F)
#LabelStd <- read.table(file="./parameters.csv",stringsAsFactors=FALSE,sep=",")
LabelStd<- read.csv(file=studyInfo,stringsAsFactors = FALSE)
LabelStd<-read.csv(file=studyInfo,stringsAsFactors=FALSE)
print(std)
LabelStd[,4]= as.numeric(std)
write.csv(LabelStd,file=file.path(outDir,"parameters.csv"),row.names = FALSE)
write.csv(LabelStd,file=studyInfo,row.names = FALSE)
LabelStd[,4]=as.numeric(std)
write.csv(LabelStd,file=file.path(outDir,"parameters.csv"),row.names=FALSE)
write.csv(LabelStd,file=studyInfo,row.names=FALSE)

206
workflow/qhtcp-workflow
View File

@@ -43,7 +43,10 @@
# TODO:
# * Scripts should be made modular enough that they can be stored in the same dir
# * Don't cd in scripts
# * If you must, do it in a subshell at least!
# * Pass variables
# * Pass options
# * Pass arguments
# * Variable scoping is horrible right now
# * I wrote this sequentially and tried to keep track the best I could
# * Local vars have a higher likelihood of being lower case, global vars are UPPER
@@ -250,6 +253,19 @@ random_words() {
done
printf "%s_" "${arr[@]}" | sed 's/_$//'
}
# @description Backup one or more files to an incremented .bk file
backup() {
debug "Running: ${FUNCNAME[0]}" "$@"
for f in "$@"; do
[[ -f $f ]] || continue
count=1
while [[ -f $f.bk.$count ]]; do
count=$((count+1))
done
debug "rsync -a $f $f.bk.$count"
rsync -a "$f" "$f.bk.$count"
done
}
# @section Modules
@@ -260,6 +276,8 @@ random_words() {
# * Grouping multiple submodules (and modules) into a larger task
# * Dictating the ordering of multiple submodules
# * Modules should competently handle pushd and popd for their submodules if they do not reside in the SCANS/PROJECT_DIR
# * Apps and submodules should avoid changing directories
# * Pass input data from somewhere and output data somewhere
module install_dependencies
# @section Install dependencies
@@ -701,7 +719,6 @@ qhtcp() {
(( next_study_num=next_study_num+1 ))
done
# Now this is tricker than first appears
# Use initials from project not whoami
# Best I can do is first two letters of username
# See TODO in markdown
@@ -811,22 +828,33 @@ module remc
# * Don't cd within scripts, it's confusing
# * Use arguments to pass configuration variables
# * 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]}" "$@"
# Enter REMc directory to run the scripts there
pushd "$QHTCP_PROJECT_DIR/REMc" || return 1
# Run modules
# If any modules fail the rest will not run, this is fundamental to module design
# Remove trailing && to run regardless
# If any submodules fail the rest will not run, this is fundamental to module design
# Remove leading && to run regardless
# TODO can this be
r_join_interactions "$QHTCP_PROJECT_DIR/1-join_interactions" 2 "$1" "${@:2}" &&
java_extract &&
r_add_shift_values &&
r_heat_maps_zscores &&
r_heat_maps_homology &&
popd || return 1
r_join_interactions \
"$QHTCP_PROJECT_DIR/out" # output directory
2 \ % sd value
"$1" # studyInfo file
"${@:2}" \
&& java_extract \
"$QHTCP_PROJECT_DIR/out/" \
&& r_add_shift_values \
"$QHTCP_PROJECT_DIR/REMcRdy_lm_only.csv-finalTable.csv" \
"$QHTCP_PROJECT_DIR/Shift_only.csv" \
"$1" \
"$QHTCP_PROJECT_DIR/REMcWithShift.csv" \
&& r_create_heat_maps \
"$QHTCP_PROJECT_DIR/REMcWithShift.csv" \
"$QHTCP_PROJECT_DIR/out" \
&& r_heat_maps_homology \
"$QHTCP_PROJECT_DIR/REMcRdy_lm_only.csv-finalTable.csv" \
"$APPS_DIR/r/170503_DAmPs_Only.txt" \
"$APPS_DIR/r/Yeast_Human_Homology_Mapping_biomaRt_18_0920.csv" \
"$QHTCP_PROJECT_DIR/out/homology"
}
@@ -835,22 +863,27 @@ module gtf
# @description GTF module for QHTCP
gtf() {
debug "Running: ${FUNCNAME[0]}"
process_dir="GTF/Process"
function_dir="GTF/Function"
component_dir="GTF/Component"
out_dir="REMcRdy_lm_only"
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"
py_gtf_dcon "$process_dir" "$out_dir"
py_gtf_dcon \
"$process_dir" \
"$out_dir"
# Perform operations in each directory in parallel
# Perform operations on each directory in parallel
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
rsync -a "$process_dir/$out_dir" "$d"/
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"
done
r_compile_gtf
r_compile_gtf "$out_dir"
}
@@ -1144,9 +1177,9 @@ submodule r_join_interactions
r_join_interactions() {
debug "Running: ${FUNCNAME[0]}" "$@"
script="$APPS_DIR/r/joinInteractExps.R"
debug "$RSCRIPT $script"
debug "$RSCRIPT $script $*"
"$RSCRIPT" "$script" "$@"
local out_files=("REMcRdy_lm_only.csv" "Shift_only.csv" "parameters.csv")
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)
done
@@ -1156,23 +1189,28 @@ 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
# @arg $1 string The output directory
java_extract() {
debug "Running: ${FUNCNAME[0]}"
classpath="jingyuJava_1_7_extractLib.jar"
out_file="REMcRdy_lm_only.csv-finalTable.csv"
classpath="$APPS_DIR/java/javaExtract.jar"
out_file="$1/REMcRdy_lm_only.csv-finalTable.csv"
# backup REMcRdy_lm_only.csv-finalTable.csv
[[ -f $out_file ]] && mv "$out_file" "$out_file.bk"
if ! backup "$out_file"; then
ask "Backup of $out_file failed, continue?" || return 1
fi
java_cmd=(
"$JAVA" -Xms512m -Xmx2048m -Dfile.encoding=UTF-8 -classpath "$classpath" ExecMain
"REMcRdy_lm_only.csv"
"GeneByGOAttributeMatrix_nofiltering-2009Dec07.tab"
"ORF_List_Without_DAmPs.txt" 1 true true
"$QHTCP_PROJECT_DIR/REMcRdy_lm_only.csv"
"$APPS_DIR/java/GeneByGOAttributeMatrix_nofiltering-2009Dec07.tab"
"$APPS_DIR/java/ORF_List_Without_DAmPs.txt" 1 true true
)
debug "${java_cmd[@]}"
"${java_cmd[@]}"
debug "pushd && ${java_cmd[*]} && popd"
pushd "$1" && "${java_cmd[@]}" && popd || return 1
[[ -f $out_file ]] || (echo "$out_file does not exist"; return 1)
}
@@ -1180,23 +1218,31 @@ 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 $2 string Shift_only.csv
# @arg $3 string REMcWithShift.csv
# @arg $4 string The sd value
r_add_shift_values() {
debug "Running: ${FUNCNAME[0]}"
out_file="REMcHeatmaps/REMcWithShift.csv"
debug "$RSCRIPT AddShiftVals2.R"
"$RSCRIPT" AddShiftVals2.R
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
[[ -f $out_file ]] || (echo "$out_file does not exist"; return 1)
}
submodule r_heat_maps_zscores
# @description Execute REMcHeatmaps_zscores.R
r_heat_maps_zscores() {
debug "Running: ${FUNCNAME[0]}"
out_file="REMcHeatmaps/compiledREMcHeatmaps.pdf"
debug "$RSCRIPT REMcHeatmaps_zscores.R"
"$RSCRIPT" REMcHeatmaps_zscores.R
submodule r_create_heat_maps
# @description Execute createHeatMaps.R
# @arg $1 string The final shift table (REMcWithShift.csv)
# @arg $2 string The output directory
r_create_heat_maps() {
debug "Running: ${FUNCNAME[0]}" "$@"
script="$APPS_DIR/r/createHeatMaps.R"
out_file="$QHTCP_PROJECT_DIR/compiledREMcHeatmaps.pdf"
debug "$RSCRIPT $script $*"
"$RSCRIPT" "$script" "$@"
pdfs=(REMcHeatmaps/*.pdf)
debug "pdftk ${pdfs[*]} output $out_file"
pdftk "${pdfs[@]}" output "$out_file"
@@ -1205,32 +1251,26 @@ r_heat_maps_zscores() {
submodule r_heat_maps_homology
# @description Execute REMcHeatmaps_Z_lm_wDAmPs_andHomology_221212.R
# TODO this R script needs configurable output so we don't need to push/pop
# @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 $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]}"
work_dir="REMcHeatmapsWithHomology"
source_file="REMcHeatmaps/REMcWithShift.csv"
target_file="$work_dir/REMcWithShift.csv"
out_file="$work_dir/Homology/compiledREMcHomologyHeatmaps.pdf"
debug "rsync --archive $source_file $target_file"
rsync --archive "$source_file" "$target_file"
debug "Running: ${FUNCNAME[0]}" "$@"
script="$APPS_DIR/r/createHeatMapsHomology.R"
out_file="$4/compiledREMcHomologyHeatmaps.pdf"
# Clean old output
rm "$work_dir/Homology/"*.{pdf,csv}
rm "$4/"*.{pdf,csv}
pushd "$work_dir" || return 1
# TODO This
"$RSCRIPT" \
REMcHeatmaps_Z_lm_wDAmPs_andHomology_221212.R \
"$RSCRIPT" "$script" \
REMcWithShift.csv \
Homology \
17_0503_DAmPs_Only.txt \
"$APPS_DIR/r/170503_DAmPs_Only.txt" \
Yeast_Human_Homology_Mapping_biomaRt_18_0920.csv
popd || return 1
pdfs=("$work_dir"/Homology/*.pdf)
pdfs=("$work_dir"/homology/*.pdf)
pdftk "${pdfs[@]}" output "$out_file"
[[ -f $out_file ]] || (echo "$out_file does not exist"; return 1)
@@ -1242,11 +1282,11 @@ submodule py_gtf_dcon
# @arg $1 string Directory to process
# @arg $2 string Output directory name
py_gtf_dcon() {
debug "Running: ${FUNCNAME[0]}"
in_file="REMcRdy_lm_only.csv-finalTable.csv"
out_file="$1/$2/1-0-0-finaltable.csv"
debug "$PYTHON DconJG2.py $in_file $1/"
"$PYTHON" DconJG2.py "$in_file" "$1/"
debug "Running: ${FUNCNAME[0]}" "$@"
script="$APPS_DIR/python/DconJG2.py"
debug "$PYTHON $SCRIPT $1 $2/"
"$PYTHON" "$SCRIPT" "$1" "$2/"
out_file="$2/1-0-0-finaltable.csv"
[[ -f $out_file ]] || (echo "$out_file does not exist"; return 1)
}
@@ -1257,21 +1297,16 @@ submodule pl_gtf
# @arg $2 string Output directory name to look for txt files
pl_gtf() {
debug "Running: ${FUNCNAME[0]}" "$@"
pushd "$1" || return 1
set1="ORF_List_Without_DAmPs.txt"
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_analyze $set1 $s2"
debug "pl_gtf_analyze $set1 $s2"
pl_gtf_analyze "$set1" "$s2"
debug "pl_terms2tsv $s2"
pl_gtf_terms2tsv "$s2"
done
popd || return 1
}
@@ -1284,21 +1319,22 @@ submodule pl_gtf_analyze
# @arg $2 string Set 2 TODO naming
pl_gtf_analyze() {
debug "Running: ${FUNCNAME[0]}"
script="analyze_v2.pl"
an="gene_association.sgd"
out_file="gene_ontology_edit.obo"
debug "$PERL $script -an $an -as P -o $out_file -b $1 $2"
"$PERL" "$script" -an "$an" -as P -o "$out_file" -b "$1" "$2"
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"
}
submodule pl_gtf_terms2tsv
# @description Perl terms2tsv submodule
# Probably should be translated to shell/python
#
# @arg $1 string Terms file TODO naming
pl_gtf_terms2tsv() {
debug "Running: ${FUNCNAME[0]}"
script="terms2tsv_v4.pl"
debug "Running: ${FUNCNAME[0]}" "$@"
script="$APPS_DIR/perl/terms2tsv.pl"
debug "$PERL $script $1.terms > $1.tsv"
"$PERL" "$script" "$1.terms" > "$1.tsv"
}
@@ -1313,23 +1349,21 @@ submodule py_gtf_concat
# @arg $3 string output file
py_gtf_concat() {
debug "Running: ${FUNCNAME[0]}"
pushd "$1" || return 1
script="Concatenate_GTF_results.py"
script="$APPS_DIR/python/concatGTFResults.py"
debug "$PYTHON $script $2/ $3"
"$PYTHON" "$script" "$2/" "$3"
[[ -f $3 ]] || (echo "$3 does not exist"; return 1)
popd || return 1
}
submodule r_compile_gtf
# @description Compile GTF in R
# @noargs
# @arg $1 string gtf output directory
r_compile_gtf() {
debug "Running: ${FUNCNAME[0]}"
script="CompileGTF.R"
debug "$RSCRIPT $script"
"$RSCRIPT" "$script"
script="$APPS_DIR/r/CompileGTF.R"
debug "$RSCRIPT $script $1"
"$RSCRIPT" "$script" "$1"
}

71
workflow/templates/qhtcp/REMc/GTF/Concatenate_GTF_results.py
View File

@@ -1,71 +0,0 @@
#!/usr/bin/env python
# This code is to concatenate the batch GO Term Finder results (.tsv) generated from batch GTF perl code(Chris Johnson, U of Tulsa) into a list table
import sys, os, string, glob
# return the file list
def list_all_files(Path):
list_all_files = []
list_all_files = glob.glob(Path +'/*.txt.tsv')
return list_all_files
# Main function
try:
data_file_Path = sys.argv[1]
output_file_Path = sys.argv[2]
except:
print ('Usage: python Concatenate_GTF_results.py /datasetPath /outputFilePath_and_Name')
print ('Data file not found, error in given directory')
sys.exit(1)
# Open the output file
try:
output = open(output_file_Path, 'w')
except OSError:
print ('output file error')
# get all the GTF result files in given directory
File_list = []
File_list = list_all_files(data_file_Path)
File_list.sort()
i = 0
for file in File_list:
#parse the file names given in absolute path
file_name = file.strip().split('/')[-1]
file_name = file_name.rstrip('.txt.tsv')
# function to read tsv files from a given directory
#open the file
data = open(file,'r')
#reading the label line
labelLine = data.readline()
label = labelLine.strip().split('\t')
#write the label
#updates2010July26: update following label writing code
if i == 0:
# output.write('cluster origin')
for element in label:
output.write(element)
output.write('\t')
i = i + 1
#updates2010July26 End
#switch to the next line
output.write('\n')
#read the GO terms
GOTermLines = data.readlines()
for GOTerm in GOTermLines:
GOTerm = GOTerm.strip().strip('\t')
if GOTerm != '':
#updates2010July26: remove the code to write the first column 'REMc cluster ID'
#output.write(file_name)
#output.write('\t')
##updates2010July26 update end
output.write(GOTerm + '\n')
#output.write('\n')
output.close()

55
workflow/templates/qhtcp/REMc/GTF/terms2tsv_v4.pl
View File

@@ -1,55 +0,0 @@
#!/usr/bin/perl
use strict;
use warnings;
use diagnostics;
use File::Map qw(map_file);
my $infile = shift;
my $input;
map_file $input, $infile;
{
local $_ = $input;
(my $f = $infile) =~ s/(.*\/)?(.*)(\.[^\.]*){2}/$2/;
my %orfgene = (/(Y\w+)\s+(\w+)\n/g);
my @indices = (/\Q-- \E(\d+) of \d+\Q --\E/g);
my @ids = (/GOID\s+GO:(\d+)/g);
my @terms = (/TERM\s+(.*?)\n/g);
my @pvalues = (/\nCORRECTED P-VALUE\s+(\d.*?)\n/g);
my @clusterf = (/NUM_ANNOTATIONS\s+(\d+ of \d+)/g);
my @bgfreq = (/, vs (\d+ of \d+) in the genome/g);
my @orfs = (/The genes annotated to this node are:\n(.*?)\n/g);
s/, /:/g for @orfs;
my @genes;
for my $orf (@orfs) {
my @otmp = split /:/, $orf;
my @gtmp = map { $orfgene{$_} } @otmp;
push @genes, (join ':', @gtmp);
}
&header();
for my $i (0 .. (@ids - 1)) {
&report($f, $ids[$i], $terms[$i], $pvalues[$i], $clusterf[$i], $bgfreq[$i], $orfs[$i], $genes[$i]);
}
}
sub header {
print "REMc ID\tGO_term ID\tGO-term\tCluster frequency\tBackground frequency\tP-value\tORFs\tGenes\n";
}
sub report {
my ($f, $id, $term, $p, $cfreq, $bgfreq, $orfs, $genes) = @_;
$cfreq =~ /(\d+) of (\d+)/;
$cfreq = sprintf "%d out of %d genes, %.1f%%", $1, $2, (100*$1/$2);
$bgfreq =~ /(\d+) of (\d+)/;
$bgfreq = sprintf "%d out of %d genes, %.1f%%", $1, $2, (100*$1/$2);
print "$f\t$id\t$term\t$cfreq\t$bgfreq\t$p\t$orfs\t$genes\n";
}

6310
workflow/templates/qhtcp/REMc/GeneByGOAttributeMatrix_nofiltering-2009Dec07.tab
View File

File diff suppressed because one or more lines are too long

264
workflow/templates/qhtcp/REMc/JoinInteractExps3dev.R
View File

@@ -1,264 +0,0 @@
#JoinInteractExps3dev.R
#User prompt for std Value
cat("Enter a Standard Deviation value to filter data to be used by REMc ?\n")
inp <- readLines(file("stdin"), n = 1L)
cat(paste("Standard Deviation Value is", inp, "\n"))
inp<- as.numeric(inp)
#set std deviation multiplier default if no user entry
if(is.double(inp) & !is.na(inp)){
std= inp
}else{std= 2}
#The input files should be entered in order from the greatest number of rows(Orfs) to least.
#Args <- commandArgs(TRUE)
#if(length(Args)==0){
# std=0
#}else{
# std=Args[1]
#}
print(paste("SD=",std))
#Wstudy= getwd()
#dir.create("../JoinFiles") #(paste0(Wstudy,"/JoinFiles"))
#outDir <- "../JoinFiles" #paste0(Wstudy,"/JoinFiles")
#dir.create("../REMc") #(paste0(Wstudy,"/JoinFiles"))
outDir <- "./" #../REMc" #paste0(Wstudy,"/JoinFiles")
print(outDir)
#Args= 'asdf'
#Args[1]= "../Exp1/ZScores/ZScores_Interaction.csv" #paste0(Wstudy,"/Exp1/ZScores/ZScores_Interaction.csv")
#Args[2]= "../Exp2/ZScores/ZScores_Interaction.csv" #paste0(Wstudy,"/Exp2/ZScores/ZScores_Interaction.csv")
#Args[3]= "../Exp3/ZScores/ZScores_Interaction.csv" #paste0(Wstudy,"/Exp1/ZScores/ZScores_Interaction.csv")
#Args[4]= "../Exp4/ZScores/ZScores_Interaction.csv" #paste0(Wstudy,"/Exp2/ZScores/ZScores_Interaction.csv")
#Args[3]= "../Exp3/ZScores/ZScores_Interaction.csv" #paste0(Wstudy,"/Exp3/ZScores/ZScores_Interaction.csv")
inputFile= 'asdf'
if (file.exists('../Exp1/ZScores/ZScores_Interaction.csv')){
inputFile <- '../Exp1/ZScores/ZScores_Interaction.csv'
}
if (file.exists('../Exp2/ZScores/ZScores_Interaction.csv')){
inputFile[2] <- '../Exp2/ZScores/ZScores_Interaction.csv'
}
if (file.exists('../Exp3/ZScores/ZScores_Interaction.csv')){
inputFile[3] <- '../Exp3/ZScores/ZScores_Interaction.csv'
}
if (file.exists('../Exp4/ZScores/ZScores_Interaction.csv')){
inputFile[4] <- '../Exp4/ZScores/ZScores_Interaction.csv'
}
#Args= inputFile
#outDir <- ArgsJoin[1] #Output Directory
print(length(inputFile)) #display the number of arguments on terminal
#open required library for the join function (libraries must already be install on R)
library(plyr)
library(sos)
library(dplyr)
#read in the files for your experiment and
#join the two files at a time as a function of how many inputFile, list the larger file first ? in this example X2 has the larger number of genes.
#if X1 has a larger number of genes, switch the order of X1 and X2
if(length(inputFile)==2) {
X1 <- read.csv(file= inputFile[1],stringsAsFactors = FALSE)
X2 <- read.csv(file= inputFile[2],stringsAsFactors = FALSE)
X <- join(X1,X2,by="OrfRep")
OBH=X[,order(colnames(X))] #OrderByHeader
headSel= select(OBH, contains('OrfRep'), matches('Gene'), contains('Z_lm_K'), contains('Z_Shift_K'),contains('Z_lm_L'), contains('Z_Shift_L'))
headSel= select(headSel, -'Gene.1') #remove 'Gene.1 column
headSel2 = select(OBH, contains('OrfRep'), matches('Gene')) #Frame for interleaving Z_lm with Shift colums
headSel2 = select(headSel2, -'Gene.1') #remove 'Gene.1 column #Frame for interleaving Z_lm with Shift colums
}else if(length(inputFile)==3){
X1 <- read.csv(file= inputFile[1],stringsAsFactors = FALSE) #exp1File,stringsAsFactors = FALSE)
X2 <- read.csv(file= inputFile[2],stringsAsFactors = FALSE) #exp2File,stringsAsFactors = FALSE)
X3 <- read.csv(file= inputFile[3],stringsAsFactors = FALSE) #exp3File,stringsAsFactors = FALSE)
X <- join(X1,X2,by="OrfRep")
X <- join(X,X3,by="OrfRep")
OBH=X[,order(colnames(X))] #OrderByHeader
headSel= select(OBH, contains('OrfRep'), matches('Gene'), contains('Z_lm_K'), contains('Z_Shift_K'),contains('Z_lm_L'), contains('Z_Shift_L'))
headSel= select(headSel, -'Gene.1',-'Gene.2')
headSel2 = select(OBH, contains('OrfRep'), matches('Gene'))
headSel2 = select(headSel2, -'Gene.1',-'Gene.2')
}else if(length(inputFile)==4){
X1 <- read.csv(file= inputFile[1],stringsAsFactors = FALSE) #exp1File,stringsAsFactors = FALSE)
X2 <- read.csv(file= inputFile[2],stringsAsFactors = FALSE) #exp2File,stringsAsFactors = FALSE)
X3 <- read.csv(file= inputFile[3],stringsAsFactors = FALSE) #exp3File,stringsAsFactors = FALSE)
X4 <- read.csv(file= inputFile[4],stringsAsFactors = FALSE) #exp4File,stringsAsFactors = FALSE)
X <- join(X1,X2,by="OrfRep")
X <- join(X,X3,by="OrfRep")
X <- join(X,X4,by="OrfRep")
OBH=X[,order(colnames(X))] #OrderByHeader
headSel= select(OBH, contains('OrfRep'), matches('Gene'), contains('Z_lm_K'), contains('Z_Shift_K'),contains('Z_lm_L'), contains('Z_Shift_L'))
headSel= select(headSel, -'Gene.1',-'Gene.2',-'Gene.3')
headSel2 = select(OBH, contains('OrfRep'), matches('Gene'))
headSel2 = select(headSel2, -'Gene.1',-'Gene.2',-'Gene.3')
}
#headSel$contains('Z_Shift') %>% replace_na(0.001)
headers<-colnames(headSel)
i=0
for(i in 1:length(headers)){
if(grepl("Shift",headers[i])) {
headSel[headers[i]][is.na(headSel[headers[i]])] = 0.001
}
if(grepl("Z_lm_",headers[i])) {
headSel[headers[i]][is.na(headSel[headers[i]])] = 0.0001
}
}
#2SD option code to exclude Z_lm values less than 2 standard Deviations
#+++++++++++++
REMcRdy= select(headSel, contains('OrfRep'), matches('Gene'), contains('Z_lm_'))
shiftOnly= select(headSel, contains('OrfRep'), matches('Gene'), contains('Z_Shift'))
#Code to replace the numeric (.1 .2 .3) headers with experiment names from StudyInfo.txt
Labels <- read.csv(file= "../Code/StudyInfo.csv",stringsAsFactors = FALSE,sep= ",")
#Using Text search grepl to relabel headers+++++++++++++++++++++++++++++++++++++++++
REMcRdyHdr= colnames(REMcRdy)
REMcRdyLabels= 'asdf'
shftHdr= colnames(shiftOnly)
shiftLabels='asdf'
shiftLabels[1:2]<-shftHdr[1:2]
REMcRdyLabels[1:2]<-REMcRdyHdr[1:2]
for(i in 3:(length(shftHdr))){
if(i==3){
shiftLabels[3]<-paste0(Labels[1,2],".",shftHdr[3])
REMcRdyLabels[3]<-paste0(Labels[1,2],".",REMcRdyHdr[3]) }
if(i==5){
shiftLabels[5]<-paste0(Labels[1,2],".",shftHdr[5])
REMcRdyLabels[5]<-paste0(Labels[1,2],".",REMcRdyHdr[5])
}
if(i==7){
shiftLabels[7]<-paste0(Labels[1,2],".",shftHdr[7])
REMcRdyLabels[7]<-paste0(Labels[1,2],".",REMcRdyHdr[7])
}
if(grepl(".1",shftHdr[i],fixed=true)){
shiftLabels[i]<-paste0(Labels[2,2],".",shftHdr[i])
REMcRdyLabels[i]<-paste0(Labels[2,2],".",REMcRdyHdr[i])}
if (grepl(".2",shftHdr[i],fixed=true)){
shiftLabels[i]<-paste0(Labels[3,2],".",shftHdr[i])
REMcRdyLabels[i]<-paste0(Labels[3,2],".",REMcRdyHdr[i])}
if(grepl(".3",shftHdr[i],fixed=true)){
shiftLabels[i]<-paste0(Labels[4,2],".",shftHdr[i])
REMcRdyLabels[i]<-paste0(Labels[4,2],".",REMcRdyHdr[i])}
}
for(i in 3:(length(REMcRdyLabels))){
j=as.integer(i)
REMcRdyLabels[j]<- gsub("[.]", "_", REMcRdyLabels[j])
shiftLabels[j]<- gsub("[.]", "_", shiftLabels[j])
}
colnames(shiftOnly)<- shiftLabels
colnames(REMcRdy)<- REMcRdyLabels
#+++++++++++++++++++++++
combI= headSel2 #Starting Template orf, Genename columns
#headersRemc<-colnames(REMcRdy)
#Reoder columns to produce an interleaved set of Z_lm and Shift data for all the cpps.
for(i in 3:length(colnames(REMcRdy))){
combI=cbind.data.frame(combI, shiftOnly[i])
combI=cbind.data.frame(combI, REMcRdy[i])
}
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Vec1= NA
Vec2= NA
Vec3= NA
Vec4= NA
Vec5= NA
Vec6= NA
Vec7= NA
Vec8= NA
if(length(REMcRdy)== 6){
Vec1= abs(REMcRdy[,3])>=std
Vec2= abs(REMcRdy[,4])>=std
Vec3= abs(REMcRdy[,5])>=std
Vec4= abs(REMcRdy[,6])>=std
bolVec= Vec1 | Vec2 |Vec3 |Vec4
REMcRdyGT2= REMcRdy[bolVec,1:2]
REMcRdyGT2[ ,3:6]= REMcRdy[bolVec,3:6]
shiftOnlyGT2= shiftOnly[bolVec,1:2]
shiftOnlyGT2[ ,3:6]= shiftOnly[bolVec,3:6]
}
if(length(REMcRdy)== 8){
Vec1= abs(REMcRdy[,3])>=std
Vec2= abs(REMcRdy[,4])>=std
Vec3= abs(REMcRdy[,5])>=std
Vec4= abs(REMcRdy[,6])>=std
Vec5= abs(REMcRdy[,7])>=std
Vec6= abs(REMcRdy[,8])>=std
bolVec= Vec1 | Vec2 |Vec3 | Vec4 |Vec5 |Vec6
REMcRdyGT2= REMcRdy[bolVec,1:2]
REMcRdyGT2[ ,3:8]= REMcRdy[bolVec,3:8]
shiftOnlyGT2= shiftOnly[bolVec,1:2]
shiftOnlyGT2[ ,3:8]= shiftOnly[bolVec,3:8]
}
if(length(REMcRdy)== 10){
Vec1= abs(REMcRdy[,3])>=std
Vec2= abs(REMcRdy[,4])>=std
Vec3= abs(REMcRdy[,5])>=std
Vec4= abs(REMcRdy[,6])>=std
Vec5= abs(REMcRdy[,7])>=std
Vec6= abs(REMcRdy[,8])>=std
Vec7= abs(REMcRdy[,9])>=std
Vec8= abs(REMcRdy[,10])>=std
bolVec= Vec1 | Vec2 |Vec3 |Vec4|Vec5|Vec6|Vec7|Vec8
REMcRdyGT2= REMcRdy[bolVec,1:2]
REMcRdyGT2[ ,3:10]= REMcRdy[bolVec,3:10]
shiftOnlyGT2= shiftOnly[bolVec,1:2]
shiftOnlyGT2[ ,3:10]= shiftOnly[bolVec,3:10]
}
if(std!=0){
REMcRdy= REMcRdyGT2 #[,2:length(REMcRdyGT2)]
shiftOnly= shiftOnlyGT2 #[,2:length(shiftOnlyGT2)]
}
if(std==0){
REMcRdy= REMcRdy #[,2:length(REMcRdy)]
shiftOnly= shiftOnly #[,2:length(shiftOnly)]
}
#R places hidden "" around the header names. The following
# is intended to remove those quote so that the "" donot blow up the Java REMc.
#Use ,quote=F in the write.csv statement to fix R output file.
print(paste("SD=",std))
print(getwd())
#write.csv(combI,file = file.path(outDir,"CombinedKLzscores.csv"),row.names = FALSE)
write.csv(REMcRdy,file = file.path(outDir,"REMcRdy_lm_only.csv"),row.names = FALSE, quote=F)
write.csv(shiftOnly,file = file.path(outDir,"Shift_only.csv"),row.names = FALSE, quote=F)
#LabelStd <- read.table(file= "./Parameters.csv",stringsAsFactors = FALSE,sep= ",")
pwd=getwd()
print(getwd)
LabelStd<- read.csv(file= "../Code/StudyInfo.csv",stringsAsFactors = FALSE)
print(std)
LabelStd[,4]= as.numeric(std)
write.csv(LabelStd,file="../Code/Parameters.csv",row.names = FALSE)
write.csv(LabelStd,file="../Code/StudyInfo.csv",row.names = FALSE)
cat(paste("Standard Deviation Value was set as", std, "\n"))

5150
workflow/templates/qhtcp/REMc/ORF_List_Without_DAmPs.txt
View File

File diff suppressed because it is too large Load Diff

BIN
workflow/templates/qhtcp/REMc/jingyuJava_1_7_extractLib.jar
View File

Binary file not shown.