diff --git a/workflow/script-run-workflow b/workflow/script-run-workflow index 2c2ec5ed..2bf4013b 100755 --- a/workflow/script-run-workflow +++ b/workflow/script-run-workflow @@ -247,7 +247,8 @@ install_dependencies() { depends_brew=(graphiz pandoc gd pdftk-java shdoc nano rsync) depends_perl=(File::Map ExtUtils::PkgConfig GD GO::TermFinder) depends_r=(BiocManager ontologyIndex ggrepel tidyverse sos openxlsx ggplot2 - plyr extrafont gridExtra gplots stringr plotly ggthemes pandoc rmarkdown) + plyr extrafont gridExtra gplots stringr plotly ggthemes pandoc rmarkdown + plotly htmlwidgets) depends_bioc=(org.Sc.sgd.db) [[ $1 == "--get-depends" ]] && return 0 # if we just want to read the depends vars @@ -341,6 +342,10 @@ init_project() { module easy # @section EASY # @description Start an EASY analysis +# Eliminated EstartConsole.m +# TODO Don't create output in the scans folder, put it in an output directory +# TODO The !!Results output files need standardized naming +# TODO Don't perform directory operations in EASY # * The QHTCPImageFolders and 'MasterPlateFiles' folder are the inputs for image analysis with EASY software. # * EASY will automatically generate a 'Results' directory (within the ExpJobs/'ExperimentJob' folder) w/ timestamp and an optional short description provided by the user (Fig.2). # * The 'Results' directory is created and entered, using the "File >> New Experiment" dropdown in EASY. @@ -400,9 +405,34 @@ module easy # * When finished, the '!!ResultsStd_.txt' will be about the same file size and it should be used in the following StudiesQHTCP analysis. # * 'NoGrowth_.txt', and 'GrowthOnly_.txt' files will be generated in the 'PrintResults' folder. # +# Issues: +# * We need full documentation for all of the current workflow. There are different documents that need to be integrated. This will need to be updated as we make improvements to the system. +# * MasterPlate_ file must have ydl227c in orf column, or else it Z_interaction.R will fail, because it can’t calculate shift values. +# * Make sure there are no special characters; e.g., (), “, ‘, ?, etc.; dash and underscore are ok as delimiters +# * Drug_Media_ file must have letter character to be read as ‘text’. +# * MasterPlate_ file and DrugMedia_ are .xlsx or .xls, but !!Results_ is .txt. +# * In Z_interactions.R, does it require a zero concentration/perturbation (should we use zero for the low conc, even if it’s not zero), e.g., in order to do the shift correctly. +# * Need to enable all file types (not only .xls) as the default for GenerateResults (to select MP and DM files as .xlsx). +# * Explore differences between the ELR and STD files - 24_0414; John R modified Z script to format ELR file for Z_interactions.R analysis. +# * To keep time stamps when transferring with FileZilla, go to the transfer drop down and turn it on, see https://filezillapro.com/docs/v3/advanced/preserve-timestamps/ +# * Could we change the ‘MasterPlateFiles’ folder label in EASY to ‘MasterPlate_DrugMedia’ (since there should be only one MP and there is also a DM file required? +# * I was also thinking of adding a ‘MasterPlateFilesOnly’ folder to the QHTCP directory template where one could house different MPFiles (e.g., with and without damps, with and without Refs on all MPs, etc; other custom MPFiles, updated versions, etc) +# * Currently updated files are in ‘23_1011_NewUpdatedMasterPlate_Files’ on Mac (yeast strains/23_0914…/) +# * For EASY to report cell array positions (plate_row_column) to facilitate analyzing plate artifacts. The MP File in Col 3 is called ‘LibraryLocation’ and is reported after ‘Specifics’ in the !!Results. +# * Can EASY/StudiesQ-HTCP be updated at any time by rerunning with updated MP file (new information for gene, desc, etc)- or maybe better to always start with a new template? +# * Need to be aware of file formatting to avoid dates (e.g., with gene names like MAY24, OCT1, etc, and with plate locations 1E1, 1E2, etc)- this has been less of a problem. +# * In StudiesQHTCP folders, remember to annotate Exp1, Exp2, in the StudyInfo.csv file. +# * Where are gene names called from for labeling REMc heatmaps, TSHeatmaps, Z-interaction graphs, etc? Is this file in the QHTCP ‘code’ folder, or is it in the the results file (and thus ultimately the MP file)? +# * Is it ok for a MasterPlate_ file to have multiple sheets (e.g., readme tab- is only the first tab read in)? +# * What are the rules for pulling information from the MasterPlateFile to the !!Results_ (e.g., is it the column or the Header Name, etc that is searched? Particular cells in the DrugMedia file?). +# * Modifier, Conc are from DM sheet, and refer to the agar media arrays. OrfRep is from MasterPlate_ File. ‘Specifics’ (Last Column) is experiment specific and accommodate designs involving differences across the multi-well liquid arrays. ‘StrainBkGrd’ (now ‘Library location’) is in the 3rd column and reported after ‘Specifics’ at the last col of the ‘!!Results..’ file. +# * Do we have / could we make an indicator- work in progress or idle/complete with MP/DM and after gen-report. Now, we can check for the MPDMmat.mat file, or we can look in PrintResults, but would be nice to know without looking there. +# * File>>Load Experiment wasn’t working (no popup to redirect). Check this again. + easy() { debug "Running: ${FUNCNAME[0]}" - EASY="/mnt/data/EASY/EasyDev2024/BU/EASY240430AppExported/EstartConsole.m" + #EASY="/mnt/data/EASY/EasyDev2024/BU/EASY240430AppExported/EstartConsole.m" + EASY="/mnt/data/EASY/EasyDev2024/BU/EASY240505AppExported/EASYConsole.m" pushd "$SCAN_DIR" || return 1 @@ -447,6 +477,8 @@ easy() { ! ((YES)) && ask "Start EASY in MATLAB? This requires a GUI." && matlab -nosplash -r "$EASY" # glob EASY output and make sure it exists + # currently this is just for informative purposes of how to glob some of the EASY output files + # The EASY output files need to be standardized shopt -s nullglob EASY_RESULTS_DIRS=( Results* ) shopt -u nullglob @@ -508,7 +540,11 @@ module qhtcp # * Do not double-click on the file from the directory. # * When prompted, navigate to the ExpJobs folder and the PrintResults folder within the correct job folder. # * Repeat this for every Exp# folder depending on how many experiments are being performed. +# * Note: Before doing this, it’s a good idea to compare the ref and non-ref CPP average and median values. If they are not approximately equal, then may be helpful to standardize Ref values to the measures of central tendency of the Non-refs, because the Ref CPPs are used for the z-scores, which should be centered around zero. +# * This script will copy the !!ResultsStd file (located in /PrintResults in the relevant job folder in /ExpJobs **rename this !!Results file before running front end; we normally use the ‘STD’ (not the ‘ELR’ file) chosen to the Exp# directory as can be seen in the “Current Folder” column in MATLAB, and it updates ‘StudiesDataArchive.txt’ file that resides in the /StudiesQHTCP folder. ‘StudiesDataArchive.txt’ is a log of file paths used for different studies, including timestamps. # +# Do this to document the names, dates and paths of all the studies and experiment data used in each study. Note, one should only have a single ‘!!Results…’ file for each /Exp_ to prevent ambiguity and confusion. If you decide to use a new or different ‘!!Results…’ sheet from what was used in a previous “QHTCP Study”, remove the one not being used. NOTE: if you copy a ‘!!Results…’ file in by hand, it will not be recorded in the ‘StudiesDataArchive.txt’ file and so will not be documented for future reference. If you use the ExpFrontend.m utility it will append the new source for the raw !!Results… to the ‘StudiesDataArchive.txt’ file. +# As stated above, it is advantageous to think about the comparisons one wishes to make so as to order the experiments in a rational way as it relates to the presentation of plots. That is, which results from sheets and selected ‘interaction … .R’, user modified script, is used in /Exp1, Exp2, Exp3 and Exp4 as explained in the following section. # TODO MUST CLEAN UP QHTCP TEMPLATE DIRECTORY # # Code/Directory Structure: @@ -827,6 +863,7 @@ qhtcp() { QHTCP_TEMPLATE_DIR="$SCRIPT_DIR/templates/qhtcp" STUDY_TEMPLATE_DIR="$QHTCP_TEMPLATE_DIR/ExpTemplate" OUTPUT_DIR="/mnt/data/StudiesQHTCP" + STUDIES_ARCHIVE="$OUTPUT_DIR/StudiesDataArchive.txt" QHTCP_DIR="$OUTPUT_DIR/$PROJECT" if [[ -d $QHTCP_DIR ]]; then @@ -843,23 +880,22 @@ qhtcp() { fi # Print current studies - STUDY_INFO_FILE="$QHTCP_DIR/Code/StudyInfo.csv" - [[ -f $STUDY_INFO_FILE ]] && - echo "Current studies from $STUDY_INFO_FILE" && - cat "$STUDY_INFO_FILE" + STUDY_INFO="$QHTCP_DIR/Code/StudyInfo.csv" + [[ -f $STUDY_INFO ]] && + echo "Current studies from $STUDY_INFO" && + cat "$STUDY_INFO" - # Ask user edit STUDY_INFO_FILE - if ! ((YES)) && ask "Would you like to edit $STUDY_INFO_FILE to add or modify studies?"; then + # Ask user to edit STUDY_INFO + if ! ((YES)) && ask "Would you like to edit $STUDY_INFO to add or modify studies?"; then cat <<-EOF Give each experiment the labels you wish to be used for the plots and specific files. Enter the desired Experiment names and order them in the way you want them to appear in the REMc heatmaps EOF - nano "$STUDY_INFO_FILE" + nano "$STUDY_INFO" fi - # Sets STUDIES_NUM - get_studies "$STUDY_INFO_FILE" + get_studies "$STUDY_INFO" # Initialize missing dirs for s in "${STUDIES_NUM[@]}"; do @@ -872,43 +908,26 @@ qhtcp() { fi done - # MATLAB stuff - cat <<-EOF - ExpFrontend.m was made for recording into a spreadsheet - ('StudiesDataArchive.txt') the date and files used (i.e., directory paths to the - !!Results files used as input for Z-interaction script) for each multi-experiment study. - - Run the front end MATLAB programs in the correct order (e.g., run front end in "exp1" - folder to call the !!Results file for the experiment you named as exp1 in the StudyInfo.csv file) - The GTA and pairwise, TSHeatmaps, JoinInteractions and GTF Heatmap scripts use this table - to label results and heatmaps in a meaningful way for the user and others. - The BackgroundSD and ZscoreJoinSD fields will be filled automatically according to user - specifications, at a later step in the QHTCP study process. - - Open MATLAB and in the application navigate to each specific /Exp folder, - call and execute ExpFrontend.m by clicking the play icon. - Use the "Open file" function from within Matlab. - Do not double-click on the file from the directory. - When prompted, navigate to the ExpJobs folder and the PrintResults folder within the correct job folder. - Repeat this for every Exp# folder depending on how many experiments are being performed. - The Exp# folder must correspond to the StudyInfo.csv created above. - EOF - if ! ((YES)) && - ask "Start MATLAB to run ExpFrontend.m? This requires a GUI."; then - matlab -nosplash - fi + mat_exp_frontend - # Enter REMc directory to run the scripts there - pushd "$QHTCP_DIR/REMc" || return 1 - # Run modules - r_join_interact - java_extract - r_add_shift_values - r_heat_maps_zscores - r_heat_maps_homology - popd || return 1 + for s in "${STUDIES_NUM[@]}"; do + study_dir="$QHTCP_DIR/Exp$s" + # Z_InteractionTemplate.R + r_interactions "$study_dir" "!!Results"* + done + + # Enter REMc directory to run the scripts there + pushd "$QHTCP_DIR/REMc" || return 1 + + # Run modules + r_join_interact + java_extract + r_add_shift_values + r_heat_maps_zscores + r_heat_maps_homology + popd || return 1 } @@ -961,6 +980,95 @@ gta() { # * Functions you do not want to perform by default (submodules should be called modules) # * Should not call cd or pushd (let module dictate) +submodule mat_exp_frontend +# @description Run the ExpFrontend.m program +mat_exp_frontend() { + debug "Running: ${FUNCNAME[0]}" + # MATLAB stuff + cat <<-EOF + ExpFrontend.m was made for recording into a spreadsheet + ('StudiesDataArchive.txt') the date and files used (i.e., directory paths to the + !!Results files used as input for Z-interaction script) for each multi-experiment study. + + Run the front end MATLAB programs in the correct order (e.g., run front end in "exp1" + folder to call the !!Results file for the experiment you named as exp1 in the StudyInfo.csv file) + The GTA and pairwise, TSHeatmaps, JoinInteractions and GTF Heatmap scripts use this table + to label results and heatmaps in a meaningful way for the user and others. + The BackgroundSD and ZscoreJoinSD fields will be filled automatically according to user + specifications, at a later step in the QHTCP study process. + + Open MATLAB and in the application navigate to each specific /Exp folder, + call and execute ExpFrontend.m by clicking the play icon. + Use the "Open file" function from within Matlab. + Do not double-click on the file from the directory. + When prompted, navigate to the ExpJobs folder and the PrintResults folder within the correct job folder. + Repeat this for every Exp# folder depending on how many experiments are being performed. + The Exp# folder must correspond to the StudyInfo.csv created above. + EOF + if ! ((YES)) && + ask "Start MATLAB to run ExpFrontend.m? This requires a GUI."; then + matlab -nosplash + fi + [[ -f $STUDIES_ARCHIVE ]] || (err "$STUDIES_ARCHIVE missing"; return 1) +} + + +submodule r_interactions +# @description Run the R interactions analysis +# TODO don't want to rename Z_InteractionTemplate.R because that will break logic, just edit in place instead +# Is this script interactive +# @arg $1 string The current working directory +r_interactions() { + debug "Running: ${FUNCNAME[0]}" + + cat <<-EOF +In each /Exp# folder, rename the Z_InteractionTemplate.R script according to the experiment focus +Example: Interaction, Experimenter Initials, Experiment Focus --> ‘int_RM_2PE.R’ +5. Open the renamed interaction script, and edit each one beginning at the ‘++BEGIN USER DATA SELECTION++’ +This is designed so that the data of interest for each experiment is appropriately selected from the !!Results…txt file +The user can edit, step through, and test the R script without running through the whole routine by observing the resultant data table created in RStudio. +The Z_InteractionTemplate.R script has a collection of code lines that have been used for prior analyses (generally to select data from various !!Results…txt files), which may be commented out (if not relevant), reused as needed, and/or modified for a new study. These include lines associated with the removal of ‘dAmps’, specific concentrations, and items described in the ‘Specifics’ and ‘Media’, i.e., information specific to a particular experiment design. There are also code lines to replace gene names ‘OCT1/YKL134C’ /’MAY24/YPR153W’ and that get converted to date format in excel, by using only the ORF name and to remove data rows with ‘Blank’ listed; these lines of code convenient to reuse. Hopefully, these code lines can be used, commented out, or adapted to aid the user in modifying this section to the specific data requirements of the study. As a new user data filter code is developed for each ‘Study’ (and vetted), those lines can be added to the InteractionTemplate230119.R code in the /StudyTemplate folders to aid in future studies. +6. Open a terminal, navigate to each /Exp# folder, and execute the (customized) ‘Z_InteractionTemplate_…” script by using the command line below: + + +Rscript RenamedInteractionTemplate.R \!\!Results… .txt + +**need to change wording to choose SD of Delta_Background to exclude Data from analysis. +[1] "Be sure to enter Background noise filter standard deviation i.e., 3 or 5 per Sean" +Enter a Standard Deviation value to noise filter >> + +[1] Enter Standard deviation value for removing data for cultures due to high background (e.g., contaminated cultures). Generally set this very high (e.g., ‘20’) on the first run in order NOT to remove data, e.g. ‘20’. Review QC data and inspect raw image data to decide if it is desirable to remove data, and then rerun analysis. +Enter a Background SD threshold for EXCLUDING culture data from further analysis: + + +The script will request for the user to input a ‘Background Standard Deviation Value’. This Background value removes data where there is high pixel intensity in the background regions of a spot culture (i.e., suspected contamination). 5 is a minimum recommended value, because lower values result in more data being removed, and often times this is undesirable if contamination occurs late after the carrying capacity of the yeast culture is reached. This is most often “trial and error”, meaning there is a ‘Frequency_Delta_Background.pdf’ report in the /Exp_/ZScores/QC/ folder to evaluate whether the chosen value was suitable (and if not the analysis can simply be rerun with a more optimal choice). In general, err on the high side, with BSD of 10 or 12…. One can also use EZview to examine the raw images and individual cultures potentially included/excluded as a consequence of the selected value. Background values are reported in the results sheet and so could also be analyzed there.. + + (For new terminal users, directory navigation tips are described below) +To navigate to the directory one can use the directory GUI (in X2Go, use the GUI to navigate to desired operating directory and then from the ‘File’ menu, choose “Open in Terminal’) +Alternatively, navigate there through the terminal window: ‘pwd’ “prints the current working directory”, ‘ls’ “lists” the subfolders in the current directory. ‘cd’’ followed by the name of the ‘subdirectory’ will move down into it. “cd .. “ changes to the parent directory +The tab key can be used to autofill unique characters after typing the initial letters of a folder or file you wish to call. + + +The template structure above assists the user with organization and management of Q-HTCP files and provides a uniform directory structure to streamline reference across different users and experiments. +Since we are systematically comparing perturbations, most Q-HTCP studies will consist of either 2 or 4 experiment subfolders. +The Zscores files are used for subsequent analyses, including REMc, GTA and Term Specific Heatmaps. These further analyses are described below and can be completed in any order and/or concurrently from separate terminals. +**Annotate Files produced and comment out code that produces files that are obsolete or clutter. + + EOF + + script="$1/Z_InteractionTemplate.R" + + debug "$RSCRIPT $script" + "$RSCRIPT" "$script" + + # 1. Path to input file +# 2. /output/ directory +# 3. Path to StudyInfo.csv +# 4. Standard deviation value + +} + + submodule r_join_interact # @description JoinInteractExps3dev.R creates REMcRdy_lm_only.csv and Shift_only.csv r_join_interact() { diff --git a/workflow/templates/easy/DMPexcel2mat.m b/workflow/templates/easy/DMPexcel2mat.m new file mode 100755 index 00000000..716d9391 --- /dev/null +++ b/workflow/templates/easy/DMPexcel2mat.m @@ -0,0 +1,246 @@ +%% CALLED BY EASYconsole.m %% +w = pwd; +%{ +cd .. +if exist('ExpJobs/') + cd('ExpJobs/'); +else + pwd +end +%} + +try + cd(ExpPath) +catch + cd(w) +end +pwd + +numOfMPs=0; +%GUI input for selecting a MasterPlate Excel file +questdlg('\fontsize{20} Select Master Plate File','File Selection','OK', struct('Default','OK','Interpreter','tex')) +[Scanfiles, pathname]=uigetfile('*.*','MultiSelect','off'); +if ispc + MPdir= fullfile(pathname,'\'); +else + MPdir=fullfile(pathname,'/'); +end + + +infile= Scanfiles(1,:); +cd(MPdir) +%fid=fopen(infile)%('exp23PrintTimes.xls'); % textread puts date and time sequentially into vector +if ispc + [num, txt, raw] = xlsread(infile); %,'Yor1HitsMPsetFinal'); + fields= {txt(2,1:15)}; %or 1:17 for later but dont wish to exceed and cause error ? if used + +else + clear MPtbl + opts = detectImportOptions(infile); + MPtbl = readtable(infile,opts); + MPtbl= readtable(infile); + fields= {opts.VariableNames}; %? if used anywhere although 'saved' to MPDMmat + MPcell= readcell(infile); +end + + +cd(w) +numb= 0; +clear MP; +try + + if ispc + excLnNum=3; + while (isequal(txt{excLnNum,1},'###')) +numb=numb+1; + MP(numb).head = {raw(excLnNum,2:6)}; + MP(numb).recNum = {raw((excLnNum+1):(excLnNum+384),1)}; + MP(numb).orf = {raw((excLnNum+1):(excLnNum+384),2)}; + MP(numb).strain = {raw((excLnNum+1):(excLnNum+384),3)}; + MP(numb).genename= {raw((excLnNum+1):(excLnNum+384),12)}; + MP(numb).drug= {raw((excLnNum+1):(excLnNum+384),8)}; + MP(numb).media= {raw((excLnNum+1):(excLnNum+384),7)}; +if size(raw,2)>15 + MP(numb).orfRep= {raw((excLnNum+1):(excLnNum+384),16)}; %added 12_1005 to specify replicates Orfs in MP + MP(numb).specifics= {raw((excLnNum+1):(excLnNum+384),17)}; %added 12_1008 to specify replicates Specific details in MP +else + MP(numb).orfRep= ' '; + MP(numb).specifics= ' '; +end +%{ +if size(raw,2)>17 ; Future MP field + MP(numb).specifics2= {raw((excLnNum+1):(excLnNum+384),18)}; %added 12_1008 to specify strain Bkground in MP +else + MP(numb).specifics2= ' '; +end +%} +excLnNum=excLnNum+385; + +msg=strcat('NumberOfMP = ',num2str(numb), ' lastLineNo. = ',num2str(excLnNum)); + + end % end if ispc + else %if not ispc + excLnNum=1; + while (isequal(MPcell{(excLnNum+2),1},'###')) + numb=numb+1; + MP(numb).head = {MPtbl(excLnNum,2:6)}; + MP(numb).head{1}= table2cell(MP(numb).head{1}); + MP(numb).recNum = {MPtbl((excLnNum+1):(excLnNum+384),1)}; + MP(numb).recNum{1}= table2cell(MP(numb).recNum{1}); + MP(numb).orf = {MPtbl((excLnNum+1):(excLnNum+384),2)}; + MP(numb).orf{1}= table2cell(MP(numb).orf{1}); + MP(numb).strain = {MPtbl((excLnNum+1):(excLnNum+384),3)}; + MP(numb).strain{1}= table2cell(MP(numb).strain{1}); + MP(numb).genename= {MPtbl((excLnNum+1):(excLnNum+384),12)}; + MP(numb).genename{1}= table2cell(MP(numb).genename{1}); + MP(numb).media= {MPtbl((excLnNum+1):(excLnNum+384),7)}; + MP(numb).media{1}= table2cell(MP(numb).media{1}); + if size(MPtbl,2)>15 + MP(numb).orfRep= {MPtbl((excLnNum+1):(excLnNum+384),16)}; %added 12_1005 to specify replicates Orfs in MP + MP(numb).orfRep{1}= table2cell(MP(numb).orfRep{1}); + MP(numb).specifics= {MPtbl((excLnNum+1):(excLnNum+384),17)}; %added 12_1008 to specify replicates Specific details in MP + MP(numb).specifics{1}= table2cell(MP(numb).specifics{1}); + else + MP(numb).orfRep= ' '; + MP(numb).specifics= ' '; + end + +excLnNum=excLnNum+385; +msg=strcat('NumberOfMP = ',num2str(numb), 'lastLineNo. = ',num2str(excLnNum)) + end + + end %end of if ispc for MP while loop choice PC or Linux-other (ln~48) + +catch ME + h = msgbox(msg,'Check Number of Master Plates and Excel Lines') + uiwait(h); +end %end for try MP excel sheet input + +%********************************************************************* +%DMupload +%Drug and Media Plate setup Upload from Excel +%scan().plate().LocInd().gene +cd(MPdir); +excLnNum=1; +%w=pwd +%cd('D:\EASY\Experiments\'); %cd(ExpPath); % +numOfDrugs=0; +numOfMedias=0; +%GUI input for selecting a MasterPlate Excel file +questdlg('\fontsize{20} Select DrugMedia File','File Selection','OK', struct('Default','OK','Interpreter','tex')); +[Scanfiles, pathname]=uigetfile('*.*', 'MultiSelect','off'); +if ispc + DMdir= fullfile(pathname,'\'); +else + DMdir=fullfile(pathname,'/'); +end + +clear infile; +infile= Scanfiles(1,:); +cd(DMdir) +%fid=fopen(infile)%('exp23PrintTimes.xls'); % textread puts date and time sequentially into vector +if ispc + [num, txt, raw] = xlsread(infile); %,'Yor1HitsMPsetFinal'); + fields= {txt(2,1:5)}; + Linked= num(1,1); +else + opts = detectImportOptions(infile); + DMtbl= readtable(infile,opts); + fields= opts.VariableOptions; + Linked= DMtbl{1,1}; + DMcell= readcell(infile); +end +cd(w) +numb=0; + +%fields= {txt(2,1:5)}; +%Linked= num(1,1); + + +if isequal(Linked,1) %Drugs and Media are linked 1 to 1; else they are combinatorial +clear DM; +%try +excLnNum=2; +if ispc +while (~isequal(txt{excLnNum,2},'###')) + numb=numb+1; + DM.drug(numb) = {raw(excLnNum,2)}; + DM.conc(numb) = {raw(excLnNum,3)}; + DM.media(numb) = {raw(excLnNum,4)}; + DM.mod1(numb) = {raw(excLnNum,5)}; + DM.conc1(numb) = {raw(excLnNum,6)}; + DM.mod2(numb) = {raw(excLnNum,7)}; + DM.conc2(numb) = {raw(excLnNum,8)}; + excLnNum=excLnNum+1; + msg=strcat('NumberOf1:1DrugMediaPlates = ',num2str(numb), ' lastLineNo. = ',num2str(excLnNum)) +end +else + clear DM + excLnNum=1; + while (~isequal(DMcell{excLnNum+1,2},'###')) + numb=numb+1; + DM.drug(numb) = {DMtbl(excLnNum,2)}; + DM.drug(numb)= table2cell(DM.drug{numb}); + DM.conc(numb) = {DMtbl(excLnNum,3)}; + DM.conc(numb)= table2cell(DM.conc{numb}); + DM.media(numb) = {DMtbl(excLnNum,4)}; + DM.media(numb)= table2cell(DM.media{numb}); + DM.mod1(numb) = {DMtbl(excLnNum,5)}; + DM.mod1(numb)= table2cell(DM.mod1{numb}); + DM.conc1(numb) = {DMtbl(excLnNum,6)}; + DM.conc1(numb)= table2cell(DM.conc1{numb}); + DM.mod2(numb) = {DMtbl(excLnNum,7)}; + DM.mod2(numb)= table2cell(DM.mod2{numb}); + DM.conc2(numb) = {DMtbl(excLnNum,8)}; + DM.conc2(numb)= table2cell(DM.conc2{numb}); + excLnNum=excLnNum+1; + msg=strcat('NumberOf1:1DrugMediaPlates = ',num2str(numb), ' lastLineNo. = ',num2str(excLnNum)) + end +end + + + +end +%Legacy contengency -Not ever used!! +if isequal(Linked,0) %0 indicates Drugs and Media are combinatorial +clear DM; +excLnNum=2; +drgCnt=0; +medCnt=0; +if ispc + while (~isequal(txt{excLnNum,2},'###')) + drgCnt=drgCnt+1; + DM.drug(drgCnt) = {raw(excLnNum,2)}; + DM.conc(drgCnt) = {raw(excLnNum,3)}; + excLnNum=excLnNum+1; + end + + while (~isequal(txt{excLnNum,4},'###')) + medCnt=medCnt+1; + DM.media(medCnt) = {raw(excLnNum,4)}; + excLnNum=excLnNum+1; + end +else %else if not PC (Then linux or other) + excLnNum=1; + while (~isequal(DMcell{excLnNum+1,2},'###')) + drgCnt=drgCnt+1; + DM.drug(drgCnt) = {DMtbl(excLnNum,2)}; + DM.conc(drgCnt) = {DMtbl(excLnNum,3)}; + excLnNum=excLnNum+1; + end + + while (~isequal(DMcel{excLnNum+1,4},'###')) + medCnt=medCnt+1; + DM.media(medCnt) = {DMtbl(excLnNum,4)}; + excLnNum=excLnNum+1; + end +end + +msg=strcat('NumberOfDrugs = ',num2str(drgCnt), ' NumberOfMedias = ',num2str(medCnt) ) + +end + +save (fullfile(MPdir,'MPDMmat'), 'fields','MP','DM','Linked'); + +cd(w) +msgbox(['Drug-Media-MasterPlate Annotation File Generation Complete']) \ No newline at end of file diff --git a/workflow/templates/easy/DgenResults.m b/workflow/templates/easy/DgenResults.m new file mode 100755 index 00000000..9481be1f --- /dev/null +++ b/workflow/templates/easy/DgenResults.m @@ -0,0 +1,727 @@ +%% CALLED BY EASYconsole.m %% +%to 'Fotos' from 'PTmaps' +%Console globals******* +global ExpPath +global resDir +global ExpOutmat +global ImParMat + +w=pwd; +ln=1; + + +%******Version compatability fixes*******assoc'd with diff and solve v11a='7.12.0.635 (R2011a)'; +if verLessThan('matlab','8.3') %original work 23_0227 updated 23_0525 (8.4 changed to 8.3) + fd4=diff(sym('K / (1 + exp(-r* (t - l )))'),4); + sols= solve(fd4); +else %accomodate new matlab changes after 2014a fix 2nd update 23_0227 + syms t K r l; + fd4= diff(K / (1 + exp(-r* (t - l ))),t,4); + sols= solve(fd4); + tmpswap= sols(1); + sols(1)= sols(3); + sols(3)= tmpswap; +end + +%MPnum=ImParMat(1); +%opt = questdlg('Print Results Only (RES), DB Only (DB), or Both','Results Printout Options','Res','DB','Both','Both'); +opt='Res'; + +if ~exist('ImParMat','var') || isempty(ImParMat) + load ImParameters +end +destPerMP=ImParMat(2); +load (ExpOutmat); +load (fullfile(ExpPath,'MasterPlateFiles','MPDMmat.mat')) +numOfDrgs= length(DM.drug); +numOfMeds= length(DM.media); +destPerMP=numOfDrgs; + +%Determine the longest set of intensity(tPts) for the Experiment Data Set +maxNumIntens=0; +for n=1:size(scan,2) + for m=1:size(scan(n).plate,2) + maxNumIntens= max(maxNumIntens,size(scan(n).plate(m).intens,2)); + end +end +if length(ExpPath)== max(strfind(ExpPath,'\')) + localExpPath=ExpPath(1:end-1); +else + localExpPath=ExpPath; +end +if ispc %Linux accommodation +expNm=localExpPath(max(strfind(localExpPath,'\'))+1:end); +else + if isunix,expNm=localExpPath(max(strfind(localExpPath,'/'))+1:end);end + if ismac,expNm=localExpPath(max(strfind(localExpPath,'/'))+1:end);end +end +if ispc +drivePos=min(strfind(localExpPath,'\')); +else + drivePos= min(strfind(localExpPath,'/')); +end +drive= localExpPath(1:(drivePos-1)); + +DBupload=fullfile(drive,'EZdbFiles','DBupLOADfiles'); +%DBupload=['G:\EZdbFiles\DBupLOADfiles'] +%resultsFilename= strcat(resDir,'\PrintResults\!!Results_Output.txt') %'\PrintResults\zDevelCFitOutIndxSelectedRawDatTimesDrugMedia.txt' +%DBfilename= strcat(resDir,'\PrintResults\!!Dbase_Output.txt') % Print !Dbase file in PrintResults dir +%resultsFilename= strcat(resDir,'\PrintResults\!!Results_',expNm,'.txt') +%DBfilename= strcat(resDir,'\PrintResults\!!Dbase_',expNm,'.txt') + +%Added to allow backward compatability************************************ +%Test for CFoutStd as indication of 2018version with r_refined fit code;If +%earlier version with only a standard composite fite, Print results to +%!!ResultsStd_...txt only +%*********************** +try + scan(1).plate(1).CFoutStd(1,1); + resultsFilename= fullfile(resDir,'PrintResults', strcat('!!ResultsELr_',expNm,'.txt')); + DBfilename= fullfile(resDir,'PrintResults', strcat('!!DbaseELr_',expNm,'.txt')); +catch + resultsFilename= fullfile(resDir,'PrintResults', strcat('!!ResultsStd_',expNm,'.txt')); + DBfilename= fullfile(resDir,'PrintResults', strcat('!!DbaseStd_',expNm,'.txt')); +end +%************************************************************************** + +if isequal(opt,'Res')||isequal(opt,'Both'),fid = fopen(resultsFilename,'w');end +if isequal(opt,'DB')||isequal(opt,'Both'),fid2 = fopen(DBfilename,'w');end %121012 Combo + +if isequal(opt,'Res')||isequal(opt,'Both') %print Results + fprintf(fid,'%d\t',ln); %Results header + fprintf(fid,'%s\t\n',localExpPath); + ln=ln+1; + fprintf(fid,'%d\t',ln); +end + +mpCnt=0; +totPlCnt=0; +drgCnt=0; +medCnt=0; +%load(fullfile(resDir,'PTmats','Nbdg')) %Convolute scan array data into plates +try +load(fullfile(resDir,'Fotos','Nbdg')) %Convolute scan array data into plates +catch +load(fullfile(resDir,'PTmats','Nbdg')) %Convolute scan array data into plates +end +for s=1:size(scan,2) + %Convolute scan array data into plates DconB for DBcombo + clear Diag + try + Diag(:,:,:,1)= sbdg{s}(1:1:24,16:-1:1,:); + catch + sbdg{s}; + end + %***************************************************** + for p=1:size((scan(s).plate),2) + totPlCnt=totPlCnt+1; + + if destPerMP>1 &&rem(totPlCnt,destPerMP)==1, mpCnt=mpCnt+1; end + if destPerMP==1,mpCnt=mpCnt+1; end + pertCnt= rem(totPlCnt,destPerMP); + if pertCnt==0, pertCnt= destPerMP;end + pert= strcat('Perturb_',num2str(pertCnt)); + s + p + clear outCmat + outCmat=scan(s).plate(p).CFout; + %Print Time Point HEADER for each plate for newly added intensity data + if isequal(opt,'Res')||isequal(opt,'Both') + fprintf(fid, '\n'); + ln=ln+1; + fprintf(fid,'%d\t',ln); + fprintf(fid,'Scan\tPlate\tRow\tCol\t'); + try + asd=cell2mat(scan(s).plate(1).CFparameters(1)); + aucEndPt= strcat('AUC',num2str(asd(9))); + catch + asd=cell2mat(scan(s).plate(1).CFparameters{1,1}(1,384)); + aucEndPt= strcat('AUC',num2str(asd(9))); + end + + fprintf(fid, 'Num.\tDiagnostics\tDrug\tConc\tMedia\tModifier1\tConc1\tModifier2\tConc2\tORF\tGene'); + fprintf(fid, '\t %s',aucEndPt); + fprintf(fid, '\triseTm\tK\tr\tl\tR-squared\tK-lower\tK-upper\tr-lower\tr-upper\tl-lower\tl-upper\tArea\tLastInten\tSplineMaxRateTm\tLastFitTm\t1stFitTm\tMedianBase\tFitBase\tMinTm\tThreshTm'); + + if size(outCmat,2)==27 + fprintf(fid, '\ttc11Cut\ttc12Cut\ttc21Cut\ttc22Cut'); %'\tEarly1\tEarly2\tLate1\tLate2'); 17_0629 MinBaseIntens update for MedianBase label + end + fprintf(fid, '\tTotFitPts\tPostThreshFitPts\t1stBackgrd\tLstBackgrd\t1stMeanTotBackgrd\tLstMeanTotBackgrd'); + + + end + clear outTseries + outTseries=[]; + outTseries=scan(s).plate(p).tSeries; + TseriesSize= size(outTseries,1); + + %*************************************************************** + %clear outCmat + %outCmat=scan(s).plate(p).CFout; + clear outIntens + outIntens=[]; + RawIntens=[]; + RawIntens=scan(s).plate(p).intens; + RawIntensSize= size(RawIntens,2) + + clear Ag; %Ag is Growth Area + Ag=scan(s).plate(p).Ag; + AgSize= size(Ag); +%*************************************** +dataLength= min(TseriesSize,RawIntensSize); + %for j=1:length(outTseries)%(size(outTseries,2) + if isequal(opt,'Res')||isequal(opt,'Both') + for j=1:dataLength + fprintf(fid, '\t%.5f', outTseries(j)); + end + end +%numBlkColTm=(maxNumIntens - size(outTseries,1)); + numBlkCol=(maxNumIntens - dataLength); %size(outTseries,1)); + %for nn=1:numBlkColTm %extend to col beyond longest rawDataSet + if isequal(opt,'Res')||isequal(opt,'Both') + for nn=1:numBlkCol %extend to col beyond longest rawDataSet + fprintf(fid, '\t'); + end + + fprintf(fid,'\tOrfRep'); + fprintf(fid,'\tSpecifics'); + fprintf(fid,'\tStrainBkGrd'); + fprintf(fid, '\n'); + ln=ln+1; + fprintf(fid,'%d\t',ln); + end +%*****************************************Begin Data Section + + n=0; + for r=1:16 + + for c=1:24 + n=n+1; + clear selcode; + Kval=outCmat(n,3); + rSq=outCmat(n,6); + lval=outCmat(n,5); + if Kval>160, selcode='K_Hi'; else selcode=' ';end + if Kval<40, selcode=strcat(selcode,' K_Lo'); end %12_1030 + if rSq<.97 && rSq>0, selcode=strcat(selcode,' rSqLo'); end + if lval>(0.85*(max(outTseries))), selcode=strcat(selcode,' late');end + if isnan(outCmat(n,7))||isnan(outCmat(n,8))||isnan(outCmat(n,9))... + ||isnan(outCmat(n,10))||isnan(outCmat(n,11))... + ||isnan(outCmat(n,12)), selcode=strcat(selcode,' NaN'); + end + + %*************RiseTime Calculation***************************************** + K= (outCmat(n,3)); + R= (outCmat(n,4)); + L= (outCmat(n,5)); + + if R>0 && L>0 && K>0 + rr=R; ll=L; + tc1= eval(sols(2)); + tc2= eval(sols(3)); + LL= eval(sols(1)); + riseTm= LL-tc1; + else + riseTm=0; + end +%************************************************************************* + + if Ag(n)< .30*(scan(s).Awindow),selcode=strcat(selcode,' smArea'); end %12_1030 + if outCmat(n,3)==0,selcode=strcat('0 ',selcode); end + + orf=cell2mat(MP(mpCnt).orf{1}(n)); + gene=cell2mat(MP(mpCnt).genename{1}(n)); + orfRep=cell2mat(MP(mpCnt).orfRep{1}(n)); + specifics=cell2mat(MP(mpCnt).specifics{1}(n)); + strain=cell2mat(MP(mpCnt).strain{1}(n)); + drug= char(DM.drug{pertCnt}); + conc= char(DM.conc{pertCnt}); + media= char(DM.media{pertCnt}); + try + mod1= char(DM.mod1{pertCnt}); + conc1= char(DM.conc1{pertCnt}); + catch + mod1= ' '; + conc1= ' '; + end + try + mod2= char(DM.mod2{pertCnt}); + conc2= char(DM.conc2{pertCnt}); + catch + mod2= ' '; + conc2= ' '; + end + + if ~isempty(outCmat) + if isequal(opt,'Res')||isequal(opt,'Both') + fprintf(fid,'%d\t %d\t %d\t %d\t %d\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t',s,p,r,c,n,selcode,drug,conc,media,mod1,conc1,mod2,conc2,orf,gene); + fprintf(fid, '%.5f\t %.5f\t %.5f\t %.5f\t %.5f\t %.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f',... + outCmat(n,1),riseTm,outCmat(n,3),outCmat(n,4),... + outCmat(n,5),outCmat(n,6),outCmat(n,7),outCmat(n,8),... + outCmat(n,9),outCmat(n,10),outCmat(n,11),outCmat(n,12),... + outCmat(n,13),outCmat(n,14),outCmat(n,15),outCmat(n,16)); + + fprintf(fid, '\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f',... + outCmat(n,17),outCmat(n,18),outCmat(n,19),... + outCmat(n,20),outCmat(n,21)); + %failPt= [s n r c] + + %Added for data cut times used in 'r'optomized method 06/14/2018 + if size(outCmat,2)==27 + fprintf(fid, '\t%.5f\t%.5f\t%.5f\t%.5f',... + outCmat(n,24),outCmat(n,25),outCmat(n,26),outCmat(n,27)); + end + fprintf(fid, '\t%d\t%d\t%d\t%d\t%d\t%d',... + outCmat(n,22),outCmat(n,23),Diag(c,r,1,p),Diag(c,r,2,p),Diag(c,r,3,p),Diag(c,r,4,p)); %,Diag(r,c,3,p),Diag(r,c,4,p)); + end %Results print + %DBfile******************************************* + %riseTm calc***** + if isequal(opt,'DB')||isequal(opt,'Both') + dbRsq= 0;dbKup= 0; dbKlo= 0; dbrup= 0; dbrlo= 0; dbLlo= 0; dbLup= 0; + if isnumeric(outCmat(n,6)), dbRsq= outCmat(n,6);end + if isnumeric(outCmat(n,7)), dbKup= outCmat(n,7);end + if isnumeric(outCmat(n,8)), dbKlo= outCmat(n,8);end + if isnumeric(outCmat(n,9)), dbrup= outCmat(n,9);end + if isnumeric(outCmat(n,10)), dbrlo= outCmat(n,10);end + if isnumeric(outCmat(n,11)), dbLlo= outCmat(n,11);end + if isnumeric(outCmat(n,12)), dbLup= outCmat(n,12);end + end + + %**************** + if isequal(opt,'DB')||isequal(opt,'Both') + fprintf(fid2,'%s\t %d\t %d\t %d\t %d\t %d\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t',expNm,s,p,r,c,n,selcode,drug,conc,media,mod1,conc1,mod2,conc2,orf,gene); + fprintf(fid2, '%.5f\t %.5f\t %.5f\t %.5f\t %.5f\t %.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f',... + outCmat(n,1),riseTm,outCmat(n,3),outCmat(n,4),... + outCmat(n,5),dbRsq,dbKup,dbKlo,... + dbrup,dbrlo,dbLlo,dbLup); %\t%.5f\t%.5f\t%.5f\t%.5f + %... outCmat(n,13),outCmat(n,14),outCmat(n,15),outCmat(n,16)); + + end +%*****DB analysis data end************************************************ + %Add Intensities series to end of curve fit data + outIntens=[]; + outIntens=zeros(384,dataLength); + intensBlob=''; + tmBlob =''; + + for j=1:dataLength %size(RawIntens,2) %size(outTseries,1) + if Ag(n)==0,Ag(n)=scan(s).Awindow;end + outIntens(n,j)= RawIntens(n,j)/Ag(n); + if isequal(opt,'Res')||isequal(opt,'Both') + fprintf(fid, '\t%.5f', outIntens(n,j)); % Results print Intens + end + if isequal(opt,'DB')||isequal(opt,'Both') + if j1 &&rem(totPlCnt,destPerMP)==1, mpCnt=mpCnt+1; end + if destPerMP==1,mpCnt=mpCnt+1; end + pertCnt= rem(totPlCnt,destPerMP); + if pertCnt==0, pertCnt= destPerMP;end + pert= strcat('Perturb_',num2str(pertCnt)); + s + + %Print Time Point HEADER for each plate for newly added intensity data + if isequal(opt,'Res')||isequal(opt,'Both') + fprintf(fid, '\n'); + ln=ln+1; + fprintf(fid,'%d\t',ln); + fprintf(fid,'Scan\tPlate\tRow\tCol\t'); + + try + asd=cell2mat(scan(s).plate(1).CFparameters(1)); + aucEndPt= strcat('AUC',num2str(asd(9))); + catch + asd=cell2mat(scan(s).plate(1).CFparameters{1,1}(1,384)); + aucEndPt= strcat('AUC',num2str(asd(9))); + end + fprintf(fid, 'Num.\tDiagnostics\tDrug\tConc\tMedia\tModifier1\tConc1\tModifier2\tConc2\tORF\tGene') + fprintf(fid, '\t %s',aucEndPt) + fprintf(fid, '\triseTm\tK\tr\tl\tR-squared\tK-lower\tK-upper\tr-lower\tr-upper\tl-lower\tl-upper\tArea\tLastInten\tSplineMaxRateTm\tLastFitTm\t1stFitTm\tMedianBase\tFitBase\tMinTm\tThreshTm\tTotFitPts\tPostThreshFitPts\t1stBackgrd\tLstBackgrd\t1stMeanTotBackgrd\tLstMeanTotBackgrd'); %17_0629 MinBaseIntens update for MedianBase label + end + clear outTseries + outTseries=[]; + outTseries=scan(s).plate(p).tSeries; + TseriesSize= size(outTseries,1); + + %*************************************************************** + clear outCmat + outCmat=scan(s).plate(p).CFoutStd; + clear outIntens + outIntens=[]; + RawIntens=[]; + RawIntens=scan(s).plate(p).intens; + RawIntensSize= size(RawIntens,2) + clear Ag; %Ag is Growth Area + Ag=scan(s).plate(p).Ag; + AgSize= size(Ag); +%*************************************** +dataLength= min(TseriesSize,RawIntensSize); + if isequal(opt,'Res')||isequal(opt,'Both') + for j=1:dataLength + fprintf(fid, '\t%.5f', outTseries(j)); + end + end + + numBlkCol=(maxNumIntens - dataLength); %size(outTseries,1)); + + if isequal(opt,'Res')||isequal(opt,'Both') + for nn=1:numBlkCol %extend to col beyond longest rawDataSet + fprintf(fid, '\t'); + end + + fprintf(fid,'\tOrfRep'); + fprintf(fid,'\tSpecifics'); + fprintf(fid,'\tStrainBkGrd'); + fprintf(fid, '\n'); + ln=ln+1; + fprintf(fid,'%d\t',ln); + end +%*****************************************Begin Data Section + + n=0; + for r=1:16 + + for c=1:24 + n=n+1; + clear selcode; + Kval=outCmat(n,3); + rSq=outCmat(n,6); + lval=outCmat(n,5); + if Kval>160, selcode='K_Hi'; else selcode=' ';end + if Kval<40, selcode=strcat(selcode,' K_Lo'); end %12_1030 + if rSq<.97 && rSq>0, selcode=strcat(selcode,' rSqLo'); end + if lval>(0.85*(max(outTseries))), selcode=strcat(selcode,' late');end + if isnan(outCmat(n,7))||isnan(outCmat(n,8))||isnan(outCmat(n,9))... + ||isnan(outCmat(n,10))||isnan(outCmat(n,11))... + ||isnan(outCmat(n,12)), selcode=strcat(selcode,' NaN'); + end + + %*************RiseTime Calculation***************************************** + K= (outCmat(n,3)); + R= (outCmat(n,4)); + L= (outCmat(n,5)); + if R>0 && L>0 && K>0 + rr=R; ll=L; + tc1= eval(sols(2)); + tc2= eval(sols(3)); + LL= eval(sols(1)); + riseTm= LL-tc1; + else + riseTm=0; + end + +%************************************************************************* + + if Ag(n)< .30*(scan(s).Awindow),selcode=strcat(selcode,' smArea'); end %12_1030 + if outCmat(n,3)==0,selcode=strcat('0 ',selcode); end + + orf=cell2mat(MP(mpCnt).orf{1}(n)); + gene=cell2mat(MP(mpCnt).genename{1}(n)); + orfRep=cell2mat(MP(mpCnt).orfRep{1}(n)); + specifics=cell2mat(MP(mpCnt).specifics{1}(n)); + strain=cell2mat(MP(mpCnt).strain{1}(n)); + drug= char(DM.drug{pertCnt}); + conc= char(DM.conc{pertCnt}); + media= char(DM.media{pertCnt}); + try + mod1= char(DM.mod1{pertCnt}); + conc1= char(DM.conc1{pertCnt}); + catch + mod1= ' '; + conc1= ' '; + end + try + mod2= char(DM.mod2{pertCnt}); + conc2= char(DM.conc2{pertCnt}); + catch + mod2= ' '; + conc2= ' '; + end + + if ~isempty(outCmat) + if isequal(opt,'Res')||isequal(opt,'Both') + fprintf(fid,'%d\t %d\t %d\t %d\t %d\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t',s,p,r,c,n,selcode,drug,conc,media,mod1,conc1,mod2,conc2,orf,gene); + fprintf(fid, '%.5f\t %.5f\t %.5f\t %.5f\t %.5f\t %.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f',... + outCmat(n,1),riseTm,outCmat(n,3),outCmat(n,4),... + outCmat(n,5),outCmat(n,6),outCmat(n,7),outCmat(n,8),... + outCmat(n,9),outCmat(n,10),outCmat(n,11),outCmat(n,12),... + outCmat(n,13),outCmat(n,14),outCmat(n,15),outCmat(n,16)); + + fprintf(fid, '\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%d\t%d',... + outCmat(n,17),outCmat(n,18),outCmat(n,19),... + outCmat(n,20),outCmat(n,21),outCmat(n,22),outCmat(n,23)); + + fprintf(fid, '\t%d\t%d\t%d\t%d',... + Diag(c,r,1,p),Diag(c,r,2,p),Diag(c,r,3,p),Diag(c,r,4,p)); %,Diag(r,c,3,p),Diag(r,c,4,p)); + end %Results print + %DBfile******************************************* + %riseTm calc***** + if isequal(opt,'DB')||isequal(opt,'Both') + dbRsq= 0;dbKup= 0; dbKlo= 0; dbrup= 0; dbrlo= 0; dbLlo= 0; dbLup= 0; + if isnumeric(outCmat(n,6)), dbRsq= outCmat(n,6);end + if isnumeric(outCmat(n,7)), dbKup= outCmat(n,7);end + if isnumeric(outCmat(n,8)), dbKlo= outCmat(n,8);end + if isnumeric(outCmat(n,9)), dbrup= outCmat(n,9);end + if isnumeric(outCmat(n,10)), dbrlo= outCmat(n,10);end + if isnumeric(outCmat(n,11)), dbLlo= outCmat(n,11);end + if isnumeric(outCmat(n,12)), dbLup= outCmat(n,12);end + end + + %**************** + if isequal(opt,'DB')||isequal(opt,'Both') + fprintf(fid2,'%s\t %d\t %d\t %d\t %d\t %d\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t',expNm,s,p,r,c,n,selcode,drug,conc,media,mod1,conc1,mod2,conc2,orf,gene); + fprintf(fid2, '%.5f\t %.5f\t %.5f\t %.5f\t %.5f\t %.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f',... + outCmat(n,1),riseTm,outCmat(n,3),outCmat(n,4),... + outCmat(n,5),dbRsq,dbKup,dbKlo,... + dbrup,dbrlo,dbLlo,dbLup); %\t%.5f\t%.5f\t%.5f\t%.5f + %... outCmat(n,13),outCmat(n,14),outCmat(n,15),outCmat(n,16)); + + end +%*****DB analysis data end************************************************ + %Add Intensities series to end of curve fit data + outIntens=[]; + outIntens=zeros(384,dataLength); + intensBlob=''; + tmBlob =''; + + for j=1:dataLength %size(RawIntens,2) %size(outTseries,1) + if Ag(n)==0,Ag(n)=scan(s).Awindow;end + outIntens(n,j)= RawIntens(n,j)/Ag(n); + if isequal(opt,'Res')||isequal(opt,'Both') + fprintf(fid, '\t%.5f', outIntens(n,j)); % Results print Intens + end + if isequal(opt,'DB')||isequal(opt,'Both') + if j + if exist(fullfile(resDir,'PTmats','ImParameters.mat')) + load(fullfile(resDir,'PTmats','ImParameters.mat')); + else + curDir=pwd; + returnHome + load ImParameters.mat + cd(curDir) + end + end + + matDir= fullfile(openExppath,'\'); + mkdir(fullfile(matDir,'BkUp')); + + % create supporting files + if ~exist(fullfile(resDir,'PrintResults')) + mkdir(fullfile(resDir,'PrintResults')); + end + if ~exist(fullfile(resDir,'figs')) + mkdir(fullfile(resDir,'figs')); + end + if ~exist(fullfile(resDir,'CFfigs')) + mkdir(fullfile(resDir,'CFfigs')); + end + if ~exist(fullfile(resDir,'PTmats')) + mkdir(fullfile(resDir,'PTmats')); + end + if ~exist(fullfile(resDir,'Fotos')) + mkdir(fullfile(resDir,'Fotos')); + end + + catch + returnHome + end + clear scan + + if exist('resDir','var')&&~isempty(resDir) + fhconsole= gcf; + set(fhconsole,'Name',strcat('EASYconsole- ',char(resDir))) + else + set(fhconsole,'Name','EASYconsole -Exp. Analysis NOT selected.') + end + +%% CALLBACKS %% +% callback for the 'Run' in the dropdown menu +function run_Callback(~, ~, ~) + returnHome + +function runPlateMapPintool_Callback(~, ~, ~) + returnHome + try + NImapPT + catch ME + returnHome + EASYconsole + end + +function NImCFcombo_Callback(~, ~, ~) + returnHome + try + par4Gbl_Main8c + catch + returnHome + EASYconsole + end + +function runPlateImAnal_Callback(~, ~, ~) + returnHome + try + NImStartupOnly + catch ME + returnHome + EASYconsole + end + +function PlateCFit_Callback(~, ~, ~) + returnHome + try + NCstart + catch ME + returnHome + end + +function GenPrintouts_Callback(~, ~, ~) + +function uploadExcelMP2DB_Callback(~, ~, ~) + +function runDMPexcel_Callback(~, ~, ~) + try + DMPexcel2mat + catch ME + returnHome + EASYconsole + end + +function runResults_DBcombo_Callback(~, ~, ~) + try + DgenResults %similar but semicolons removed to restore so cmdLine display info. + %DgenResults %par4global -convert 1x1cell of 384cells to be like previous 1x384 cells CFparameter + catch ME + disp('Error in DgenResults') + returnHome + EASYconsole + end + +function Tools_Callback(~, ~, ~) + +function runOverlayPlots_Callback(~, ~, ~) + returnHome + try + DoverlayPlots2 + returnHome + EASYconsole + catch ME + returnHome + EASYconsole + end + +function runFotoStrip_Callback(~, ~, ~) + returnHome + try + F_NImStartup_CentCir + returnHome + EASYconsole + catch ME + returnHome + EASYconsole + end + +function runDisplayFig_Callback(~, ~, ~) + returnHome + try + UfigDisplay + catch ME + returnHome + EASYconsole + end + +function runViewParameters_Callback(~, ~, ~) + returnHome + try + % This was blank in original + catch ME + returnHome + EASYconsole + end + +function QkviewN_Callback(~, ~, ~) + global ExpPath + returnHome + try + try + cd fullfile(ExpPath) + catch + if ispc cd c:\ + else + cd (fullfile('~')); + end + end + QkviewImages + catch ME + returnHome + EASYconsole + end + +function CFdisplay_Callback(~, ~, ~) + returnHome + try + NCsingleDisplay + returnHome + EASYconsole + catch ME + returnHome + EASYconsole + end diff --git a/workflow/templates/easy/NCdisplayGui.m b/workflow/templates/easy/NCdisplayGui.m new file mode 100755 index 00000000..dc74bafa --- /dev/null +++ b/workflow/templates/easy/NCdisplayGui.m @@ -0,0 +1,497 @@ +%% CALLED WHEN ACCESSING 'CurveFit Display' %% +function [scLst, row, col] = NCdisplayGui(eDir) %(ExpPath) +%global ExpPath +%global expDir + + xPos=0.05; + btnWid=0.10; + btnHt=0.05; + spacing=0.02;% Spacing between the button and the next command's label + +%hSpot= figure; +%figure +%==================================== + % The ADD Groups button + btnNumber=1; + yPos=0.85-(btnNumber-1)*(btnHt+spacing); + btnPos=[xPos yPos-spacing btnWid btnHt]; +% +row=1; +hedit = uicontrol(... + 'Style', 'edit',... + 'String',row,... + 'Units','normalized',... + 'Position', btnPos,... % [.002 .70 .08 .10],... + 'callback',{@editRowNum}); % 'Position', [5 100 60 20]) + + function editRowNum(source,~) +user_entry = str2double(get(source,'string')); + if (isnan(user_entry)||(user_entry<0)||(user_entry>17) ) + errordlg('Enter a Row between 1 and 16','Bad Input','modal') + return + end + row= user_entry; + + end +%-------------------222222----------- +% +btnNumber=2; + yPos=0.85-(btnNumber-1)*(btnHt+spacing); + btnPos=[xPos yPos-spacing btnWid btnHt]; + + +col= 1; +hedit = uicontrol(... + 'Style', 'edit',... + 'String',col,... + 'Units','normalized',... + 'Position', btnPos,... % [.002 .70 .08 .10],... + 'callback',{@entryColNum}); % 'Position', [5 100 60 20]) + + function entryColNum(source,~) +user_entry = str2double(get(source,'string')); + if (isnan(user_entry)||(user_entry<0)||(user_entry>25)) + errordlg('Enter a Column between 1 and 24','Bad Input','modal') + return + end + col= user_entry; + + end +%} + +%************************************************************* + +%Ncode 12_0120 for reading in numeric folder names +nlist=dir(fullfile(eDir,'*')); %(ExpPath,'*')); +nnn=0; +for n=1:size(nlist,1) + if (~isempty(str2num(nlist(n).name))) + nnn=nnn+1; + PnumLst(nnn)= (str2num(nlist(n).name)); + slst(nnn,1)={(nlist(n).name)}; + end +end + +hListbox = uicontrol(... + 'Style', 'listbox',... + 'String',sort(slst),... + 'value',[],... + 'max',1000,... + 'min',1,... + 'Units','normalized',... + 'Position', [.40 .40 .10 .60],... + 'callback',{@load_listbox}); %'uiresume(gcbf)'); 'Position', [5 100 60 20]) + + function load_listbox(source,~) + userIndx = (get(source,'value')); + userStr = (get(source,'string')); + %scLstIndx= str2num(char(strrep(userStr(userIndx), 'Scan', ''))) + + user_entry=userStr(userIndx); + scLst= user_entry; + end + +%***************************************************************** +%-------------------10 10 10 10----------- + +btnNumber=10; + yPos=0.85-(btnNumber-1)*(btnHt+spacing); + btnPos=[xPos yPos-spacing btnWid btnHt]; + +hedit8 = uicontrol(... + 'Style', 'pushbutton',... + 'String',{'Continue'},... + 'Units','normalized',... + 'Position', btnPos,... + 'callback','uiresume(gcbf)'); + + +%------------------------------ + +%*************************************************************** +% { +%********LABELS****************************** + % The Labels + xLPos=0.175; + yPos=0; + btnWid=0.20; + + + % Removed Not needed for Ncode ImRobot + lblNumber=1; + yPos=0.85-(lblNumber-1)*(btnHt+spacing); + btnPos=[xLPos yPos-spacing btnWid btnHt]; + +htxt = uicontrol(... + 'Style', 'text',... + 'String','Row',... + 'Units','normalized',... + 'Position', btnPos); +%-------------------222222----------- + lblNumber=2; + yPos=0.85-(lblNumber-1)*(btnHt+spacing); + btnPos=[xLPos yPos-spacing btnWid btnHt]; + +htxt = uicontrol(... + 'Style', 'text',... + 'String','Column',... + 'Units','normalized',... + 'Position', btnPos); + % Not needed for Ncode ImRobot + + uiwait(gcf); +end %function end $$$$$ +%} + %{ +%-------------------333333----------- + + lblNumber=3; + yPos=0.85-(lblNumber-1)*(btnHt+spacing); + btnPos=[xLPos yPos-spacing btnWid btnHt]; + +htxt = uicontrol(... + 'Style', 'text',... + 'String','BG Threshold (%above) Detection',... + 'Units','normalized',... + 'Position', btnPos); + +%-------------------4----------- + + lblNumber=4; + yPos=0.85-(lblNumber-1)*(btnHt+spacing); + btnPos=[xLPos yPos-spacing btnWid btnHt]; + +htxt = uicontrol(... + 'Style', 'text',... + 'String','SpotDetThres(1-60%)',... + 'Units','normalized',... + 'Position', btnPos); + +%-------------------55555----------- + lblNumber=5; + yPos=0.85-(lblNumber-1)*(btnHt+spacing); + btnPos=[xLPos yPos-spacing btnWid btnHt]; + +htxt = uicontrol(... + 'Style', 'text',... + 'String','Radius',... %'String','Width',... + 'Units','normalized',... + 'Position', btnPos); +%-------------------66666----------- +lblNumber=6; + yPos=0.85-(lblNumber-1)*(btnHt+spacing); + btnPos=[xLPos yPos-spacing btnWid btnHt]; + +htxt = uicontrol(... + 'Style', 'text',... + 'String','Dither',... + 'Units','normalized',... + 'Position', btnPos); +%-------------------77777----------- + +lblNumber=7; + yPos=0.85-(lblNumber-1)*(btnHt+spacing); + btnPos=[xLPos yPos-spacing btnWid btnHt]; + +htxt = uicontrol(... + 'Style', 'text',... + 'String','SearchRange',... + 'Units','normalized',... + 'Position', btnPos); + +%-------------------88888----------- +%{ +lblNumber=8; + yPos=0.85-(lblNumber-1)*(btnHt+spacing); + btnPos=[xLPos yPos-spacing btnWid btnHt]; + +htxt = uicontrol(... + 'Style', 'text',... + 'String','Blank2',... + 'Units','normalized',... + 'Position', btnPos); +%---------------------------------------- + %} +%} + + + +%{ + + +%-------------------66666----------- +btnNumber=6; + yPos=0.85-(btnNumber-1)*(btnHt+spacing); + btnPos=[xPos yPos-spacing btnWid btnHt]; +srcExtendFactor=ImParMat(7); +hedit = uicontrol(... + 'Style', 'edit',... + 'String',srcLoIntensThres,... + 'Units','normalized',... + 'Position', btnPos,... + 'callback',{@entryExtendFactor}); + + function entryExtendFactor(source,eventdata) +user_entry = str2double(get(source,'string')); + if (isnan(user_entry)||(user_entry<1.8)||(user_entry>4.0)) + errordlg('You must enter a numeric value between 1.8 and 2.1','Bad Input','modal') + return + end + ExtendFactor=user_entry + ImParMat(7)= ExtendFactor + ExtendFactor + end +%} + + +%{ +%-------------------333333----------- + + btnNumber=3; + yPos=0.85-(btnNumber-1)*(btnHt+spacing); + btnPos=[xPos yPos-spacing btnWid btnHt]; +srcBGthreshold=ImParMat(3); +hedit = uicontrol(... + 'Style', 'edit',... + 'String',srcBGthreshold,... + 'Units','normalized',... + 'Position', btnPos,... % [.002 .70 .08 .10],... + 'callback',{@entryBGthreshold}); % 'Position', [5 100 60 20]) + + function entryBGthreshold(source,eventdata) +user_entry = str2double(get(source,'string')); + if (isnan(user_entry)||(user_entry<1)||(user_entry>100)) + errordlg('Enter a numeric value between 1 and 100 percent to produce a Background Threshold value as a percent above the time series average background for each spot.','Bad Input','modal') + return + end + BGthresInput=user_entry + ImParMat(3)= BGthresInput +BGthresInput + end + +%-------------------444444----------- + +btnNumber=4; %Enter spot detection threshold (lock-in Image frame) + yPos=0.85-(btnNumber-1)*(btnHt+spacing); + btnPos=[xPos yPos-spacing btnWid btnHt]; +srcSpotThres=ImParMat(4); +hedit = uicontrol(... + 'Style', 'edit',... + 'String',srcSpotThres,... + 'Units','normalized',... + 'Position', btnPos,... + 'callback',{@entrySpotThres}); + + function entrySpotThres(source,eventdata) +user_entry = str2double(get(source,'string')); + if (isnan(user_entry)||(user_entry<1)||(user_entry>60)) + errordlg('You must enter a numeric value between 1 and 60','Bad Input','modal') + return + end + spotThres=user_entry + ImParMat(4)= spotThres + spotThres + end +% +%---------555555 Radius Entry After Sept.2014---------------------------------** + +btnNumber=5; + yPos=0.85-(btnNumber-1)*(btnHt+spacing); + btnPos=[xPos yPos-spacing btnWid btnHt]; +srcRadius=ImParMat(10); +hedit = uicontrol(... + 'Style', 'edit',... + 'String',srcRadius,... + 'Units','normalized',... + 'Position', btnPos,... % [.002 .70 .08 .10],... + 'callback',{@entryRadius}); % 'Position', [5 100 60 20]) + + function entryRadius(source,eventdata) +user_entry = str2double(get(source,'string')); + if (isnan(user_entry)||(user_entry<12)||(user_entry>17)) + errordlg('You must enter a numeric value between 12 and 17','Bad Input','modal') + return + end + Radius=user_entry + ImParMat(10)= Radius + Radius + end + +%---------555555 Width Entry prior the Sept.2014---------------------------------** +%{ +btnNumber=5; + yPos=0.85-(btnNumber-1)*(btnHt+spacing); + btnPos=[xPos yPos-spacing btnWid btnHt]; +srcWidth=ImParMat(5); +hedit = uicontrol(... + 'Style', 'edit',... + 'String',srcWidth,... + 'Units','normalized',... + 'Position', btnPos,... % [.002 .70 .08 .10],... + 'callback',{@entryWidth}); % 'Position', [5 100 60 20]) + + function entryWidth(source,eventdata) +user_entry = str2double(get(source,'string')); + if (isnan(user_entry)||(user_entry<5)||(user_entry>41)) + errordlg('You must enter a numeric value between 5 and 40','Bad Input','modal') + return + end + Width=user_entry + ImParMat(5)= Width + Width + end +%} +%-------------------66666 Dither unnecessary after Sept.2014----------- + + btnNumber=6; + yPos=0.85-(btnNumber-1)*(btnHt+spacing); + btnPos=[xPos yPos-spacing btnWid btnHt]; +srcDither= ImParMat(6); +hedit = uicontrol(... + 'Style', 'edit',... + 'String',srcDither,... + 'Units','normalized',... + 'Position', btnPos,... + 'callback',{@entryDither}); + + function entryDither(source,eventdata) +user_entry = str2double(get(source,'string')); + if (isnan(user_entry)||(user_entry<0)||(user_entry>5)) + errordlg('You must enter a numeric value between 1 and 4','Bad Input','modal') + return + end + Dither=user_entry + ImParMat(6)= Dither + Dither + end +%-------------------77777----------- Added July 7,2015 to allow Search Range constraint +btnNumber=7; + yPos=0.85-(btnNumber-1)*(btnHt+spacing); + btnPos=[xPos yPos-spacing btnWid btnHt]; + try + CSrchRange=ImParMat(12); + CSrchRng=ImParMat(12) + catch %Legacy default value was 18 before being made a user input variable (ImParMat(12)). A preferable value now might be 12 or 14. + CSrchRange=18; + ImParMat(12)=18 + CSrchRng=ImParMat(12) + end + + %{ + if size(scLst)>1 + CSrchRange=ImParMat(12); + else + try + CSrchRange=CSearchRange(str2double(scLst)) + catch + CSrchRange=ImParMat(12); + end + end + %} +hSearchRange = uicontrol(... + 'Style', 'edit',... + 'String',CSrchRange,... + 'Units','normalized',... + 'Position', btnPos,... + 'callback',{@CsearchRange}); + + function CsearchRange(source,eventdata) + user_entry = str2double(get(source,'string')); + if (isnan(user_entry)||(user_entry<1)||(user_entry>50)) %originally 18; 19_0729 increased + errordlg('You must enter a numeric value between 1 and 18 12->18 recommended. (ImParMat(12)))','Bad Input','modal') + return + end + CSrchRng=user_entry + + end + +%-------------------77777----------- +%{ +btnNumber=7; + yPos=0.85-(btnNumber-1)*(btnHt+spacing); + btnPos=[xPos yPos-spacing btnWid btnHt]; +srcExtend=ImParMat(7); + +hedit = uicontrol(... + 'Style', 'edit',... + 'String',srcExtend,... + 'Units','normalized',... + 'Position', btnPos,... + 'callback',{@entryExtend}); + + function entryExtend(source,eventdata) +user_entry = str2double(get(source,'string')); + if (isnan(user_entry)||(user_entry<-0.10)||(user_entry>0.4)) + errordlg('You must enter a numeric value between 0 and 0.4. 0.10 recommended','Bad Input','modal') + return + end + extend=user_entry + ImParMat(7)= extend + extend + end +%-------------------888888----------- + +btnNumber=8; + yPos=0.85-(btnNumber-1)*(btnHt+spacing); + btnPos=[xPos yPos-spacing btnWid btnHt]; +%ImParMat(8)=1; + srcpointExtend=ImParMat(8); +hedit = uicontrol(... + 'Style', 'edit',... + 'String',srcpointExtend,... + 'Units','normalized',... + 'Position', btnPos,... + 'callback',{@entrypointExtend}); + + function entrypointExtend(source,eventdata) +user_entry = str2double(get(source,'string')); +user_entry= floor(user_entry); + if (isnan(user_entry)||(user_entry<-3)||(user_entry>5)) + errordlg('You must enter an integer value between 0 and 5. 1 recommended','Bad Input','modal') + return + end + pointExtend=user_entry + ImParMat(8)= pointExtend + pointExtend + end +%} +%-------------------999999----------- + +btnNumber=9; + yPos=0.85-(btnNumber-1)*(btnHt+spacing); + btnPos=[xPos yPos-spacing btnWid btnHt]; + +hedit = uicontrol(... + 'Style', 'popupmenu',... + 'String',{'GrowthArea','FixedArea'},... + 'Units','normalized',... + 'Position', btnPos,... + 'callback',{@grwArea}); + + function grwArea(source,eventdata) + str = get(source, 'String'); + val = get(source,'Value'); + % Set current data to the selected data set. + switch str{val}; + case 'GrowthArea' ;% User selects Peaks. + SWgrowthArea = 1 + case 'FixedArea' % User selects Membrane. + SWgrowthArea = 0 + end + end +%} +%global SWsingleSc +%global SWgrowthArea +%global selScan +%global scan +%global scLst +%global ImParMat +%global CSearchRange +%global CSrchRng +%global defImParMat +%global fhImRun +%global fhconsole +%global resDir +%global ExpOutmat +%global numRows; +%global numCols; diff --git a/workflow/templates/easy/NCfitImCFparforFailGbl2.m b/workflow/templates/easy/NCfitImCFparforFailGbl2.m new file mode 100755 index 00000000..d5aaa3d5 --- /dev/null +++ b/workflow/templates/easy/NCfitImCFparforFailGbl2.m @@ -0,0 +1,275 @@ +%% CALLED BY par4GblFnc8c.m %% +function [par4scanselIntensStd,par4scanselTimesStd,par4scanTimesELr,par4scanIntensELr,par4scanCFparameters,par4scanCFdate,outC,outCstd]= ... + NCfitImCFparforFailGbl2(parMat,times, values, timeOffsets, fileSuffix, AUCfinalTime, ~, spotAreas, outputDirectory, ~,~, sols, ~) %,scan) + + +%Preallocation for parfor loop********************************** +st(1,1:size(times,2))= 1111; +resMat(1,1:27)= 0; +resMatStd= resMat; + outC= zeros(384,27); + outCstd= zeros(384,27); +for m=1:384 + pa{m}= st; + + par4scanCFparameters{m}= parMat; + par4scanCFdate{m}= datestr((now),31); +end + par4scanselTimesStd= pa; + par4scanselIntensStd= pa; + par4scanTimesELr= pa; + par4scanIntensELr= pa; + par4resMat= zeros(384,27); + par4resMatStd= zeros(384,27); +%**************************************************************** + +%*****************Begin the Spot(cultures) loop**************************** +for ii= 1:384 %startSpot:numCultures +%parfor ii= 1:384 %startSpot:numCultures + ii; %%%db print out the culture number + + timepts= []; + currValues=[]; + currSpotAreas=[]; + currSpotArea=[]; + + dataMatrix=[]; + selTimesStd=[]; %191024 parfor + selIntensStd=[]; %191024 parfor + FiltTimesELr=[]; %191024 parfor + NormIntensELr=[]; %191024 parfor + + % add offset...1 offset PER PLATE + timepts = times + timeOffsets; %(floor((ii-1)/arrayFormat) + 1); + currValues = values(ii,:); %change values(spotNum,:); + % get spot areas for this culture + currSpotArea = spotAreas(:,ii); + + % just use the area at the last time point + %currSpotArea = currSpotAreas(1); + + %-------------------------------------------------------------- + + %Preallocate to accomodate parfor loop ************************** + resMatStd= zeros(1,27); + resMat= zeros(1,27); + currNormIntens = currValues/currSpotArea; + tmpx=find(currNormIntens>5); %15jh % 2.3); + validSpot=true; + if(isempty(tmpx) || (length(tmpx)<3)) + validSpot= false; + normIntens= currNormIntens; + filterTimes= timepts; %filterTimes; %currTimes; + + selTimesStd= timepts; + selIntensStd= currNormIntens; + FiltTimesELr= timepts; + NormIntensELr= currNormIntens; + + else + +%******************************************************* +%NCfilImCF +%NCfilImCF.m called from NCfitImCF.m line 119******************* +%Preallocate incase something bails in NCscurImCFparfor +resMatStd= zeros(1,27); +resMat= zeros(1,27); + + hold off %09_0704 +dataMatrix=[]; + K=0;r=0;l=0;Klow=0;Kup=0;rlow=0;rup=0;llow=0;lup=0;AUC=0;MSR=0; rsquare=0; + bl=0; + Tpt1=0; numFitTpts =0;thresGT2=0;minTime=0;fitbl=0; %diagnostic outputs only + +timepts= timepts; %timepts= currTimes; parfor +normIntens= currNormIntens; + + dataMatrix= []; %parfor move clear from NCfitImCF...m + loIntensThres=parMat(4); + stdLoIntLim=parMat(5); + +%******Begin basic filter********** +%[loInten Thres,stdbased Trim before Scurve start, and dropout detection] + if(max(normIntens) > 2.29) + threshold = loIntensThres; %1.9; %Increase this value to reduce low data point (flag=2) + else + threshold = 0; + end + + dropThreshold= -0.0001* max(normIntens); +%Initialize dataMatrix********************************** + dataMatrix(1,:)=timepts; + dataMatrix(2,:)=normIntens; + dataMatrix(3,:)= ones; + dataMatrix(4,:)=normIntens; + +%Determine a mean Intensity Index point and assoc'd TimePt +a= min(normIntens(normIntens>=0)); %(find(normIntens>=0))); +b= max(normIntens(normIntens>=0)); %(find(normIntens>=0))); +c= 0.5*(b-a); +d= b-c; + +meanIntIndPt= find(normIntens>d,1); +meanInt= normIntens(meanIntIndPt); +meanTime= times(meanIntIndPt); +%********************************************************** + +%************************************************************************** +% NCLoIntstdTrim +%NCLoSstdTrim.m called by NCfilImCF and NCfil.m +flg1=0; +loScurvLim=stdLoIntLim; % +loStimeN=1; +last2n=1; + +stdDev= []; +nrmIntens0= normIntens; +for n=1:meanIntIndPt + if nrmIntens0(n)<=0 + nrmIntens0(n)=0; + end + + if(nrmIntens0(n)0,dataMatrix(3,1:(n-2))= 2; %add to lowIntensity cull flags, the pre S cull data + else dataMatrix(3,1:n)= 2;end + dataMatrix(3,1:(n-2))= 2; + last2n=n; + end + if n<(length(nrmIntens0)-3) + x=nrmIntens0(n:(n+3)); + stdDev(n)=std(x); + if (stdDev(n)6 ,flg1=1;end %detect S curve rise->stop stdLoInt detection + + end + +end %end for + +if (loStimeN-2)>0,dataMatrix(3,1:(loStimeN-2))= 2; %add to lowIntensity cull flags, the pre S cull data +else dataMatrix(3,1:(loStimeN-2))= 2;end + +qcutoff=2; +qind=find(normIntens>2); %,:,'first'); +if ~isempty(qind(3)), qcutoff=qind(3);end +[minInt,I]= min(normIntens(2:qcutoff)); +bl=minInt; +minTime=dataMatrix(1,I); %diagnostic output only + +if (length(qind)>5)&&I>1,dataMatrix(3,1:(I-1))=5; end + +tGT2Flg=0; +for n=1:length(normIntens) + dataMatrix(4,n)=normIntens(n)-bl; + if n>I && dataMatrix(4,n)>=2 && tGT2Flg==0 + thresGT2=n;thresGT2Tm=dataMatrix(1,n);tGT2Flg=1; + end %diagnostic output only +end + +resMat(18)= bl; +resMatStd(18)= bl; +resMatStd(20)= minTime; +resMat(20)= minTime; + +%&&&&&&&&&&&&&&&&&&&&&&& &&&&&&&&&&&&&&&&&&&&&&&& +%********************************************************* + %DropOut cull section(single drop points)******************* + DropOutStartPt=length(normIntens); + for n=1:length(normIntens) + if(n>1) + if(((normIntens(n)- normIntens(n-1))< dropThreshold))&& ... + (n > max(meanIntIndPt,thresGT2) ) + dataMatrix(3,n)= 6; + + end + end + end + +%??????should/could this be removed as recreated in%NCscurImCF_3parfor.m???????????????????? +%Post Stdtest cull for low intensities inclusion of additional low value points +%selTimes= [--,--] %don't know size before as it is a filtered output +tmpIndx= 0; + for n=1:length(normIntens) + if (dataMatrix(3,n)==1) + tmpIndx= tmpIndx+1; + selTimes(tmpIndx)= dataMatrix(1,n); %selTimes(nn)= dataMatrix(1,n); + selIntens(tmpIndx)= dataMatrix(4,n); %selIntens(nn)= dataMatrix(4,n); + + end + end + + +selTimes=selTimes'; +selIntens=selIntens'; +%??????????????????????????????????????????????????? +%**********End Basic filter****************** + + +%****************************************************** + + filtNormIntens=normIntens; + + dataMatrix0= dataMatrix; + filterTimes= timepts; %parfor + +%_____________________________________________________________ +%******************************************* +%NCscurImCF +%NCscurImCF_1 +%NCscurImCF_2 +%NCscurImCF_3 +%NCscurImCF_3parfor + %NCscurImCF_3parfor(dataMatrix, AUCfinalTime) + +%dataMatrix %debugging parfor gbl ok 85.7145; 126.4579,6, 124.5264 37tPt +%adsfj %debugging parfor gbl +[resMatStd, resMat, selTimesStd, selIntensStd, FiltTimesELr, NormIntensELr] =... + NCscurImCF_3parfor(dataMatrix0, AUCfinalTime, currSpotArea, sols, bl, minTime); + + end % end JUMP OVER associated with if(find(intensities>1000)<5) + %resMatStd + %asdf + %To accommodate parfor can't use global 'scan' variable 191002================ + % { + par4scanselTimesStd{ii}= selTimesStd; %timepts'; %filterTimes'; + par4scanselIntensStd{ii}= selIntensStd; %normIntens'; + par4scanTimesELr{ii}= FiltTimesELr; %19_1021 preserve for CurveDisplay and EZview + par4scanIntensELr{ii}= NormIntensELr; %19_1021 preserve for CurveDisplay and EZview + %} +%$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ +%************* striped down OutCell code was put back into NCfitImCF + + + outC(ii,:)= resMat; %{ii, par4resMat}; + outCstd(ii,:)= resMatStd; %{ii, par4resMatStd}; + +end %Multispot parfor ii loop end PARFOR LOOP END############################################################################## +%############################################################################################################################### +%############################################################################################################################### + + + +%outC + +%*********************************19_1001*********************************** +%To accomodate parfor copy par4scan thru global p4 functions inside of +%parfor loop --then outside to par4Gbl_Main8b.m + +%************************************************************************** +fileExt = '.txt'; +filePrefix = 'FitResultsComplete_'; +fileNamePlate = [filePrefix fileSuffix fileExt]; +fileName= fullfile(outputDirectory, fileNamePlate); %[outputDirectory fileNamePlate]; +fid = fopen(fileName,'w'); +fprintf(fid, 'Fit Results Complete\n'); +%fprintf(fid, 'Num.\tAUC\tMSR\tK\tr\tl\tR-squared\tK-lower\tK-upper\tr-lower\tr-upper\tl-upper\tl-lower\tArea\tLastInten\tSpineMaxRateTimePt\tLastFitTimePt\n'); +fclose(fid); + + +end %function end + +%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& +%########################################################################## + diff --git a/workflow/templates/easy/NCscurImCF_3parfor.m b/workflow/templates/easy/NCscurImCF_3parfor.m new file mode 100755 index 00000000..ce9b241d --- /dev/null +++ b/workflow/templates/easy/NCscurImCF_3parfor.m @@ -0,0 +1,437 @@ +%% CALLED BY NCfitImCFparforFailGbl2.m %% +function [resMatStd, resMat, selTimesStd, selIntensStd, FiltTimesELr, NormIntensELr] =... + NCscurImCF_3parfor(dataMatrix, AUCfinalTime, currSpotArea, sols, bl, minTime) +%Major revision for Early-Late data cuts to improve accuracof 'r'. Removed legacy iterative method. +%Significant Modification for Parfor +%*************************************************************** + +%########################################################################## + %******************************************* New Stage 1*************************************************************** +%Preallocate +resMatStd= zeros(1,27); +resMat= zeros(1,27); +%Set internal variables sent to matlab fit function +me=200; +meL=750; +mi=25; %50 +miL=250; +%*********************************** + rmsStg1=0; + rmsStg1I(1)= 0; + slps=1; + + filterTimes=[]; + normIntens=[]; + nn=1; + numFitTpts=0; +%Build filterTimes and normIntens from spot dataMatrix selection codes produced in filter section +for n=1:size(dataMatrix,2) + if (((dataMatrix(3,n)==1))||(dataMatrix(3,n)==3)||(dataMatrix(3,n)==2)... + ||(dataMatrix(3,n)==0)) + filterTimes(nn)= dataMatrix(1,n); + normIntens(nn)=dataMatrix(4,n); + nn=nn+1; + end +end + +%------------------------------------------------------------------ + +%++++++++++++++++++++++++++++++++++++++++ + +filterTimes=filterTimes'; +selTimesStd=filterTimes; +normIntens=normIntens'; +selIntensStd=normIntens; + +%normIntens %debugging parfor gbl 200330 good values +%afgj + + lastTptUsed= 1; + lastIntensUsed= 1; + thresGT2TmStd=0; + try + + lastTptUsed=max(filterTimes); + lastIntensUsed=normIntens(length(normIntens)); + lastIntensUsedStd= lastIntensUsed; + + lastTptUsedStd= lastTptUsed; + Tpt1Std= filterTimes(1); + numFitTptsStd=nnz((normIntens(:)>=0)==1); + thresGT2 = find(((normIntens(:)>2)==1), 1); + if isempty(thresGT2) + thresGT2TmStd=0; + else + thresGT2TmStd = filterTimes(thresGT2); + end + numTptsGT2Std = nnz((normIntens(:)>=2)==1); %nnz(filterTimes(find(filterTimes>=thresGT2Tm))); + K_Guess = max(normIntens); + numTimePts = length(filterTimes); + opts = fitoptions('Method','Nonlinear','Robust','On',... + 'DiffMinChange',1.0E-11,'DiffMaxChange',0.001,... + 'MaxFunEvals',me, 'MaxIter', mi, 'TolFun', 1.0E-12, 'TolX', 1.0E-10, 'Lower', [K_Guess*0.5,0,0], 'StartPoint', [K_Guess,filterTimes(floor(numTimePts/2)),0.30], 'Upper', [K_Guess*2.0,max(filterTimes),1.0],'Display','off'); + ftype = fittype('K / (1 + exp(-r* (t - l )))','independent','t','dependent',['K','r','l'],'options',opts); + + % carry out the curve fitting process + [fitObject, errObj] = fit(filterTimes,normIntens,ftype); + coeffsArray = coeffvalues(fitObject); + rmsStg1=errObj.rsquare; + rmsStg1I(slps)= errObj.rsquare; + sDat(slps,1)=errObj.rsquare; + K = coeffsArray(1); sDat(slps,2)= coeffsArray(1); % Carrying Capacity + l = coeffsArray(2); sDat(slps,3)= coeffsArray(2); %lag time + r = coeffsArray(3); sDat(slps,4)= coeffsArray(3); % rateS + + % integrate (from first to last time point) + numVals = size(filterTimes); + numVals = numVals(1); + t_begin = 0; + t_end = AUCfinalTime; + AUC = (K/r*log(1+exp(-r*(t_end-l)))-K/r*log(exp(-r*(t_end-l)))) - (K/r*log(1+exp(-r*(t_begin-l)))-K/r*log(exp(-r*(t_begin-l)))); + MSR = r; + + rsquare = errObj.rsquare; + confObj = confint(fitObject,0.9); % get the 90% confidence + + NANcond=0; stdNANcond=0; %stdNANcond added to relay not to attempt ELr as there is no curve to find critical point + confObj_filtered=confObj; + Klow = confObj(1,1); sDat(slps,5)=confObj(1,1); + Kup = confObj(2,1); sDat(slps,6)=confObj(2,1); + llow = confObj(1,2); sDat(slps,7)=confObj(1,2); + lup = confObj(2,2); sDat(slps,8)=confObj(2,2); + rlow = confObj(1,3); sDat(slps,9)=confObj(1,3); + rup = confObj(2,3); sDat(slps,10)=confObj(2,3); + if(isnan(Klow)||isnan(Kup)||isnan(llow)||isnan(lup)||isnan(rlow)||isnan(rup)) + NANcond=1; stdNANcond=1; %stdNANcond added to relay not to attempt ELr as there is no curve to find critical point + end + + %rup %debugging parfor gbl 200330 + %Klow + %asdfjj114 +% { +catch %ME + % if no data is given, return zeros + AUC = 0;MSR = 0;K = 0;r = 0;l = 0;rsquare = 0;Klow = 0;Kup = 0; + rlow = 0;rup = 0;lup = 0;llow = 0; + NANcond=1; stdNANcond=1; %stdNANcond added to relay not to attempt ELr as there is no curve to find critical point +end %end Try +%} + if (exist('K','var')&& exist('r','var') && exist('l','var')) + t = (0:1:200); + Growth = K ./ (1 + exp(-r.* (t - l ))); + fitblStd= min(Growth); %jh diag + end + + cutTm(1:4)= 1000; %-1 means cuts not used or NA +%{ +l %debugging parfor +K +r +Klow +k +%} +%***Preserve for ResultsStd+++++++++++++++++++++++++++++ + +resMatStd(1)= AUC; +resMatStd(2)= MSR; +resMatStd(3)= K; +resMatStd(4)= r; +resMatStd(5)= l; +resMatStd(6)= rsquare; +resMatStd(7)= Klow; +resMatStd(8)= Kup; +resMatStd(9)= rlow; +resMatStd(10)= rup; +resMatStd(11)= llow; +resMatStd(12)= lup; +resMatStd(13)= currSpotArea; +resMatStd(14)= lastIntensUsedStd; %filtNormIntens(length(filtNormIntens)); + +%spline fit unneccessary and removed;therefore No maxRateTime assoc'd w/spline fit +%try +%resMatStd(15)= maxRateTime; +%catch + maxRateTime= 0; %[]; %Std shows []; ELr shows 0; %parfor + resMatStd(15)= 0; %maxRateTimestdMeth; +%end + + +resMatStd(16)= lastTptUsedStd; %filterTimes(length(filterTimes)); + if isempty(Tpt1Std) + Tpt1Std= 777; %0.000002; % + end +resMatStd(17)= Tpt1Std; +resMatStd(18)= bl; %perform in the filter section of NCfitImCFparfor +resMatStd(19)= fitblStd; %Taken from NCfil... and not affected by NCscur...changes +resMatStd(20)= minTime; %Not affected by changes made in NCscur...for refined 'r' +resMatStd(21)= thresGT2TmStd; +resMatStd(22)= numFitTptsStd; +resMatStd(23)= numTptsGT2Std; + +resMatStd(24)= 999; %The Standard method has no cuts .:.no cutTm +resMatStd(25)= 999; +resMatStd(26)= 999; +resMatStd(27)= 999; + + +%if SwitchEvsEL==3 %Remove 'SwitchEvsEL==...' temporary SWITCH when Hartman decides what he wants +%********************************************************************************* +%ELr New Experimental data through L+deltaS Logistic fit for 'Improved r' Fitting +%********************************************************************************* +FiltTimesELr= []; %{ii}= filterTimes; +NormIntensELr= []; %{ii}= normIntens; +normIntens= selIntensStd; +filterTimes= selTimesStd; +stdIntens= selIntensStd; +tmpIntens= selIntensStd; +stdTimes= selTimesStd; + +if stdNANcond==0 +%Determine critical points and offsets for selecting Core Data based on +%Standard curve fit run. Put diff into NImStartupImCF02.m calling source +%to reduce repeated execution since it doesn't change. + %fd4=diff(sym('K / (1 + exp(-r* (t - l )))'),4); + %sols=solve(fd4); +tc1= eval(sols(2)); +tc2= eval(sols(3)); +LL= l; %eval(sols(1)); %exactly the same as 'l' from std. fit method-Save time +rsTmStd= LL-tc1; %%riseTime (first critical point to L) +deltS= rsTmStd/2; + tc1Early= tc1-deltS; %AKA- tc1AdjTm %2*tc1 -LL + L_Late= LL+deltS; + + + tc1EdatPt= find(filterTimes>(tc1Early),1,'first'); + cutTm(1)= filterTimes(2); + cutDatNum(1)= 2; + cutTm(2)= tc1Early; + cutDatNum(2)= tc1EdatPt-1; + + L_LDatPt= find(filterTimes< L_Late,1,'last'); + tc2LdatPt= find(filterTimes< tc2+rsTmStd,1,'last'); + + cutTm(3)= L_Late; + cutDatNum(3)= L_LDatPt; + +%Select Core Data Set (Remove Early data before critical point) +ints=[]; +ints(1:L_LDatPt-tc1EdatPt+2)= (tmpIntens(L_LDatPt)); +ints(2:end)= tmpIntens(tc1EdatPt:L_LDatPt); +ints(1)= tmpIntens(1); +tms=[]; +tms(1:L_LDatPt-tc1EdatPt+2)= (stdTimes(L_LDatPt)); +tms(2:end)= stdTimes(tc1EdatPt:L_LDatPt); +tms(1)= stdTimes(1); +%----------------------------------------------- +%Include/Keep late data that define K ********* +if length(tmpIntens(tc2LdatPt:end))> 4 + KlastInts= tmpIntens(tc2LdatPt:end); + KlastTms= stdTimes(tc2LdatPt:end); + lengthKlast= length(tmpIntens(tc2LdatPt:end)); + ints(end:(end+ lengthKlast-1))= KlastInts; + tms(end:(end+ lengthKlast-1 ))= KlastTms; + + cutTm(4)= tc2+rsTmStd; + cutDatNum(4)= tc2LdatPt-1; + + +else + lengthKlast= length(tmpIntens(tc2LdatPt-1:end)); + if lengthKlast>1 + KlastInts= tmpIntens(end-(lengthKlast-1):end); + KlastTms= stdTimes(end-(lengthKlast-1):end); + ints(end:(end+ lengthKlast-1 ))= KlastInts; + tms(end:(end+ lengthKlast-1 ))= KlastTms; + end + +cutTm(4)= stdTimes(tc2LdatPt-1); +cutDatNum(4)= tc2LdatPt-2; %length(stdTimes(end-(lengthKlast-1):end)); + +end + + +%************************************************ +Ints=[]; +Tms=[]; +Ints= ints'; +Tms= tms'; + +try + filterTimes= Tms; filterTimes4= Tms; + normIntens= Ints; normIntens4=Ints; + +%---------------------------------------------------------------------------------------------- + + %classic symetric logistic curve fit setup restated as COMMENTS for reference convenience---------------------------- + %opts= fitoptions is the same as for Std and so is redundant + %opts = fitoptions('Method','Nonlinear','Robust','On',... + % 'DiffMinChange',1.0E-11,'DiffMaxChange',0.001,... + % 'MaxFunEvals',me, 'MaxIter', mi, 'TolFun', 1.0E-12, 'TolX', 1.0E-10, 'Lower', [K_Guess*0.5,0,0], 'StartPoint', [K_Guess,filterTimes(floor(numTimePts/2)),0.30], 'Upper', [K_Guess*2.0,max(filterTimes),1.0]); + + ftype = fittype('K / (1 + exp(-r* (t - l )))','independent','t','dependent',['K','l','r'],'options',opts); + fitObject=[]; errObj=[]; + % carry out the curve fitting process + [fitObject, errObj] = fit(Tms,Ints,ftype); + coeffsArray = coeffvalues(fitObject); + r3 = coeffsArray(3); %sDat(slps,4)= coeffsArray(3); % rateS + + if (exist('K','var')&& exist('r','var') && exist('l','var')) + t = (0:1:200); + GrowthELr = K ./ (1 + exp(-r.* (t - l ))); + fitblELr= min(GrowthELr); %jh diag + end + + catch ME + % if no data is given, return zeros + AUC = 0;MSR = 0;K = 0;r = 0;l = 0;rsquare = 0;Klow = 0;Kup = 0; + rlow = 0;rup = 0;lup = 0;llow = 0; %normIntens=[]; + +end %end Try + +end %if stdNANcond=0 + %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + %Update values if r is better(higher) with removal of early data +try + if r3>r && stdNANcond==0 + r = r3; sDat(slps,4)= sDat(slps,4); % rateS + K = coeffsArray(1); sDat(slps,2)= coeffsArray(1); % Carrying Capacity + l = coeffsArray(2); sDat(slps,3)= coeffsArray(2); %lag time + + coeffsArray = coeffvalues(fitObject); + rmsStg1=errObj.rsquare; + rmsStg1I(slps)= errObj.rsquare; + sDat(slps,1)=errObj.rsquare; + + %jh diagnostics************* + numFitTpts=nnz((normIntens(:)>=0)==1); + thresGT2 = find(((normIntens(:)>2)==1), 1); + thresGT2Tm = filterTimes(thresGT2); + numTptsGT2 = nnz((normIntens(:)>=2)==1); + numTimePts = length(filterTimes); + + AUC = (K/r*log(1+exp(-r*(t_end-l)))-K/r*log(exp(-r*(t_end-l)))) - (K/r*log(1+exp(-r*(t_begin-l)))-K/r*log(exp(-r*(t_begin-l)))); + MSR = r3; + %*************************** + + rsquare = errObj.rsquare; + confObj = confint(fitObject,0.9); % get the 90% confidence + NANcond=0; + + confObj_filtered=confObj; + Klow = confObj(1,1); sDat(slps,5)=confObj(1,1); + Kup = confObj(2,1); sDat(slps,6)=confObj(2,1); + llow = confObj(1,2); sDat(slps,7)=confObj(1,2); + lup = confObj(2,2); sDat(slps,8)=confObj(2,2); + rlow = confObj(1,3); sDat(slps,9)=confObj(1,3); + rup = confObj(2,3); sDat(slps,10)=confObj(2,3); + if(isnan(Klow)||isnan(Kup)||isnan(llow)||isnan(lup)||isnan(rlow)||isnan(rup)) + NANcond=1; + end + + filterTimes= Tms; + normIntens= Ints; + resMat(17)= .00002; + resMat(18)= bl; + resMat(19)= fitblELr; + resMat(20)= minTime; + else % r is better than r3 so use the Std data in the ELr result sheet + filterTimes=selTimesStd; + normIntens=selIntensStd; + + lastTptUsed= lastTptUsedStd; %Reinstall Std values for jh diags + + Tpt1=filterTimes(1); + try + if isempty(Tpt1) + Tpt1= 0.00002; %777; + end + catch + Tpt1= 0.00002; %777; + end + +resMat(17)= Tpt1; + numFitTpts= numFitTptsStd; + numTptsGT2= numTptsGT2Std; + thresGT2Tm = thresGT2TmStd; + + cutTm(1:4)= 1000; %-1 means cuts not used or NA + + resMat(18)= bl; %only applicable to Std curve Fit; ELr superceeds and makes meaningless + resMat(19)= fitblStd; %only applicable to Std curve Fit; ELr superceeds and makes meaningless + resMat(20)= minTime; %only applicable to Std curve Fit; ELr superceeds and makes meaningless + end %if r3>r1 + + %rup + %asdf352 +% { +catch ME + % if no data is given, return zeros + AUC = 0;MSR = 0;K = 0;r = 0;l = 0;rsquare = 0;Klow = 0;Kup = 0; + rlow = 0;rup = 0;lup = 0;llow = 0; %normIntens=[]; + +end %end Try +%} +resMat(1)= AUC; +resMat(2)= MSR; +resMat(3)= K; +resMat(4)= r; +resMat(5)= l; +resMat(6)= rsquare; +resMat(7)= Klow; +resMat(8)= Kup; +resMat(9)= rlow; +resMat(10)= rup; +resMat(11)= llow; +resMat(12)= lup; +resMat(13)= currSpotArea; +resMat(14)= lastIntensUsed; %filtNormIntens(length(filtNormIntens)); + +%spline fit unneccessary and removed therfor no max spline rate time->set 0 + maxRateTime=0; %ELr will show 0; Std will show [] + resMat(15)= maxRateTime; + +resMat(16)= lastTptUsed; %filterTimes(length(filterTimes)); + + +try %if Std fit used no cuts .:. no cutTm +resMat(24)= cutTm(1); +resMat(25)= cutTm(2); +resMat(26)= cutTm(3); +resMat(27)= cutTm(4); +catch +resMat(24)= 999; %if Std fit used no cuts .:. no cutTm +resMat(25)= 999; +resMat(26)= 999; +resMat(27)= 999; +end + +FiltTimesELr= filterTimes; +NormIntensELr= normIntens; + +%********************************************************************************** +%########################################################################## + + lastTptUsed=max(filterTimes); + lastIntensUsed=normIntens(length(normIntens)); + + if (exist('K','var')&& exist('r','var') && exist('l','var')) + t = (0:1:200); + Growth = K ./ (1 + exp(-r.* (t - l ))); + fitbl= min(Growth); %jh diag + end + %} + try + if isempty(thresGT2Tm),thresGT2Tm=0;end %jh diag + catch + thresGT2Tm= 0; + numTptsGT2= 0; + end + +resMat(21)= thresGT2Tm; +resMat(22)= numFitTpts; +resMat(23)= numTptsGT2; + + +end %function end + + \ No newline at end of file diff --git a/workflow/templates/easy/NCsingleDisplay.m b/workflow/templates/easy/NCsingleDisplay.m new file mode 100755 index 00000000..6d90d290 --- /dev/null +++ b/workflow/templates/easy/NCsingleDisplay.m @@ -0,0 +1,110 @@ +%% CALLED WHEN ACCESSING 'CurveFit Display' %% +global expDir +global scLst +global ExpPath +global scan + +eDir= ExpPath; +hf= figure; +%**************Parameter Entry****************** +[scLst, row, col] = NCdisplayGui(eDir); %(ExpPath); %Ncode 122111replaced removed ,numOfPrtTimes) +%*********************************************** + +close(hf) + + +selSpot= (row-1)*24 + col; +for iPlate=1:length(scLst) + scanPltNum= str2double(scLst(iPlate)); + %scanPltNum(iPlate); +%----------------------------------------------------------------------- +K= scan(scanPltNum).plate(1).CFout((selSpot),3); +r= scan(scanPltNum).plate(1).CFout((selSpot),4); +l= scan(scanPltNum).plate(1).CFout((selSpot),5); + +suffix= strcat('Scan-Plate', scLst(iPlate)); %char(QspLst(n)); +%fileSpotSuffix = strcat('-Spot#',num2str(selSpot),'-Row=',selSpotRC(1),'-Col=',selSpotRC(2),'-FitData','-L=',num2str(l),'-r=',num2str(r),'-K=',num2str(K)); +fileSpotSuffix = strcat('-Spot#',num2str(selSpot),'-Row=',num2str(row),'-Col=',num2str(col),'-FitData','-L=',num2str(l),'-r=',num2str(r),'-K=',num2str(K)); +filenameNoExt= [suffix fileSpotSuffix]; + +timeArr= scan(scanPltNum).plate(1).tSeries; +rawIntens= scan(scanPltNum).plate(1).intens((selSpot),:)/scan(scanPltNum).plate(1).Ag((selSpot)); +try +filterTms= scan(scanPltNum).plate(1).filterTimes{(selSpot)}; +normInts= scan(scanPltNum).plate(1).normIntens{(selSpot)}; +catch + end +%******************************** + + if (exist('K','var')&& exist('r','var') && exist('l','var')) + t = (0:1:200); + Growth = K ./ (1 + exp(-r.* (t - l ))); + end +if length(scLst)>1 + figure +else + cla +end + +hold on +plot(timeArr,rawIntens,'g*'); +try +plot(filterTms,normInts,'o'); +catch +end +hold on; +title(strcat(filenameNoExt,'')); +xlabel('Hours'); +ylabel('Intensities Normalized by Area'); +grid on; + +if (exist('K','var')&& exist('r','var') && exist('l','var')) +plot(t, Growth,'b-'); +%*********Plot L on curvefit figure*********** + grL=Growth(round(l)); %growthCurve timePT for l in hours + plot(l,0:grL,'-b') %to display position of l + plot(l,grL,'^b') %to display l position on curve as an arrowhead + +%***Plot Arbitrary User Entry AUC "finalTimePt" + % plot(finalTimePt,0,'+m') + % plot(0:finalTimePt,bl,'-m') + +end +end +%end %end for SWsingle display/plot figure + +%############################################################ + + + +%Spot entry form------------------------------------------------------ +%{ +prompt = {'Enter spot to analyse:'}; +dlg_title = 'Input spot to curve fit'; +num_lines = 1; +def = {'1'}; +selSpot = inputdlg(prompt,dlg_title,num_lines,def); +K= scan(scanPltNum).plate(1).CFout(str2double(selSpot),3); +r= scan(scanPltNum).plate(1).CFout(str2double(selSpot),4); +l= scan(scanPltNum).plate(1).CFout(str2double(selSpot),5); + +suffix= strcat('Scan-Plate', scLst); %char(QspLst(n)); +fileSpotSuffix = strcat('-Spot#',selSpot,'-FitData','-L=',num2str(l),'-r=',num2str(r),'-K=',num2str(K)); +filenameNoExt= [suffix fileSpotSuffix]; + +timeArr= scan(scanPltNum).plate(1).tSeries; +rawIntens= scan(scanPltNum).plate(1).intens(str2double(selSpot),:)/scan(scanPltNum).plate(1).Ag(str2double(selSpot)); +filterTms= scan(scanPltNum).plate(1).filterTimes{str2double(selSpot)}; +normInts= scan(scanPltNum).plate(1).normIntens{str2double(selSpot)}; + +%} +%----------------------------------------------------------------------- +%{ +prompt = {'Enter Spot row:','Enter Spot column:'}; +dlg_title = 'Input spot to curve fit'; +num_lines = 2; +def = {'1','1'}; +selSpotRC = inputdlg(prompt,dlg_title,num_lines,def); +row= str2double(selSpotRC(1)); col=str2double(selSpotRC(2)); +%} +%row= cell2mat(row); diff --git a/workflow/templates/easy/NIcircle.m b/workflow/templates/easy/NIcircle.m new file mode 100755 index 00000000..4443c459 --- /dev/null +++ b/workflow/templates/easy/NIcircle.m @@ -0,0 +1,62 @@ +%% CALLED BY par4Gbl_Main8c.m %% +%{ +%Imaging ToolBox method +r=14; +A=zeros(70,70); %(fIntsc(refPtR:(refPtRExt),refPtC:(refPtCExt))) +m= {40,40}; +A(m{:})=1; +B= imdilate(A,strel('disk',r,0) ); +imshow(B) + +area=pi*r^2 + +clear all +%} + +%without Image Proc. Toolbox +%r=14; +%A=zeros(70,70); +%A=zeros(r,r); +%P=[40,40]; +%center= [refPtR+ round(.5*width), refPtC+ round(.5*width)]; +%A=zeros(70,70); +%--------------------------------------------------------------------- +%radius=14; +diaExt=2*(radius+1); +circBoxA=zeros(diaExt,diaExt); +center= [radius+2, radius+2]; +[m, n ]=size(circBoxA); +X = bsxfun(@plus,(1:m)', zeros(1,n)); +Y = bsxfun(@plus,(1:n), zeros(m,1)); +cirMask= sqrt(sum(bsxfun(@minus,cat(3,X,Y),reshape(center,1,1,[])) .^2,3))<=radius; +area=pi*radius^2; + +cirPixA= nnz(cirMask); + +optCirMask= double(cirMask); +optCirMask(optCirMask==0)= 0.8; + +%+++++++++++Foto Circle Fram(e)+++++++++++++++++++++ +expansion=2; +radExpan= radius+expansion; +FdiaExt=2*(radExpan); +circBoxA=zeros(FdiaExt,FdiaExt); +center= [radExpan+1, radExpan+1]; +[m, n ]=size(circBoxA); +X = bsxfun(@plus,(1:m)', zeros(1,n)); +Y = bsxfun(@plus,(1:n), zeros(m,1)); +FcirMask= sqrt(sum(bsxfun(@minus,cat(3,X,Y),reshape(center,1,1,[])) .^2,3))<=radExpan; + +%FcirPixA= nnz(cirMask); +FoptCirMask= double(FcirMask); +FoptCirMask(FoptCirMask==1)= 2; +%FoptCirMask(FoptCirMask==0)= 1; +%********Combine Masks to create circular boundry************ +padOptCirMask= padarray(optCirMask,[expansion-1 expansion-1],0.8); +FoptCirMask= FoptCirMask .* padOptCirMask; +FoptCirMask(FoptCirMask==1.6)= 0.8; +FoptCirMask(FoptCirMask==0)= 1; +FoptCirMask(FoptCirMask==2)= 1; + +%--------------------------------------------------- +%imagesc(cirMask) \ No newline at end of file diff --git a/workflow/templates/easy/NImParamRadiusGui.m b/workflow/templates/easy/NImParamRadiusGui.m new file mode 100755 index 00000000..aebac776 --- /dev/null +++ b/workflow/templates/easy/NImParamRadiusGui.m @@ -0,0 +1,281 @@ +%% CALLED BY par4Gbl_Main8c.m %% +function NImParamRadiusGui(expDir) %, numOfPrtTimes) 122111replace removed +global SWsingleSc +global SWgrowthArea +%global selScan +global scan +global scLst +global ImParMat +%global CSearchRange +global CSrchRng +global defImParMat +global fhImRun +global fhconsole +global resDir +global ExpOutmat +global numRows; +global numCols; +global scanSize +global scanMax +defImParMat=[1, 1, 15, 34, 24, 1,0,0,1,14,1,18]; %Ncode ImRobot adaptation +if ImParMat(3)==0 || ImParMat(4)==0 ||ImParMat(5)==0 || ImParMat(10)==0 ||ImParMat(11)==0 + ImParMat=defImParMat; +end +if size(ImParMat,2)<12 + ImParMat(12)=18; %Default before user input CsearchRange value + msg= 'Data made before SearchRange user entry added (ImParMat(12). 18 was the set value and the current default.)'; +end +%ImParMat=defImParMat; %Activate for INITIAL USE only +MPnum=1; +destPerMP=1; +selScan=1; +SWgrowthArea=1; + + if exist(fullfile(resDir,'PTmats','NImParameters.mat')) + load(fullfile(resDir,'PTmats','NImParameters')); + else + load NImParameters + end + ImParMat; +%end +%if ~exist('CSearchRange','var') || isempty(CSearchRange) + if ~isequal(exist(fullfile(resDir,'Fotos','CSearchRange.mat')),0) + load(fullfile(resDir,'Fotos','CSearchRange')) + CSearchRange; + end + + +% yInitPos=0.30; + xPos=0.05; + btnWid=0.10; + btnHt=0.05; + spacing=0.02;% Spacing between the button and the next command's label + + +%==================================== + % The ADD Groups button + btnNumber=1; + yPos=0.85-(btnNumber-1)*(btnHt+spacing); + btnPos=[xPos yPos-spacing btnWid btnHt]; +fhImParm=gcf; +if exist('resDir','var')&& ~isempty(resDir) + set(fhImParm,'NumberTitle','off') + set(fhImParm,'Name',strcat('ImageAnalysis- ',char(resDir))) +else + set(fhImParm,'NumberTitle','off') + set(fhImParm,'Name','EASYconsole -Exp. Analysis NOT selected.') +end + + + +btnNumber=5; + yPos=0.85-(btnNumber-1)*(btnHt+spacing); + btnPos=[xPos yPos-spacing btnWid btnHt]; +srcRadius=ImParMat(10); +hedit = uicontrol(... + 'Style', 'edit',... + 'String',srcRadius,... + 'Units','normalized',... + 'Position', btnPos,... % [.002 .70 .08 .10],... + 'callback',{@entryRadius}); % 'Position', [5 100 60 20]) + + function entryRadius(source,~) +user_entry = str2double(get(source,'string')); + if (isnan(user_entry)||(user_entry<12)||(user_entry>17)) + errordlg('You must enter a numeric value between 12 and 17','Bad Input','modal') + return + end + Radius=user_entry; + ImParMat(10)= Radius; + Radius; + end + + + btnNumber=6; + yPos=0.85-(btnNumber-1)*(btnHt+spacing); + btnPos=[xPos yPos-spacing btnWid btnHt]; +srcDither= ImParMat(6); +hedit = uicontrol(... + 'Style', 'edit',... + 'String',srcDither,... + 'Units','normalized',... + 'Position', btnPos,... + 'callback',{@entryDither}); + + function entryDither(source,~) +user_entry = str2double(get(source,'string')); + if (isnan(user_entry)||(user_entry<0)||(user_entry>5)) + errordlg('You must enter a numeric value between 1 and 4','Bad Input','modal') + return + end + Dither=user_entry; + ImParMat(6)= Dither; + Dither; + end +%-------------------77777----------- Added July 7,2015 to allow Search Range constraint +btnNumber=7; + yPos=0.85-(btnNumber-1)*(btnHt+spacing); + btnPos=[xPos yPos-spacing btnWid btnHt]; + try + CSrchRange=ImParMat(12); + CSrchRng=ImParMat(12); + catch %Legacy default value was 18 before being made a user input variable (ImParMat(12)). A preferable value now might be 12 or 14. + CSrchRange=18; + ImParMat(12)=18; + CSrchRng=ImParMat(12); + end + +hSearchRange = uicontrol(... + 'Style', 'edit',... + 'String',CSrchRange,... + 'Units','normalized',... + 'Position', btnPos,... + 'callback',{@CsearchRange}); + + function CsearchRange(source,~) + user_entry = str2double(get(source,'string')); + if (isnan(user_entry)||(user_entry<1)||(user_entry>50)) %originally 18; 19_0729 increased + errordlg('You must enter a numeric value between 1 and 18 12->18 recommended. (ImParMat(12)))','Bad Input','modal') + return + end + CSrchRng=user_entry; + + end + +%Ncode 12_0120 for reading in numeric folder names +nlist=dir(fullfile(expDir,'*')); +nnn=0; +for n=1:size(nlist,1) + if (~isempty(str2num(nlist(n).name))) + nnn=nnn+1; + PnumLst(nnn)= (str2num(nlist(n).name)); + sl(nnn,1)={(nlist(n).name)}; + end +end +scanSize= size(sl,1); +scanMax= max(str2double(sl)); + +hListbox = uicontrol(... + 'Style', 'listbox',... + 'String',sort(sl),... + 'value',[],... + 'max',1000,... + 'min',1,... + 'Units','normalized',... + 'Position', [.70 .40 .10 .60],... + 'callback',{@load_listbox}); %'uiresume(gcbf)'); 'Position', [5 100 60 20]) + + function load_listbox(source,~) + %global CSrchRng + %global CSrearchRange +userIndx = (get(source,'value')); +userStr = (get(source,'string')); +%scLstIndx= str2num(char(strrep(userStr(userIndx), 'Scan', ''))) + +user_entry=userStr(userIndx); + scLst=user_entry; + if size(scLst,1)>1 + %CSrchRng=num2str(ImParMat(12)) + set(hSearchRange,'string',num2str(ImParMat(12))) + else + try + CSrchRng=CSearchRange(str2double(scLst)); + set(hSearchRange,'string',CSearchRange(str2double(scLst))) + catch + % CSrchRng=num2str(ImParMat(12)) + % set(hSearchRange,'string',num2str(ImParMat(12))) + end + end + + %ImParMat(7)=MPnum + + end + scLst; + +%***************************************************************** +%-------------------10 10 10 10----------- + +btnNumber=10; + yPos=0.85-(btnNumber-1)*(btnHt+spacing); + btnPos=[xPos yPos-spacing btnWid btnHt]; + +hedit8 = uicontrol(... + 'Style', 'pushbutton',... + 'String',{'Continue'},... + 'Units','normalized',... + 'Position', btnPos,... + 'callback','uiresume(gcbf)'); + + +%------------------------------ + +%*************************************************************** +%********LABELS****************************** + % The Labels + xLPos=0.175; + yPos=0; + btnWid=0.20; + + +%-------------------55555----------- + lblNumber=5; + yPos=0.85-(lblNumber-1)*(btnHt+spacing); + btnPos=[xLPos yPos-spacing btnWid btnHt]; + +htxt = uicontrol(... + 'Style', 'text',... + 'String','Radius',... %'String','Width',... + 'Units','normalized',... + 'Position', btnPos); +%-------------------66666----------- +lblNumber=6; + yPos=0.85-(lblNumber-1)*(btnHt+spacing); + btnPos=[xLPos yPos-spacing btnWid btnHt]; + +htxt = uicontrol(... + 'Style', 'text',... + 'String','Dither',... + 'Units','normalized',... + 'Position', btnPos); +%-------------------77777----------- + +lblNumber=7; + yPos=0.85-(lblNumber-1)*(btnHt+spacing); + btnPos=[xLPos yPos-spacing btnWid btnHt]; + +htxt = uicontrol(... + 'Style', 'text',... + 'String','SearchRange',... + 'Units','normalized',... + 'Position', btnPos); + + +uiwait(gcf); +for i=1:length(scLst) + CSearchRange(str2double(scLst(i)))= CSrchRng; + ImParMat(12)= CSrchRng; +end + +ImParMat; +CSearchRange; +try + save('NImParameters','ImParMat') + % save('CSearchRange','CSearchRange') +catch +save(fullfile('\','NImParameters'),'ImParMat') + %save(fullfile('\','CSearchRange'),'CSearchRange') +end +save((fullfile(resDir,'PTmats','NImParameters')), 'ImParMat'); +save((fullfile(resDir,'Fotos','CSearchRange')),'CSearchRange'); + +close + +return + + +end + + + + + diff --git a/workflow/templates/easy/NIscanIntensBGpar4GblFnc.m b/workflow/templates/easy/NIscanIntensBGpar4GblFnc.m new file mode 100755 index 00000000..9c002c7c --- /dev/null +++ b/workflow/templates/easy/NIscanIntensBGpar4GblFnc.m @@ -0,0 +1,1381 @@ +%% CALLED BY par4GblFnc8c.m %% +function [Tmpsbdg2, scanIntens, F_spots, bmtp, optomizedPos, TmpexpScanIntens2, TmpFexpScanSpots2, TmpFexpScanBMtp2, TmpanlZoneRefs2,areaOfIntensAboveBG]= ... + NIscanIntensBGpar4GblFnc(Fflg, tifFileLst, ImParMat,PTmapPos,optCirMask,diaExt,doCircle,cirPixA,numRows,numCols,ImHeigth,ImWidth,cirMask, ... + tptLength,selScan,Empsc,~, ~, ~, ~,resDir, Tmpsbdg1) + +global ExpOutmat +global CSrchRng + + +%*************************************************************** + +searchIntens=[]; intensMax=[]; detMaxPos=[]; scIntens=[]; areaOfIntensAboveBG=[]; BkgrdMat=[]; + +if Fflg==1 + %*****PreAllocation ***** +tPtLength=length(tifFileLst); + +optomizedPos= cell(24,16,tPtLength); +F_spots=cell(24,16,tPtLength); %Foto spot coordinates Preallocation +bmtp=zeros(24,16,tPtLength); %Added/used for FOTOS Preallocation + +MPnum=ImParMat(1); +destPerMP=ImParMat(2); +spotThres=ImParMat(4)/100; +%width=ImParMat(5); +CSrchRng = ImParMat(12); +width=24; %absolute fixed after Sept 2014 +ditherF=ImParMat(6); +ditherF=1; %absolute fixed after Sept 2014 + + +%Get PTmap data (detPos) and create PTimage with cirMask 2019_0904 + PTmapOnesDbl= zeros(2075,1400); + for r=1:24 + for c=1:16 + [PTrefPt]=PTmapPos{r,c}; + PTrefPtR=PTrefPt(1); + PTrefPtC=PTrefPt(2); + %******* doCircle + PTmapOnesDbl(PTrefPtR:(PTrefPtR+(diaExt-1)),PTrefPtC:(PTrefPtC+(diaExt-1)))=optCirMask; + end + end + + %********************************************************************* + +end %if Fflg==1 +%******* Zeroth Initial Search Parameters***************************************** + rRangeUpper=-12; + rRangeLower=12; + cRangeLower=-12; + cRangeUpper=12; + ditherI= width/3; %8; + +%***************************************** +searchRange= floor(width*0.75); %generally -/+18 +%searchRange=CSrchRng +rCRangeUpper= -searchRange; % -12 % -8; %Ranges should be multiples of dither +rCRangeLower= searchRange; %12; %Ranges should be multiples of dither +cCRangeLower= -searchRange; %-12 %-8; %-12; %Ranges should be multiples of dither +cCRangeUpper= searchRange; %12; %8; %Ranges should be multiples of dither +ditherC= floor((2*searchRange)/9); %gnerally /6 0r /9 -> dither=6 or 4 pixels +if ditherC==0, ditherC=1;end + + +%********************************************************************* +% +rFRangeUpper= -floor(1/2*width); %(searchRange*2/3) %-12 % -8; %Ranges should be multiples of dither +rFRangeLower= floor(1/2*width); %(searchRange*2/3); %12; %Ranges should be multiples of dither +cFRangeLower= -floor(1/2*width); %(searchRange*2/3); %-12 %-8; %-12; %Ranges should be multiples of dither +cFRangeUpper= floor(1/2*width); %(searchRange*2/3) %12; %8; %Ranges should be multiples of dither +ditherF=1; + +% +%********************************************************* +%********************************************************* +%The purpose is to reduce incidence of 'chasing' scratches and abrasions to +%the edge of a frame and also to limit late timept finds of overgrown +%spot(that are yet small). The idea is to make the threshold high enough to +%prevent the tracking of scratches and small contaminants(overgrowth) while +%not limiting detection of small slow grow cultures. Recommend value +%slightly above the expected background level. +noSpotThres=60; %If PixelIntensity for area < area*this, use the Orig.search Refpt +if doCircle==1 +noSpot= cirPixA*noSpotThres; +else +noSpot= (width^2)*noSpotThres; %Area * (a selected background pixel level i.e., 70) +end +%*************************************************************** + +%********************************************************** +%Begin Time series loop + +%preallocation +stopSearch=zeros(numRows,numCols); %set all spots 'stopSearch' to false +intens=zeros(24,16); +intensPrev=zeros(24,16); +cent = cell(24,16); + +BGTav=zeros(numRows,numCols); +BGthres(1:numRows,1:numCols)= 95; %70;before incr.to95 to accommodate dark media %50forEpson;%Initialize for first tPt run *****************************************************<< +if length(tifFileLst) >2 +%**************************************************************************** + +%******Determine a good Plate Image from all the timepoints to use for Registration +%Plate Intensity curve over time points +%clear plateImage + plateImage= {}; +% { +%******Version compatability fixes*******assoc'd with diff and solve +%{ +v11a='7.12.0.635 (R2011a)'; +v14a='8.3.0.532 (R2014a)'; +%v14b='8.4.0.150421 (R2014b)'; +v18a='9.4.0.813654 (R2018a)'; +v18b='9.5.0.944444 (R2018b)'; +pool= 2; +if isequal(v11a,version)|| isequal(v14a,version) + if matlabpool('size')==0 + matlabpool(pool) + end + +else + if parpool('local')==0 + parpool('local', pool) + end +end +%} +%parfor iii=1:length(tifFileLst) +for iii=1:length(tifFileLst) + %tifFileLst + %asdf + tifFile= char(tifFileLst(iii)); + fullplt= imread(tifFile); %Ncode actually reading a .bmp file for Ncode + fullpl= 255 - fullplt(:,:,1); + %fullpl= (fullplt(1:ImHeigth,1:ImWidth)); + pltInts(iii)= sum(sum(fullpl(:,:)))/(ImHeigth*ImWidth); + plateImage{iii}= fullpl; %190919 to use instead of re-reading in NIscanIntensBG.m +end +guessOffsetAdj= 0.5; %Half(0.5) is the midpoint intensity +guessOffset= (max(pltInts)- min(pltInts)) * guessOffsetAdj; %Could plug-in a curvefit routine to determine an 'L' value +guessVal = max(pltInts)- guessOffset; +guesTpt= find(pltInts > guessVal,1,'first'); +if isempty(guesTpt) + guesTpt= 1; +end + +%************************************************************************** + PtptIm= tifFileLst(guesTpt); + PtptImg= char(tifFileLst(guesTpt)); + ImageRd= imread(PtptImg); %Ncode actually reading a .bmp file for Ncode + PltImg= (ImageRd(1:ImHeigth,1:ImWidth)); + PltImg= 255 - PltImg; + +%Convert to image to BW then AND with PTmapOnesDbl + GraythrShift= 1; + BWimage= im2bw(PltImg, GraythrShift*graythresh(PltImg)); + matchAnded= BWimage & PTmapOnesDbl; + initPTANDed= matchAnded; + ANDValInitPT=sum(sum(initPTANDed(:,:))); + +%******************************************************************* +%*********Begin Image Registration for Pintool***************************** +%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +MaxANDVal= 0; %sum(sum(BWimage(:,:))); %seed MaxANDVal from the BWimage +regDith= 4; +TestmapOnesDbl= zeros(2075,1400); +%##################################################################### +%Create a Search Range centered on C==1 and R==1 of the PinTool map +RtopBndry= 96; RbotBndry= 176; +ClfBndry = 96; CrtBndry= 154; + [PTrefPt]=PTmapPos{1,1}; + PTrefPtR=PTrefPt(1); + PTrefPtC=PTrefPt(2); + + rngUpper= RtopBndry - PTrefPtR; + quotT= round(abs(rngUpper)/regDith); + topRang= regDith * quotT; + if rngUpper<0, topRang= -1*topRang; end + + rngLower= RbotBndry - PTrefPtR; + quotB= round(abs(rngLower)/regDith); + botRang= regDith * quotB; + if rngLower<0, botRang= -1*botRang; end + + %----------------------------- + rngLeft= ClfBndry - PTrefPtC; + quotLf= round(abs(rngLeft)/regDith); + leftRang= regDith * quotLf; + if rngLeft<0, leftRang= -1*leftRang; end + + rngRight= CrtBndry - PTrefPtC; + quotRt= round(abs(rngRight)/regDith); + rightRang= regDith * quotRt; + if rngRight<0, rightRang= -1*rightRang; end +%######################################################################### +%*************Orthagonal Registration ************************ + testPos = cell(24,16); + for shiftPosR= topRang:regDith:botRang + for shiftPosC=leftRang:regDith:rightRang + TestmapOnesDbl= zeros(2075,1400); + for r=1:24 + for c=1:16 + [PTrefPt]=PTmapPos{r,c}; + testPosR= PTrefPt(1) + shiftPosR; + testPosC= PTrefPt(2) + shiftPosC; + + a=[testPosR,testPosC]; %testPos(r,c)= {testPosR,testPosC}; + testPos(r,c)={a}; + + %******* doCircle + TestmapOnesDbl(testPosR:(testPosR+(diaExt-1)),testPosC:(testPosC+(diaExt-1)))= cirMask; %optCirMask; + + end + end + GraythrShift= 1.0; + BWtestimage= im2bw(TestmapOnesDbl, GraythrShift*graythresh(TestmapOnesDbl)); + matchAnded= BWimage & BWtestimage; + ANDVal=sum(sum(matchAnded(:,:))); + if ANDVal> MaxANDVal + MaxANDVal= ANDVal; + MaxPos= testPos; + MaxPosDbl= TestmapOnesDbl; + end + end + end + if MaxANDVal> ANDValInitPT + MaxPosDbl= MaxPosDbl; + MaxPos= MaxPos; + else + MaxPosDbl= PTmapOnesDbl; + MaxPos= PTmapPos; + end + + +%*************SKEW Registration******************** + testSkPos= MaxPos; + MaxANDValSk= 0; %MaxANDVal; + skQty= 0.05; + regDith=4; + for shiftPos= -regDith/2:regDith/2 + for k= -5:5 + sk= k* skQty; + TestmapSkDbl= zeros(2075,1400); + for r=1:24 + for c=1:16 + [skrefPt]=testSkPos{r,c}; + testSkPosR= round(skrefPt(1) + (c*sk + shiftPos)); + testSkPosC= round(skrefPt(2) + (-r*sk + shiftPos)); + + a=[testSkPosR,testSkPosC]; %testPos(r,c)= {testPosR,testPosC}; + testSkPos(r,c)={a}; + + %******* doCircle + TestmapSkDbl(testSkPosR:(testSkPosR+(diaExt-1)),testSkPosC:(testSkPosC+(diaExt-1)))= cirMask; %optCirMask; + + end + end + GraythrShift= 1; + BWtestSkimage= im2bw(TestmapSkDbl, GraythrShift*graythresh(TestmapSkDbl)); + matchAndedSk= BWimage & BWtestSkimage; + ANDValSk=sum(sum(matchAndedSk(:,:))); + if ANDValSk> MaxANDValSk + MaxANDValSk= ANDValSk; + MaxPosSk= testSkPos; + MaxPosSkDbl= TestmapSkDbl; + end + end + end + + if MaxANDValSk> MaxANDVal + MaxANDVal= MaxANDValSk; + MaxPosDbl= MaxPosSkDbl; + MaxPos= MaxPosSk; + end + +%----------------------------------------- +%*********Center of Area for Final Registration**************************** +%Note GraythrShift in NImapPT can be used to a lower or higher 'decider' +%and thus effect the centroid(center of area) that this module determines. +%GraythrShift= 1.0 normally allowing the Matlab graythresh to calc. the BW decider + +range= 0; +lastDetMaxCentDbl= zeros(2075,1400); +lastDetMaxPos= cell(24,16); +%*************CIRCLE related************************* +for r=1:24 + for c=1:16 + [refPt]=MaxPos{r,c}; + refPtR=refPt(1); + refPtC=refPt(2); + + + %>>>>>>>>>>>>>>>>>>>>>>>Calc. Centroid Position of spot>>>>>>>>>>>>>>>>>>>> + detPrLo= (MaxPos{r,c}(1)) -(range); detPrHi= (MaxPos{r,c}(1)) + (diaExt+range); + detPcLo= (MaxPos{r,c}(2)) -(range); detPcHi= (MaxPos{r,c}(2)) + (diaExt+range); + + rsum=(sum(BWimage(detPrLo:detPrHi,detPcLo:detPcHi),2)); + csum=(sum(BWimage(detPrLo:detPrHi,detPcLo:detPcHi),1)); + if sum(rsum)>1 && sum(csum)>1 + Rvec = detPrLo:detPrHi; + Cvec = detPcLo:detPcHi; + centR= round((sum(Rvec .*rsum')) /sum(rsum)); + centC= round((sum(Cvec .*csum)) /sum(csum)); + aa=[(centR-diaExt/2) (centC-diaExt/2)]; + else + rsum=(sum(MaxPosDbl(detPrLo:detPrHi,detPcLo:detPcHi),2)); + csum=(sum(MaxPosDbl(detPrLo:detPrHi,detPcLo:detPcHi),1)); + Rvec = detPrLo:detPrHi; + Cvec = detPcLo:detPcHi; + centR= round((sum(Rvec .*rsum')) /sum(rsum)); + centC= round((sum(Cvec .*csum)) /sum(csum)); + %} + end + aa=[(centR-diaExt/2) (centC-diaExt/2)]; + lastDetMaxPos(r,c)={aa}; + %**********2nd Recenter Based on Previous ***************************** + detPrLo= lastDetMaxPos{r,c}(1) -(range); detPrHi= lastDetMaxPos{r,c}(1) + (diaExt+range); + detPcLo= lastDetMaxPos{r,c}(2) -(range); detPcHi= lastDetMaxPos{r,c}(2) + (diaExt+range); + + rsum=(sum(BWimage(detPrLo:detPrHi,detPcLo:detPcHi),2)); + csum=(sum(BWimage(detPrLo:detPrHi,detPcLo:detPcHi),1)); + if sum(rsum)>1 && sum(csum)>1 + Rvec = detPrLo:detPrHi; + Cvec = detPcLo:detPcHi; + centR= round((sum(Rvec .*rsum')) /sum(rsum)); + centC= round((sum(Cvec .*csum)) /sum(csum)); + aa=[(centR-diaExt/2) (centC-diaExt/2)]; + else + rsum=(sum(MaxPosDbl(detPrLo:detPrHi,detPcLo:detPcHi),2)); + csum=(sum(MaxPosDbl(detPrLo:detPrHi,detPcLo:detPcHi),1)); + Rvec = detPrLo:detPrHi; + Cvec = detPcLo:detPcHi; + centR= round((sum(Rvec .*rsum')) /sum(rsum)); + centC= round((sum(Cvec .*csum)) /sum(csum)); + end + aa=[(centR-diaExt/2) (centC-diaExt/2)]; + lastDetMaxPos(r,c)={aa}; + + [CentPt]= lastDetMaxPos{r,c}; + CentPtR= CentPt(1); + CentPtC= CentPt(2); + %******* doCircle + lastDetMaxCentDbl(CentPtR:(CentPtR+(diaExt-1)),CentPtC:(CentPtC+(diaExt-1)))= cirMask; %optCirMask; + + %******************************************************* + + end %end c +end %end r, end for search for maximum intensity spot detection alignment +%****************************************** + + BWtestCentimage= im2bw(lastDetMaxCentDbl, GraythrShift*graythresh(lastDetMaxCentDbl)); + matchAndedCent= BWimage & BWtestCentimage; + ANDValCent=sum(sum(matchAndedCent(:,:))); + if ANDValCent>= MaxANDVal + MaxPosDbl= lastDetMaxCentDbl; + MaxPos= lastDetMaxPos; + %detPos=lastDetMaxPos; + end + %detPos=lastDetMaxPos; + detPos=MaxPos; + + + + + + +%conclusion of NIplateImagPTreg******************************************** +%**************************************************************************** + +end %if length(tifFileLst) >2 +%************************************************************************* + + + +%preallocation for speed on 2019_1121 *** +scIntens= zeros(24,16,length(tifFileLst)); +BkgrdMat= zeros(24,16,length(tifFileLst)); +%areaOfIntensAboveBG= zeros(24,16,length(tifFileLst)); +lastDetMaxPos= cell(numRows, numCols); +%*********************************************************************************************************************** +for tPt=1:length(tifFileLst) + tPt; + +tifFile= char(tifFileLst(tPt)); + + +fullsc= plateImage{tPt}; +fIntsc=fullsc; +%**********Bad Image Filter Zone *************************** +%NORMALLY turned off; This is Hand HARD CODE SPECIFIC to particular ExpJOB +specificHardFilterCode=0; +if specificHardFilterCode==1 && tPt==1 + %Plate scratch filter prototype could be added here +fIntsc(:,1105:1190)=uint8(67); %Specific blocking area set to background level here +%fIntsc(:,:)=uint8(67); % set entire image to background level if 1.bmp is +%bad. This allows acquisiton of the initial "Print Time" from .bmp time stamp. +end +%************************************************** +tp=90; bt=1935; lf=90; rt=1318; blkVal=10; %tp=90; +fIntsc(1:tp,:)=uint8(blkVal); %top of blocking frame +fIntsc(bt:ImHeigth,:)=uint8(blkVal); %bottom of blocking frame +fIntsc(:,1:lf)=uint8(blkVal); %left side of blocking frame +fIntsc(:,rt:ImWidth)=uint8(blkVal); %right side of blocking frame +bwIntsc= im2bw(fIntsc,80/255); + +%********************Collect Intensities (Absorbances)************************** +%end for search for maximum intensity spot detection alignment +if doCircle==1 +earlySrc1Thres=.31*cirPixA*255; %Initial search spot det. threshold +earlySrc2Thres= spotThres*cirPixA*255; %Course search spot det. threshold +definedSpot= (spotThres+.03)*cirPixA*255; %typically .34 ; 20000; %Intensity at which a spot is well defined for a width=24 {formerly20} % 0.5*((width^2 * (255- backgrdPixel)) + +else +earlySrc1Thres=.31*(width^2)*255; %Initial search spot det. threshold +earlySrc2Thres= spotThres*(width^2)*255; %Course search spot det. threshold +definedSpot= (spotThres+.03)*(width^2)*255; %typically .34 ; 20000; %Intensity at which a spot is well defined for a width=24 {formerly20} % 0.5*((width^2 * (255- backgrdPixel)) +end +detMaxPos = cell(numRows, numCols); +bcent = cell(numRows, numCols); + +for r=1:numRows + for c=1:numCols + [refPt]=detPos{r,c}; + refPtR=refPt(1); InitR=refPt(1); + refPtC=refPt(2); InitC=refPt(2); + %*****ZeroSearchIntens **************************************************** + %clear ('searchIntens') + searchIntens= []; + + if doCircle==1 + % cutout= fIntsc(refPtR+1:(refPtR+(diaExt)),(refPtC+1:(refPtC+diaExt))); + cutout= fIntsc(refPtR:(refPtR+(diaExt-1)),refPtC:(refPtC+(diaExt-1))); + Imcir= cutout .* uint8(cirMask); + zeroSearchIntens=sum(sum(Imcir(1:(diaExt-1),1:(diaExt-1)))); + + else %standard square analysis area + refPtRExt=refPtR+widthEx; + if (refPtRExt)> size(fIntsc,1), refPtRExt= size(fIntsc,1); end + refPtCExt=refPtC+widthEx; + if (refPtCExt)> size(fIntsc,2), refPtCExt= size(fIntsc,2); end + zeroSearchIntens=sum(sum(fIntsc(refPtR:(refPtRExt),refPtC:(refPtCExt)))); + end + %************************************************************************** + %clear ('searchIntens') + searchIntens= []; + intensMax=-1; + src1cnt=0; %zero count for oneshot calc of centroid after initial location of spot during Src1 + if stopSearch(r,c)==false + + %*******************Set Range and Dither******************** + if intensPrev(r,c)< earlySrc1Thres %First Search parameters insert + rRangeUpper=rCRangeUpper; + rRangeLower=rCRangeLower; + cRangeLower=cCRangeLower; + cRangeUpper=cCRangeUpper; + dither=ditherI; + + [searchStartRef]= detPos{r,c}; + + elseif intensPrev(r,c)> earlySrc1Thres &&intensPrev(r,c)<= earlySrc2Thres %Second Search parameters insert + rRangeUpper= -CSrchRng; %rCRangeUpper; + rRangeLower= CSrchRng; %rCRangeLower; + cRangeLower= -CSrchRng; %cCRangeLower; + cRangeUpper= CSrchRng; %cCRangeUpper; + dither=ditherC; + + %src1cnt=src1cnt+1; + %After spot defined above selected threshold, start search from + %center of the calc'd spot area + %Calc. Center of Area (black and white (1's and 0's) + %detPrLo= (detPos{r,c}(1)) -18; detPrHi= (detPos{r,c}(1)) +width+18; + %detPcLo= (detPos{r,c}(2)) -18; detPcHi= (detPos{r,c}(2)) +width+18; + detPrLo= (detPos{r,c}(1)) -CSrchRng; detPrHi= (detPos{r,c}(1)) +width+CSrchRng; %Mod 15_0707 User contrain search + detPcLo= (detPos{r,c}(2)) -CSrchRng; detPcHi= (detPos{r,c}(2)) +width+CSrchRng; %Mod 15_0707 User contrain search + + brsum=(sum(bwIntsc(detPrLo:detPrHi,detPcLo:detPcHi),2)); + bcsum=(sum(bwIntsc(detPrLo:detPrHi,detPcLo:detPcHi),1)); + bRvec= detPrLo:detPrHi; + bCvec= detPcLo:detPcHi; + bcentR= round((sum(bRvec .*brsum')) /sum(brsum)); + bcentC= round((sum(bCvec .*bcsum)) /sum(bcsum)); + bcent{r,c}= [bcentR, bcentC]; + cornerOfset= ceil((0.5*width)+CSrchRng); %ceil(0.5*(width+(2*CSrchRng))); + refCorner{r,c}= [(bcentR-cornerOfset),(bcentC-cornerOfset)]; + + %NAN blowup preventive + if sum(brsum)>5 && sum(bcsum)>5 + %[searchStartRef]= cell2mat(bcent(r,c)); + [searchStartRef]= cell2mat(refCorner(r,c)); + else + [searchStartRef]= detPos{r,c}; + end + + elseif intensPrev(r,c) >earlySrc2Thres && intensPrev(r,c)<= definedSpot + rRangeUpper=rFRangeUpper; %Fine Search parameters insert + rRangeLower=rFRangeLower; + cRangeLower=cFRangeLower; + cRangeUpper=cFRangeUpper; + dither=ditherF; + zz=lastDetMaxPos(r,c); + [searchStartRef]= zz{1,1}; + end + %*******************Fixed to Roam Switch******************** + %{ + if intensPrev(r,c)<=earlySrc1Thres + % [searchRef]= detPos{r,c}+ [rOff,cOff]; + [searchStartRef]= detPos{r,c}; + end + + if intensPrev(r,c)> earlySrc1Thres && intensMax= earlySrc2Thres, %Let Roam + zz=lastDetMaxPos(r,c); + [searchStartRef]= zz{1,1}; + + end + %} + %************************************************************* + extent=CSrchRng; + Rmin=InitR-extent; Rmax=InitR+extent; + Cmin=InitC-extent; Cmax=InitC+extent; + + %clear cutout %just found 20_0327 + cutout= []; + for rOff=rRangeUpper:dither:rRangeLower + for cOff=cRangeLower:dither:cRangeUpper + [searchRef]= searchStartRef + [rOff,cOff]; + refPtR=searchRef(1); + refPtC=searchRef(2); + if refPtRRmax,refPtR= Rmax;end + if refPtCCmax,refPtC= Cmax;end + %******Do Circle analysis****** + if doCircle==1 + % cutout= fIntsc(refPtR+1:(refPtR+(diaExt)),(refPtC+1:(refPtC+diaExt))); + try + cutout= fIntsc(refPtR:(refPtR+(diaExt-1)),refPtC:(refPtC+(diaExt-1))); + catch + tPt; + r; + c; + end + Imcir= cutout .* uint8(cirMask); + searchIntens=sum(sum(Imcir(1:(diaExt-1),1:(diaExt-1)))); + + else %standard square analysis area + refPtRExt=refPtR+widthEx; + if (refPtRExt)> size(fIntsc,1), refPtRExt= size(fIntsc,1); end + refPtCExt=refPtC+widthEx; + if (refPtCExt)> size(fIntsc,2), refPtCExt= size(fIntsc,2); end + searchIntens=sum(sum(fIntsc(refPtR:(refPtRExt),refPtC:(refPtCExt)))); + end + %**************************************** + if searchIntens > intensMax + intensMax=searchIntens; + intens(r,c)=searchIntens; + intensPrev(r,c)=searchIntens; + a=[refPtR,refPtC]; + detMaxPos(r,c)={a}; + lastDetMaxPos(r,c)=detMaxPos(r,c); + end + if intensMax>definedSpot || intensMax<2 + stopSearch(r,c)=true; + end + end %for cOff + end %for rOff + + if intensMax< noSpot + intensMax= zeroSearchIntens; + lastDetMaxPos(r,c)= {refPt}; + end + + else %elseif Spot is well defined STOPSearch strong + [refPt]=lastDetMaxPos{r,c}; + refPtR=refPt(1); + refPtC=refPt(2); + if doCircle==1 + cutout= fIntsc(refPtR:(refPtR+(diaExt-1)),refPtC:(refPtC+(diaExt-1))); + Imcir= cutout .* uint8(cirMask); + intens(r,c)=sum(sum(Imcir(1:(diaExt-1),1:(diaExt-1)))); + else %standard square analysis area + intens(r,c)=sum(sum(fIntsc(refPtR:(refPtR+widthEx),refPtC:(refPtC+widthEx)))); + end + end %stopSearch loop + %******************************************************* + %{ + if intensMax> earlySrc1Thres && intensMax< earlySrc2Thres && src1cnt==0, + + %{ + src1cnt=src1cnt+1; + detPrLo= (detPos{r,c}(1)) -12; detPrHi= (detPos{r,c}(1)) +width+12; + detPcLo= (detPos{r,c}(2)) -12; detPcHi= (detPos{r,c}(2)) +width+12; + rsum=(sum(fIntsc(detPrLo:detPrHi,detPcLo:detPcHi),2)); + csum=(sum(fIntsc(detPrLo:detPrHi,detPcLo:detPcHi),1)); + Rvec= detPrLo:detPrHi; + Cvec= detPcLo:detPcHi; + centR= round((sum(Rvec .*rsum')) /sum(rsum)); + centC= round((sum(Cvec .*csum)) /sum(csum)); + cent{r,c}= [centR, centC]; + %} + end + %} + %******************************************************* + + + end %for c +end %for r -end for search +%**************************************************************************************************************** +%[F_spots,bmtp, BGthres,bmm,rwm,lstTpt,BGsc,totBkgrd]= ... + % NIgenBkGrdDataPar4Fnc(ImParMat, tifFileLst, lastDetMaxPos, numRows, numCols, fullsc, tPt,BGthres, doCircle, cirPixA,diaExt, cirMask, BGTav); +%NIgenBkGrdDataPar4Fnc + +%function [F_spots,bmtp, BGthres,bmm,rwm,lstTpt,BGsc,totBkgrd ]= ... +% NIgenBkGrdDataPar4Fnc(ImParMat, tifFileLst, lastDetMaxPos, numRows, numCols, fullsc, tPt, BGthres, doCircle, cirPixA, diaExt, cirMask, BGTav)% Called by NIscanIntensBG.m +% +%global resDir +global ExpOutmat +%global pixsAboveBG +MPnum=ImParMat(1); +destPerMP=ImParMat(2); +spotThres=ImParMat(4)/100; %selScan=ImParMat(4); +width=ImParMat(5); +dither=ImParMat(6); +BGthresInput=ImParMat(3); % a value of percent to be added (20=>1.20*BGTav) +lstTpt=length(tifFileLst); + +LfOffset= 65; %lf; %Ncode +TopOffset= 65; %tp; +rtLim=1350; %1400; %1350; +rtNotchOffset=10; +cushion=20; + %close2edgeTol=90; +botLim= 2003; +botNotchOffset=60; + + +%**********Initialize************************** +%clear('totBkgrd') +%clear('pixsAboveBG') +%clear('lineV','pixCnt','lineSum,meanBkgrd') +%clear('rw','rwm') +totBkgrd= 0; +pixsAboveBG= 0; +lineV= 0; pixCnt= 0; lineSum= 0; meanBkgrd= 85; +rw= 0; rwm= 0; +%************************************************* +%****Start Generation of BackGround data after the 2nd Row (&& 2nd Col) +%*****Optomized spot positions determined. *********** + +%if r>=2 && c>=2, + BGsc=zeros(2075,1400); %Empsc; + widthEx=width-1; %width extention from reference point +for r=1:numRows + for c=1:numCols + refP=lastDetMaxPos{r,c}; + if r1, refPrR=lastDetMaxPos{r-1,c}; end + if c>1, refPrC=lastDetMaxPos{r,c-1}; end + %****At least NINE conditions which must be treated differently + % ++++++++++ Condition (0) ---most spots fall in this central category + if r>1 && r1 && c0 + refBG(1)=refP(1)- floor(0.5*(refP(1) - TopOffset)); + refBG(2)=refP(2)- floor(0.5*(refP(2) - LfOffset)); %Icode -0 + refBG(3)=refP(1)+widthEx+ floor(0.5*((refNxR(1)-(refP(1)+widthEx)))); + refBG(4)=refP(2)+widthEx+ floor(0.5*((refNxC(2)-(refP(2)+widthEx)))); + detBG(r,c)= {uint8(refBG)}; %Store Background Reference pt data in a Cell + +%clear('lineV','pixCnt','lineSum') +lineV= fullsc(refBG(1),refBG(2):refBG(4)-1); %mean across top rt +if tPt==1 || tPt==lstTpt,rw=lineV';end +pixCnt(1)=nnz(lineV .* uint8(lineV79 + if r==1 && c>1 && c(size(fullsc,2)-close2edgeTol)%if intens area too close to edge + if refBGcExt>(rtLim-cushion)%if intens area too close to edge + refBG(4)=rtLim-rtNotchOffset; + else + refBG(4)=refP(2)+widthEx+ floor(0.5*(rtLim-refBGcExt)); + end + detBG(r,c)= {uint8(refBG)}; %Store Background Reference pt data in a Cell + +%clear('lineV','pixCnt','lineSum') +lineV= fullsc(refBG(1),refBG(2):refBG(4)-1); % across top rt +if tPt==1 || tPt==lstTpt,rw=lineV';end +pixCnt(1)=nnz(lineV .* uint8(lineV(size(fullsc,2)-close2edgeTol) %if intens area too close to edge +if refBGcExt>(rtLim-cushion) %if intens area too close to edge + meanBkgrd = (sum(lineSum(1:4)))./ (sum(pixCnt(1:4))); + if sum(pixCnt(1:4))==0, meanBkgrd=double(BGthres(r,c));end + %meanBkgrd = (lineSum(1)+lineSum(2)+lineSum(4))./ ... + %((2*(refBG(4)-refBG(2))+(refBG(3)-refBG(1)))); %divby 2*across+1*down + if doCircle==1, totBkgrd(r,c) = meanBkgrd* cirPixA; else totBkgrd(r,c) = meanBkgrd* width^2; end + +else + %clear('lineV') %3 + lineV = fullsc(refBG(1):refBG(3)-1,refBG(4)); %mean down rt + if tPt==1 || tPt==lstTpt,rw=[rw;lineV];end + pixCnt(3)=nnz(lineV .* uint8(lineV47 farRight col edge + if r>1 && r(botLim-botNotchOffset)%if intens area not too close to bot edge else + refBG(3)=botLim + cushion; + else + refBG(3)=refP(1)+widthEx+ floor(0.5*(botLim-refBGrExt)); + end + if refBGcExt>(rtLim- cushion)%if intens area not too close to edge else + refBG(4)= rtLim -cushion; + else + refBG(4)=refP(2)+widthEx+ floor(0.5*(rtLim-refBGcExt)); + end + detBG(r,c)= {uint8(refBG)}; %Store Background Reference pt data in a Cell + +%clear('lineV','pixCnt','lineSum') %1 +lineV= fullsc(refBG(1),refBG(2):refBG(4)-1); % across top rt +if tPt==1 || tPt==lstTpt,rw=lineV';end +pixCnt(1)=nnz(lineV .* uint8(lineV(size(fullsc,2)-close2edgeTol)&&... + if refBGcExt>(rtLim-cushion)&&... + refBGrExt<(botLim-botNotchOffset) %if intens area too close to right edge only + %clear('lineV') %4 + lineV = fullsc(refBG(3),refBG(2):refBG(4)-1); %mean across bot rt + if tPt==1 || tPt==lstTpt,rw=[rw;lineV'];end + pixCnt(4)=nnz(lineV .* uint8(lineV(botLim-botNotchOffset) &&... + refBGcExt>(rtLim-cushion)%if intens area too close to bottom edge only + %clear('lineV') %3 + lineV = fullsc(refBG(1):refBG(3)-1,refBG(4)); %mean down rt + if tPt==1 || tPt==lstTpt,rw=[rw;lineV];end + pixCnt(3)=nnz(lineV .* uint8(lineV(size(fullsc,2)-close2edgeTol)&& ... +elseif refBGcExt>(rtLim-cushion)&& ... + refBGrExt>(botLim-botNotchOffset) %if intens area too close to both right edge and bottom edge + meanBkgrd = (sum(lineSum(1:4)))./ (sum(pixCnt(1:4))); + if sum(pixCnt(1:4))==0, meanBkgrd=double(BGthres(r,c));end + %meanBkgrd = (lineSum(1)+lineSum(2))./ ... + %((refBG(4)-refBG(2)+(refBG(3)-refBG(1)))); %divby 1across+1*down + if doCircle==1, totBkgrd(r,c) = meanBkgrd* cirPixA; else totBkgrd(r,c) = meanBkgrd* width^2; end + +else + %clear('lineV') %3 + lineV = fullsc(refBG(1):refBG(3)-1,refBG(4)); %mean down rt + if tPt==1 || tPt==lstTpt,rw=[rw;lineV];end + pixCnt(3)=nnz(lineV .* uint8(lineV1 && c(botLim-botNotchOffset)%if intens area not too close to bot edge else + refBG(3)=botLim + cushion; + else + refBG(3)=refP(1)+widthEx+ floor(0.5*(botLim-refBGrExt)); + end + if rem(c,numCols)==0 + refBG(4)=refP(2)+widthEx+ floor(0.5*(plateEdgeSpace)); + else + refBG(4)=refP(2)+widthEx+ floor(0.5*(refNxC(2)-(refP(2)+widthEx))); + %refBG(4)=refP(2)+widthEx+ floor(0.5*(size(fullsc,2)-refBGcExt)); + end + + detBG(r,c)= {uint8(refBG)}; %Store Background Reference pt data in a Cell + +%clear('lineV','pixCnt','lineSum') %1 +lineV= fullsc(refBG(1),refBG(2):refBG(4)-1); % across top rt +if tPt==1 || tPt==lstTpt,rw=lineV';end +pixCnt(1)=nnz(lineV .* uint8(lineV(botLim-botNotchOffset) %if intens area too close to bottom edge + meanBkgrd = (sum(lineSum(1:4)))./ (sum(pixCnt(1:4))); + if sum(pixCnt(1:4))==0, meanBkgrd=double(BGthres(r,c));end + %meanBkgrd = (lineSum(1)+lineSum(2)+lineSum(4))./ ... + %((refBG(4)-refBG(2)+(2*(refBG(3)-refBG(1))))); %divby 1*across+2*down + if doCircle==1, totBkgrd(r,c) = meanBkgrd* cirPixA; else totBkgrd(r,c) = meanBkgrd* width^2; end + +else + %clear('lineV') %4 + lineV = fullsc(refBG(3),refBG(2):refBG(4)-1); %mean across bot rt + if tPt==1 || tPt==lstTpt,rw=[rw;lineV'];end + pixCnt(4)=nnz(lineV .* uint8(lineV(botLim-botNotchOffset)%if intens area not too close to bot edge else + refBG(3)=botLim + cushion; + else + refBG(3)=refP(1)+widthEx+ floor(0.5*(botLim-refBGrExt)); + end + refBG(4)=refP(2)+widthEx+ floor(0.5*(refNxC(2)-(refP(2)+widthEx))); + + detBG(r,c)= {uint8(refBG)}; %Store Background Reference pt data in a Cell + +%clear('lineV','pixCnt','lineSum') %1 +lineV= fullsc(refBG(1),refBG(2):refBG(4)-1); % across top rt +if tPt==1 || tPt==lstTpt,rw=lineV';end +pixCnt(1)=nnz(lineV .* uint8(lineV(botLim-botNotchOffset) %if intens area too close to bottom edge + meanBkgrd = (sum(lineSum(1:4)))./ (sum(pixCnt(1:4))); + if sum(pixCnt(1:4))==0, meanBkgrd=double(BGthres(r,c));end + %meanBkgrd = (lineSum(1)+lineSum(2)+lineSum(3))./ ... + %((refBG(4)-refBG(2)+(2*(refBG(3)-refBG(1))))); %divby 1*across+2*down + if doCircle==1, totBkgrd(r,c) = meanBkgrd* cirPixA; else totBkgrd(r,c) = meanBkgrd* width^2; end + + else +%clear('lineV') %4 +lineV = fullsc(refBG(3),refBG(2):refBG(4)-1); %mean across bot rt +if tPt==1 || tPt==lstTpt,rw=[rw;lineV'];end +pixCnt(4)=nnz(lineV .* uint8(lineV1 && rBG); %calc non-background within circular analysis area + %pixsAboveBG(r,c)=nnz(fullsc(refPtRx:(refPtRx+ (diaExt-1)),refPtCx:(refPtCx+ (diaExt-1)))>BG); %BG); + else %standard square analysis area + pixsAboveBG(r,c)=nnz(fullsc(refPtRx:(refPtRx+widthEx),refPtCx:(refPtCx+widthEx))>BG); %calc non-background within square analysis area + end + +%--------------------------------------------------------- +%TimePT average +%meanBGTsum(r,c)=meanBGTsum(r,c)+meanBkgrd; +%meanBGTav(r,c)= meanBGTsum(r,c) /tPt; +bmm(r,c)=uint8(floor(meanBkgrd)); + +% For FOTOS only FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF +F_spots(r,c,tPt)={refBG}; %Foto spot coordinates +bmtp(r,c,tPt)=bmm(r,c); %Added/used for FOTOS 11_0830 +%FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF + +%TimePT Average BG intensity calc method +BGTav(r,c)=uint8(round((BGTav(r,c)*(tPt-1))+meanBkgrd)/tPt); +%BGTav(r,c)=((BGTav(r,c)*(tPt-1))+meanBkgrd)/tPt; + +%sliding BG average +%{ + if r==1 && c==1, slavg(1:5)= BG;end %slavg(1:5)= 50; %120620 Epsonvalue is 33; end + slavg(5)=slavg(4); + slavg(4)=slavg(3); + slavg(3)=slavg(2); + slavg(2)=slavg(1); + slavg(1)=BG; + slidingAvg=mean(slavg); +%} +%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +% For Visualization Purposes only +%1->2 across rt 1, 2-4 +BGsc(refBG(1),refBG(2):refBG(4))=1; +%1->4 down lf 1-3, 2 +BGsc(refBG(1):refBG(3),refBG(2))=1; +%2->3 down rt 1-3, 4 +BGsc(refBG(1):refBG(3),refBG(4))=1; +%4->3 Bot Across rt 1, 2-4 +BGsc(refBG(3),refBG(2):refBG(4))=1; +%BGsc(refBG(1):refBG(3),refBG(2):refBG(4)) %just displays matrices + + + end % Outer c (column loop) + +end % outer r (row loop) +BGthres=uint8(double(BGTav) * (1+(BGthresInput/100))); + + + +%end % Outer if r>1 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%END OF NIgenBkGrdData brought inside for GLOBAL PARALLEL +%***************************************************************************************************************** +%***************************************************************************************************************** +%***************************************************************************************************************** +%***************************************************************************************************************** +%***************************************************************************************************************** +optomizedPos(:,:,tPt)=lastDetMaxPos; +%if(isequal(SWsingleSc,1))|| ((~isequal(SWsingleSc,1))&&tPt==tptLength) %190731 change to capture last timePt PTmap if(isequal(SWsingleSc,1)) +%if(isequal(SWsingleSc,1))&& (tPt==1 || tPt== guesTpt || tPt==tptLength)|| ((~isequal(SWsingleSc,1))&&tPt==tptLength) %190919 Option for observation +%if (tPt==1 || tPt== guesTpt || tPt==tptLength) %190919 Option for observation +if (tPt==tptLength) %190919 change to capture ONLY last timePt PTmap +%*******************OptmapOnes Empsc; +%***********vis************************* +%clear OptmapOnes; +OptmapOnes= []; +%OptmapOnes=Empsc+ 1.4; +OptmapOnesDbl=double(Empsc)+ .8; +for r=1:numRows + for c=1:numCols + [refPt]=lastDetMaxPos{r,c}; + refPtR=refPt(1); + refPtC=refPt(2); + + if doCircle==1 + %OptmapOnesDbl(refPtR:(refPtR+diaExt),refPtC:(refPtC+diaExt))= cirMask; + OptmapOnesDbl(refPtR:(refPtR+(diaExt-1)),refPtC:(refPtC+(diaExt-1)))=optCirMask; + + + else + refPtRExt=refPtR+widthEx; + if (refPtRExt)> size(fullsc,1), refPtRExt= size(fullsc,1); end + refPtCExt=refPtC+widthEx; + if (refPtCExt)> size(fullsc,2), refPtCExt= size(fullsc,2); end + OptmapOnesDbl(refPtR:(refPtRExt),refPtC:(refPtCExt))= 1;%255; + end + end +end +%****Cross Optomized Map Matrix onto tifFile Matrix*************** +%***********for visualization purposes only!!!!!!!!!!!!!!! + + if tPt==1,figure;end + + resIm= uint8((double(BGsc) + OptmapOnesDbl) .* double(fullsc)); + %clear('OptmapOnesDbl') + OptmapOnesDbl= []; + %hres=imagesc(resIm); + hfscanIm=imagesc(resIm); + hfIm=gcf; + haxis=gca; + title(strcat('Scan',num2str(selScan),'timePt-',num2str(tPt))); + clf(hfscanIm,'reset'); + set(hfIm,'NumberTitle','off') + set(hfIm,'Name', char(ExpOutmat)) %strcat('EASYconsole- ',char(resDir))) + if Fflg==1 + resDir; %for debugging + fullfile(resDir,'Fotos',strcat('Scan',num2str(selScan),'_timePt-',num2str(tPt))) + hgsave(fullfile(resDir,'Fotos',strcat('Scan',num2str(selScan),'_timePt-',num2str(tPt)))); %F 14_0626 + else + hgsave(fullfile(resDir,'figs',strcat('Scan',num2str(selScan),'_timePt-',num2str(tPt)))); + end + +else + selScan; + %tifFile +end +%****************************************************************** + +BkgrdMat=totBkgrd; %BkgrdMat(:,:,tPt)=totBkgrd; %debug parfor effort + +areaOfIntensAboveBG(:,:,tPt)= pixsAboveBG; %actual area of growth for each timepoint (used for printout of grArea ?not used in previous analysis) +scanIntens(:,:,tPt)= intens-BkgrdMat; + + +%capture first and last background values +if tPt==1 + bdg(:,:,1)= bmm; + bdg(:,:,3)=rwm; + +Tmpsbdg1={bdg}; +end +if tPt==lstTpt + bdg(:,:,2)= bmm; + bdg(:,:,4)=rwm; + +Tmpsbdg1={bdg}; +end + +%-------------------------------------------------------------------- + +TmpexpScanIntens1= {scanIntens}; %par4Gbl 20_0122 +TmpFexpScanSpots1= {F_spots}; +TmpFexpScanBMtp1= {bmtp}; +TmpanlZoneRefs1= {optomizedPos}; %added to store tPt by tPt positions + +TmpexpScanIntens2= TmpexpScanIntens1; +TmpFexpScanSpots2= TmpFexpScanSpots1; +TmpFexpScanBMtp2= TmpFexpScanBMtp1; +TmpanlZoneRefs2= TmpanlZoneRefs1; + +Tmpsbdg2= Tmpsbdg1; + +%----------------------------------------------------------- +end % end for tifFile for-loop (i.e. timepoint scan)******************************************************************** + + + \ No newline at end of file diff --git a/workflow/templates/easy/PTmats/NCFparms.mat b/workflow/templates/easy/PTmats/NCFparms.mat new file mode 100755 index 00000000..99dd9927 Binary files /dev/null and b/workflow/templates/easy/PTmats/NCFparms.mat differ diff --git a/workflow/templates/easy/PTmats/NImParameters.mat b/workflow/templates/easy/PTmats/NImParameters.mat new file mode 100755 index 00000000..b258fe4d Binary files /dev/null and b/workflow/templates/easy/PTmats/NImParameters.mat differ diff --git a/workflow/templates/easy/PTmats/NPTdirectParameters.mat b/workflow/templates/easy/PTmats/NPTdirectParameters.mat new file mode 100755 index 00000000..3ec8af05 Binary files /dev/null and b/workflow/templates/easy/PTmats/NPTdirectParameters.mat differ diff --git a/workflow/templates/easy/PTmats/NPTmapDirect.mat b/workflow/templates/easy/PTmats/NPTmapDirect.mat new file mode 100755 index 00000000..aae2e043 Binary files /dev/null and b/workflow/templates/easy/PTmats/NPTmapDirect.mat differ diff --git a/workflow/templates/easy/PTmats/NPTmapSearch.mat b/workflow/templates/easy/PTmats/NPTmapSearch.mat new file mode 100755 index 00000000..0ae7b0ad Binary files /dev/null and b/workflow/templates/easy/PTmats/NPTmapSearch.mat differ diff --git a/workflow/templates/easy/PTmats/NPTsearchParameters.mat b/workflow/templates/easy/PTmats/NPTsearchParameters.mat new file mode 100755 index 00000000..5435fdc7 Binary files /dev/null and b/workflow/templates/easy/PTmats/NPTsearchParameters.mat differ diff --git a/workflow/templates/easy/PTmats/Nbdg.mat b/workflow/templates/easy/PTmats/Nbdg.mat new file mode 100755 index 00000000..92c90c8a Binary files /dev/null and b/workflow/templates/easy/PTmats/Nbdg.mat differ diff --git a/workflow/templates/easy/datatipp.m b/workflow/templates/easy/datatipp.m new file mode 100755 index 00000000..da5afa70 --- /dev/null +++ b/workflow/templates/easy/datatipp.m @@ -0,0 +1,15 @@ +%% PART OF GUI FUNCTIONALITY %% +function output_txt = datatipp(~,event_obj) +% Display the position of the data cursor +% obj Currently not used (empty) +% event_obj Handle to event object +% output_txt Data cursor text string (string or cell array of strings). + +pos = get(event_obj,'Position'); +output_txt = {['X: ',num2str(pos(1),4)],... + ['Y: ',num2str(pos(2),4)]}; + +% If there is a Z-coordinate in the position, display it as well +if length(pos) > 2 + output_txt{end+1} = ['Z: ',num2str(pos(3),4)]; +end diff --git a/workflow/templates/easy/figs/NPTdirect.fig b/workflow/templates/easy/figs/NPTdirect.fig new file mode 100755 index 00000000..ec38838f Binary files /dev/null and b/workflow/templates/easy/figs/NPTdirect.fig differ diff --git a/workflow/templates/easy/figs/searchNPTIm.fig b/workflow/templates/easy/figs/searchNPTIm.fig new file mode 100755 index 00000000..c9d6d97e Binary files /dev/null and b/workflow/templates/easy/figs/searchNPTIm.fig differ diff --git a/workflow/templates/easy/p4loop8c.m b/workflow/templates/easy/p4loop8c.m new file mode 100755 index 00000000..32a02853 --- /dev/null +++ b/workflow/templates/easy/p4loop8c.m @@ -0,0 +1,43 @@ +%% CALLED BY par4Gbl_Main8c.m %% +function[p4L4,... + TmpexpScanIntens4,TmpFexpScanSpots4,TmpFexpScanBMtp4,TmpanlZoneRefs4,Tmpsbdg4]= ... + p4loop8c(parMat,tptLength,numScans,selScanNumLst,SWsingleSc,Fflg,PTmapPos,optCirMask,diaExt,doCircle,cirPixA,cirMask,width, ... + TmpexpScanIntens00,TmpFexpScanSpots00,TmpFexpScanBMtp00,TmpanlZoneRefs00,~,tifFileLstP4,pathname,ImParMat, ... + numRows,numCols,scLst,resDir,expDir, p4L00,TmpexpScanIntens4,TmpFexpScanSpots4,TmpFexpScanBMtp4,TmpanlZoneRefs4, Tmpsbdg00, Tmpsbdg4) +p4L0=p4L00; +TmpexpScanIntens0= TmpexpScanIntens00; +TmpFexpScanSpots0= TmpFexpScanSpots00; +TmpFexpScanBMtp0= TmpFexpScanBMtp00; +TmpanlZoneRefs0= TmpanlZoneRefs00; +Tmpsbdg0= Tmpsbdg00; +%resDir %for debugging +%START PARFOR PARFOR PARFOR PARFOR PARFOR PARFOR PARFOR PARFOR PARFOR PARFOR +%************************************************************************************************************************ +%---------------Global Parfor Loop----------------------------------------- +%for scCount=1:numScans +if SWsingleSc==1 + parforArg=0; +else + parforArg= inf; +end +%for (scCount=1:numScans) +parfor (scCount=1:numScans,parforArg) + scCount %for debugging +p4L0= cell(18,1); +[p4L3,... + TmpexpScanIntens3,TmpFexpScanSpots3,TmpFexpScanBMtp3,TmpanlZoneRefs3,Tmpsbdg3]= ... + par4GblFnc8c(parMat,tptLength,numScans,selScanNumLst,SWsingleSc,Fflg,PTmapPos,optCirMask,diaExt,doCircle,cirPixA,cirMask,width, ... + TmpexpScanIntens0,TmpFexpScanSpots0,TmpFexpScanBMtp0,TmpanlZoneRefs0,scCount,tifFileLstP4,pathname,ImParMat, ... + numRows,numCols, scLst,resDir,expDir, p4L0,Tmpsbdg0); + + p4L4(:,scCount)= p4L3; %(:,scCount); + TmpexpScanIntens4(scCount)= TmpexpScanIntens3; + TmpFexpScanSpots4(scCount)= TmpFexpScanSpots3; + TmpFexpScanBMtp4(scCount)= TmpFexpScanBMtp3; + TmpanlZoneRefs4(scCount)= TmpanlZoneRefs3; + +Tmpsbdg4(scCount)= Tmpsbdg3; + %End of Global Scan loop******************************************************************************************* +%************************************************************************************************************************** +end +%---------------END OF Global Parfor Loop-----------------------------------------%************************************************************************************************************************** \ No newline at end of file diff --git a/workflow/templates/easy/par4GblFnc8c.m b/workflow/templates/easy/par4GblFnc8c.m new file mode 100755 index 00000000..074ab179 --- /dev/null +++ b/workflow/templates/easy/par4GblFnc8c.m @@ -0,0 +1,357 @@ +%% CALLED By p4loop8c.m %% +function [p4L2, ... + TmpexpScanIntens3,TmpFexpScanSpots3,TmpFexpScanBMtp3,TmpanlZoneRefs3,Tmpsbdg3]= ... + par4GblFnc8c(parMat,tptLength,~,selScanNumLst,~,Fflg,PTmapPos,optCirMask,diaExt,doCircle,cirPixA,cirMask,~,... + TmpexpScanIntens,TmpFexpScanSpots,TmpFexpScanBMtp,TmpanlZoneRefs,scCount,tifFileLstP4,~,ImParMat, ... + numRows,numCols,scLst,resDir,expDir,~, Tmpsbdg) + + + +selScan= selScanNumLst(scCount); + + tptLength0= length(tifFileLstP4{scCount}); + %tptLength= numFiles; + tifFileLst={tifFileLstP4(scCount)}; + + % Extract the Imaging time stamps from selected tif files + %clear('e','f'); %can't use clear inside a parfor loop. Preallocation + %can be larger than useable .bmp files! Therefore must be small and + %increased during for loop to maintain cell size integrity between e & f + e= cell(1,2); f= cell(1,2); %(tptLength,2); f= cell(tptLength,2); + nndx=0; + + + for tPt=1:tptLength0 %size(tifFileLst,1) %length(tifFileLst) + tifFileLstP4{scCount}(tPt) + scLst; + scLst(scCount) + char(scLst(scCount)) + char(fullfile(expDir,char(scLst(scCount)))) + expDir; + swCatch=0; + nndx=nndx+1; + +tifFile=char(fullfile(expDir,char(scLst(scCount)), tifFileLstP4{scCount}(tPt))); + try + info = imfinfo(tifFile); %('D:\jwrDevel\DevelCurveFittingJWR\ImageScans\Scan2\020hr002.tif') + + catch ME + %Note: During parallel operation Matlab will not open file (fid) + %etc. Therefore error message text will not be written.The only way + %to get it out would be to pass variable through called function to + %the ..Main8c script outside the parfor loop and print to file from + %there. Consequently for now it only prints error to file when one + %edits p4loop8c.m from 'parfor' to ordinary 'for' loop + fFail=((fullfile(resDir,'PrintResults','ImageFileReadFailure.txt'))); + fid = fopen(fFail,'a'); %create,open and append + %fprintf(fid,'%s \n',char(tifFile)) + fclose(fid) + nndx=nndx-1; + swCatch=1; + rep = getReport(ME, 'basic'); + rep=strcat('Read info failure for-',tifFile,' -', rep); + %fprintf(fid,'%s \n',rep) %See Note: + + end + tptLength= nndx; + scTmNumeric= 1; %initialize for parfor + if swCatch==0 + scTmNumeric(nndx)=datenum(info.FileModDate); + e(nndx,:)= {tifFile, scTmNumeric(nndx)}; + %newtifFileLst(nndx)= tifFileLst(tPt); + end + end % end for tPt=1:length(tifFileLst) loop + + %clear tifFileLst; %Can't use clear inside a parfor loop + tifFileLst= cell(nndx,1); + f= {sortrows(e,2)}; + tifFileLst= f{1,1}(:,1); + + areaOfIntensAboveBG= zeros(24,16,length(tifFileLst)); +% *******Calculate Time series for each Plate of Selected Scan ******** +lastPlateOnScan=1; %Ncode + + tSeriesv=[]; t0Seriesv=[]; + if tptLength>0 %added to jump over and fill data for invalid Sscan(plate runs + scTmNumv= cell2mat(f{1,1}(:,2)); %(:,2)) %120613 fix for variant length data + prtTmNumv=min(scTmNumv)-.000001; + tSeriesv=((scTmNumv-prtTmNumv)*24); + t0Seriesv=((scTmNumv-scTmNumv(1))*24); + + end + + %********************* + +%------------------------------------------------------- + +if tptLength>= 3 %added to jump over and fill data for invalid Sscan(plate runs +%*******Create Blank Scan********************* + ImHeigth=2075; + ImWidth=1400; + Empsc=zeros(ImHeigth,ImWidth); %Ncode + + +%************************************************************************** +%*********Start Scan Loop************************************************** + lastPlateOnScan=1; %Ncode +%********************************************************************************************************************************** +disp('Before call to NIscanIntens.....') +resDir; %debugging parfor +%>>>*******Execute Image conversion into Intensity Data****** +[Tmpsbdg2, scanIntens, ~, ~, ~, TmpexpScanIntens2, TmpFexpScanSpots2, TmpFexpScanBMtp2, TmpanlZoneRefs2,areaOfIntensAboveBG]= ... + NIscanIntensBGpar4GblFnc(Fflg,tifFileLst, ImParMat, PTmapPos,optCirMask,diaExt,doCircle,cirPixA,numRows,numCols,ImHeigth,ImWidth,cirMask, ... + tptLength,selScan,Empsc,TmpexpScanIntens,TmpFexpScanSpots,TmpFexpScanBMtp,TmpanlZoneRefs,resDir, Tmpsbdg); + +%<<<<******************************************************************************************************************************* + +TmpexpScanIntens3= TmpexpScanIntens2; +TmpFexpScanSpots3= TmpFexpScanSpots2; +TmpFexpScanBMtp3= TmpFexpScanBMtp2; +TmpanlZoneRefs3= TmpanlZoneRefs2; +Tmpsbdg3= Tmpsbdg2; + + +%*********************************************************************************** + +%Package Image data ->for Legacy Printing + +%clear plate %Can't use clear inside a parfor loop +CFscanIntens= zeros(384,1);%zeros(16,24); +plate= []; +plate= zeros(24,16,tptLength); + +plate(:,:,:,1)= scanIntens(1:1:24,16:-1:1,:); %TmpexpScanIntens2(1:1:24,16:-1:1,:);%Ncode Dev3Vertical Similar to below + +%******************************* +%NIcheck %Check for bad scans at time points (swapped plates etc.) +%******************************* +%set up cell arrays for storing each plate in each Scan ..(scan,plate) +SWprint=0; +if SWprint==1 +wkDir=pwd; +cd (fullfile(resDir,'PrintResults')); +end + +% Construct Legacy ...Intens.txt file +if SWprint==1 +filename=(strcat('Plate',num2str(selScan),'_Intens.txt')); +fid = fopen(filename,'w'); +end + if length(size(plate))==2 %Only two dims if only one image + numOfImages=1; + else + numOfImages= size(plate,3); + end +CFscanIntens= zeros(384,numOfImages); +locIndx=0; pl=1; +for n=1:numCols %Ncode changed to 16 for Vert + for m=1:numRows %Ncode change to 24 for Vert + locIndx=locIndx+1; + + for k=1:numOfImages + if SWprint==1 + if k==1 && numOfImages==1, fprintf(fid,'%.2f\n',plate(m,n,k,pl));end + if k==1 && numOfImages~=1, fprintf(fid,'%.2f',plate(m,n,k,pl));end + if k>1 && k>>>>>>>>>>Call CurveFit routine >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + %NCfitImCF(times, CFscanIntens, offsets, suffix, AUCfinalTime, arrayFormat, scanAreas, currDirResults, autoImCFlag, selScan,sols) %,scan)%, scPltList) outputDirectory; + %[scanTmp2]= NCfitImCFparforFailGbl(times, CFscanIntens, offsets, suffix, AUCfinalTime, arrayFormat, scanAreas, currDirResults, autoImCFlag, selScan, sols, scanTmp); %,scan)%, scPltList) outputDirectory; + + [par4scanselIntensStd,par4scanselTimesStd,par4scanTimesELr,par4scanIntensELr,par4scanCFparameters,par4scanCFdate,outC2,outCstd2]= ... + NCfitImCFparforFailGbl2(parMat,times, CFscanIntens, offsets, suffix, AUCfinalTime, arrayFormat, scanAreas, currDirResults, autoImCFlag, selScan, sols); %,scan)%, scPltList) outputDirectory; +%<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< +else %fill with default values when an invalid plate scan occurs + CFscanIntens= zeros(16,24); + Ag= zeros(1,384); + %Preallocation for parfor loop********************************** + +times= tSeriesv; +st(1,1:size(times,2))= 1111; +resMat(1,1:27)= 0; +resMatStd= resMat; + outC2= zeros(384,27); + outCstd2= zeros(384,27); +for m=1:384 + pa{m}= st; + + par4scanCFparameters{m}= parMat; + par4scanCFdate{m}= datestr((now),31); +end + par4scanselTimesStd= pa; + par4scanselIntensStd= pa; + par4scanTimesELr= pa; + par4scanIntensELr= pa; + par4resMat= zeros(384,27); + par4resMatStd= zeros(384,27); +%{ + TmpexpScanIntens00= cell(1); %cell(1,scanMax); +TmpFexpScanSpots00= cell(1); %cell(1,scanMax); +TmpFexpScanBMtp00= cell(1); %cell(1,scanMax); +TmpanlZoneRefs00= cell(1); %cell(1,scanMax); +Tmpsbdg00= cell(1); +%} + +TmpexpScanIntens3= cell(1); %TmpexpScanIntens2; +TmpFexpScanSpots3= cell(1); %TmpFexpScanSpots2; +TmpFexpScanBMtp3= cell(1); %TmpFexpScanBMtp2; +TmpanlZoneRefs3= cell(1); %TmpanlZoneRefs2; +Tmpsbdg3= cell(1); %Tmpsbdg2; +%**************************************************************** + + +end %if tptLenth>=3 line19 20_0928 + + + + +tSeriesv; %debuggin parfor +p4L1{1}= tSeriesv; +p4L1{2}= t0Seriesv; +p4L1{3}= datestr(prtTmNumv,31); +p4L1{4}= CFscanIntens; + + + locIndx=0; + for n=1:numCols %Ncode changed to 16 for Vert + for m=1:numRows %Ncode change to 24 for Vert + locIndx=locIndx+1; + rc=[n,m]; + p4rcTmp(locIndx)={rc}; + p4pIndxTmp(locIndx)=locIndx; + end + end + + p4L1{5}= p4rcTmp; + p4L1{6}= p4pIndxTmp; + p4L1{7}= Ag; + + p4L1{8}= par4scanselIntensStd; + p4L1{9}= par4scanselTimesStd; + p4L1{10}= par4scanTimesELr; + p4L1{11}= par4scanIntensELr; + + p4L1{12}= par4scanCFparameters; + p4L1{13}= par4scanCFdate; + + p4L1{14}= outC2; + p4L1{15}= outCstd2; + p4L1{16}= selScan; + +p4L1{17}= cirPixA; +p4L1{18}= datestr((now),31); + +p4L2= p4L1; +end + diff --git a/workflow/templates/easy/par4Gbl_Main8c.m b/workflow/templates/easy/par4Gbl_Main8c.m new file mode 100755 index 00000000..155c9040 --- /dev/null +++ b/workflow/templates/easy/par4Gbl_Main8c.m @@ -0,0 +1,342 @@ +%% CALLED BY EASYconsole.m %% + +%function NImStartupImCF02par4Gbl +Fflg= 1; %0; + +fclose('all'); %close all open files +clear ('plate2', 'scanIntens','Scanfiles','pathname','tifFileLstP4') + +%global SWsingleSc +global SWgrowthArea +global scLst +global ImParMat +%global CSearchRange +%global CSrchRng +global expDir +global SWnewExp +%Console globals******* +global openExpfile +global openExppath +global newExpfile +global newExppath +global ExpOutmat +global fhconsole +global ExpPath +global resDir +global wCodeDir +global matDir % + +global ImWidth +global ImHeigth +global numRows; +global numCols; +%global SWprintLeg +global scan +%global scanSize +global scanMax +global tptLength + +%global sols +%global CFmeth +numRows=24; %for Single Vertical +numCols=16; %for Single Vertical + +%global scanPar4x +%*************CIRCLE related************************* + +doCircle=1; %use Circle area analysis 14_0807 +radius=14; +ImParMat(10)=radius; +ImParMat(11)=doCircle; + +try +clf(fhconsole,'reset') +catch ME +end +close +%****************************************** +EASYconsole +%****************************************** +try +load(ExpOutmat) +copyfile(ExpOutmat,(fullfile(matDir,'BkUp',strcat((num2str(datenum(now))),'.mat')))) +catch +end + +try + load(fullfile(resDir,'Fotos','Nbdg')) %Modified to load from 'Fotos' 20_0819 +catch + %Reloacated from 'PTmats' to prevent potential overwrite when PTmats is + %copied into new job when the PT template is about the same. We also + %now have a default template if one is not made. i.e., when the images + %from the new experiment are too sketchy to make a good pintool + %template. By moving it to 'Fotos' we avoid possible issues due to + %copying the Nbdg.mat file along with the default template '.mat' files. + %A copy of Ndbg.mat is placed also saved to the 'PTmats' directory + %after each run to allow previous version of EASY to access data made + %by EASY versions after 20_0819. + + load(fullfile(resDir,'PTmats','Nbdg')) %Left in to accomodate loads of work before 20_0819 +end +%++++++++++++++Load Fotos stored data++++++++++++++++++++ +try + load(fullfile(resDir,'Fotos','Coordinates')) +catch +end +try + load(fullfile(resDir,'Fotos','BGatTpts')) +catch +end +try + load(fullfile(resDir,'Fotos','anlZones')) +catch +end +%******autoImCF******************* +try +load(fullfile(resDir,'PTmats','NCFparms')) +catch ME + load parameters +end +%******************************* +%************Get Print Times******************************* +PrintTimes=[]; + w= cd; %Capture current dir in w + if ispc %Linux accommodation + expDir=fullfile(ExpPath,'\'); + else + if isunix, expDir= ExpPath;end + if ismac , expDir=ExpPath;end + end +scLst={}; +%**************Parameter Entry****************** +NImParamRadiusGui(expDir) %Ncode 122111replaced removed ,numOfPrtTimes) +%*********************************************** +%width=ImParMat(5); +width=24; +widthEx=width-1; %width extention from reference point +dither=ImParMat(6); +radius=ImParMat(10); +%>>>>>>>>>>>>>>>>>>>>>> +NIcircle +%<<<<<<<<<<<<<<<<<<<<<<< + +%***************** Load Stuff************************* + +lastPlateOnLastScan=1; %Ncode + +if size(scLst,1)==1,SWsingleSc=1;else SWsingleSc=0;end +%**************************???????????????????????? + +dvec=datevec(datestr(floor(now))); %method to get current offset year '01-Jan-"currentyr"' +dvec(2)=1; +dvec(3)=1; +%yrOffset=datenum('01-Jan-2012');%(dvec); %('01-Jan-2009'); %time num default is currentyear jan1 +%******************+++++++++++++*****************+++++++++++++++*********** +numScans=size(scLst,1); + +if(isequal(SWsingleSc,1)) + selScan=str2double(char(scLst(1))); + %numScans=selScan; + %startScan=selScan; +else + %startScan=1; +end + +SWgrowthArea= ImParMat(9); +load(fullfile(resDir,'PTmats','NPTmapSearch')) + PTmapPos= detPos; + +%************************************************************************** +selScanNumLst= []; + +%++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +%************************************************************************** +%*********Preallocation for Scan Loop************************************************** + +Scanfiles=[]; +pathname=[]; +%++++++++++++++++++++++++++++++++ SETUP FOR GLOBAL PARFOR LOOP+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +%if SWsingleSc == 1 %if single scan else do expMultiTseries to execute multiple scans +%******Gui Select Scan -> *.tif files for Analysis***************** + w= cd; %Capture current dir in w +if ispc +cd(char(strcat(expDir,scLst(1)))) +end +if isunix +cd (fullfile(expDir,cell2mat(scLst(1)))); +end +Scanfiles=[]; +pathname=[]; + + +for ii= 1:length(scLst) + if SWsingleSc == 1 +[Scanfiles, pathname]=uigetfile('*.bmp', 'Select files','MultiSelect','on');% change '*hr*.bmp' 12/20/2011 + if ischar(Scanfiles) + scd= imread(char(Scanfiles)); + tptLength=1; + + else + scd= imread(char(Scanfiles(1,1))); + tptLength= length(Scanfiles); + end + ImHeigth= size(scd,1); + ImWidth= size(scd,2); + sc=scd(1:ImHeigth,1:ImWidth); %?Not used RefOnly + end + +cd (w); +numFiles= size(Scanfiles,2); +%tifFileLst= Scanfiles; + +%Initialize tifFilesLst for parfor loop + if ispc %Linux and mac accomodation + tifFileLst4MultiT = dir(fullfile(expDir, char(scLst(ii)), '\*.bmp')); + end + if isunix || ismac + tifFileLst4MultiT = dir(fullfile(expDir, char(scLst(ii)), '*.bmp')); + end + numFiles=length(tifFileLst4MultiT); + tptLength= numFiles; + tifFileLstP4{ii}={tifFileLst4MultiT.name}; + %{ + dateTmStampP4{ii}= {tifFileLst4MultiT.date} + dateTmNum{ii}= datenum(dateTmStampP4{1}); + info{ii} = imfinfo(tifFileLstP4{ii}); + scTmNumericP4(ii)=datenum(info{ii}.FileModDate); + %} +end + + + + %************************************************************************** + %************************************************************************** + +for jj= 1:numScans %startScan:numScans + selScan= str2double(char(scLst(jj))); + selScanNumLst(jj)= selScan; +end + selScanNumLst2= selScanNumLst; %function passthrough, passback to par4gbl_Main 20_0205 +%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +%PREALLOCATION PREALLOCATION PREALLOCATION PREALLOCATION PREALLOCATION +%****************************************************************************** +scCount= []; + +TmpexpScanIntens00= cell(1); %cell(1,scanMax); +TmpFexpScanSpots00= cell(1); %cell(1,scanMax); +TmpFexpScanBMtp00= cell(1); %cell(1,scanMax); +TmpanlZoneRefs00= cell(1); %cell(1,scanMax); +Tmpsbdg00= cell(1); + +TmpexpScanIntens4= cell(1,numScans); +TmpFexpScanSpots4= cell(1,numScans); +TmpFexpScanBMtp4= cell(1,numScans); +TmpanlZoneRefs4= cell(1,numScans); +Tmpsbdg4= cell(1,numScans); +TmpexpScanIntens5= cell(1,numScans); +TmpFexpScanSpots5= cell(1,numScans); +TmpFexpScanBMtp5= cell(1,numScans); +TmpanlZoneRefs5= cell(1,numScans); +Tmpsbdg5= cell(1,numScans); + +p4L00= cell(18,1); +p4L0= p4L00; + +p4L4= cell(18,numScans); +p4L5= p4L4; + + Ag= ones(384); + CFscanIntens= zeros(16,24); + + + [p4L4,... + TmpexpScanIntens5,TmpFexpScanSpots5,TmpFexpScanBMtp5,TmpanlZoneRefs5,Tmpsbdg5]= ... + p4loop8c(parMat,tptLength,numScans,selScanNumLst,SWsingleSc,Fflg,PTmapPos,optCirMask,diaExt,doCircle,cirPixA,cirMask,width, ... + TmpexpScanIntens00,TmpFexpScanSpots00,TmpFexpScanBMtp00,TmpanlZoneRefs00,scCount,tifFileLstP4,pathname,ImParMat, ... + numRows,numCols,scLst,resDir,expDir, p4L00,TmpexpScanIntens4,TmpFexpScanSpots4,TmpFexpScanBMtp4,TmpanlZoneRefs4, ... + Tmpsbdg00,Tmpsbdg4); + + + +for scanCnt= 1:numScans +selScan= p4L4{16,scanCnt}; %determine the actual scan in the scanCnt parfor distributed "id" + scan(selScan).plate(1).tSeries= cell2mat(p4L4(1,scanCnt)); + scan(selScan).plate(1).t0Series= cell2mat(p4L4(2,scanCnt)); + scan(selScan).plate(1).printTm= cell2mat(p4L4(3,scanCnt)); + scan(selScan).plate(1).intens= cell2mat(p4L4(4,scanCnt)); + scan(selScan).plate(1).rc= p4L4(5,scanCnt); + scan(selScan).plate(1).pIndx= cell2mat(p4L4(6,scanCnt)); + scan(selScan).plate(1).Ag= cell2mat(p4L4(7,scanCnt)); + + scan(selScan).plate(1).selIntens= p4L4(8,scanCnt); + scan(selScan).plate(1).selTimes= p4L4(9,scanCnt); + scan(selScan).plate(1).filterTimes= p4L4(10,scanCnt); + scan(selScan).plate(1).normIntens= p4L4(11,scanCnt); + + %scan(selScan).plate(1).CFparameters= p4L4(12,scanCnt); %Need to convert to a matrix form like Old versions + CFparm(1:384)=p4L4{12,scanCnt}(1:384); + scan(selScan).plate(1).CFparameters= CFparm; + scan(selScan).plate(1).CFdate= p4L4(13,scanCnt); + scan(selScan).plate(1).CFout= cell2mat(p4L4(14,scanCnt)); + scan(selScan).plate(1).CFoutStd= cell2mat(p4L4(15,scanCnt)); + + + scan(selScan).Awindow= cell2mat(p4L4(17,scanCnt)); + scan(selScan).Idate= cell2mat(p4L4(18,scanCnt)); + + + expScanIntens(selScan)= TmpexpScanIntens5(scanCnt); + FexpScanSpots(selScan)= TmpFexpScanSpots5(scanCnt); + FexpScanBMtp(selScan)= TmpFexpScanBMtp5(scanCnt); + anlZoneRefs(selScan)= TmpanlZoneRefs5(scanCnt); + if ~isempty(Tmpsbdg5{scanCnt}) + sbdg(selScan)= Tmpsbdg5(scanCnt); + else + sbdg{selScan}= uint8(zeros(24,16,4)); + end +end + +%SAVE DATA in .mat files +save(ExpOutmat,'scan'); +save((fullfile(resDir,'PTmats','Nbdg')), 'sbdg'); %legacy location can probably get rid of in time +save((fullfile(resDir,'Fotos','Nbdg')), 'sbdg'); + +save((fullfile(resDir,'Fotos','Coordinates')),'FexpScanSpots') %Saves frames at each tPt +save((fullfile(resDir,'Fotos','BGatTpts')),'FexpScanBMtp') +save((fullfile(resDir,'Fotos','anlZones')),'anlZoneRefs')%Saves anl Positions at each tPt + +%Print FitResults ******************************************************* +fileExt = '.txt'; +filePrefix = 'FitResults_'; +for scanCnt= 1:numScans +selScan= p4L4{16,scanCnt}; %determine the actual scan in the scanCnt parfor distributed "id" + +fileSuffix= strcat('Scan', num2str(selScan),'_Plate', num2str(1)); +fileNamePlate = [filePrefix fileSuffix fileExt]; +fileName= fullfile(resDir,'PrintResults', fileNamePlate); %[outputDirectory fileNamePlate]; + +%This,fprint for loop,is an very old legacy feature which slows processing. Could be +%removed but allows easy observation of how a run is progressing and can be +%used as a diagnostic tool. + + +outCprint= p4L4; + +fid = fopen(fileName,'w'); +fprintf(fid, 'Num.\tAUC\tMSR\tK\tr\tl\tR-squared\tK-lower\tK-upper\tr-lower\tr-upper\tl-upper\tl-lower\tArea\tLastInten\tSpineMaxRateTimePt\tLastFitTimePt\n'); +for n=1:384 %startCount:numCultures + fprintf(fid,'%d\t',n); + fprintf(fid, '%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\n',... + outCprint{14,scanCnt}(n,1),outCprint{14,scanCnt}(n,2),outCprint{14,scanCnt}(n,3),outCprint{14,scanCnt}(n,4),... + outCprint{14,scanCnt}(n,5),outCprint{14,scanCnt}(n,6),outCprint{14,scanCnt}(n,7),outCprint{14,scanCnt}(n,8),... + outCprint{14,scanCnt}(n,9),outCprint{14,scanCnt}(n,10),outCprint{14,scanCnt}(n,11),outCprint{14,scanCnt}(n,12),... + outCprint{14,scanCnt}(n,13),outCprint{14,scanCnt}(n,14),outCprint{14,scanCnt}(n,15),outCprint{14,scanCnt}(n,16)); +end +fclose(fid); +end + +EASYconsole; + + + %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ \ No newline at end of file diff --git a/workflow/templates/qhtcp/ExpTemplate/Z_InteractionTemplate.R b/workflow/templates/qhtcp/ExpTemplate/Z_InteractionTemplate.R index fe31ca14..97a737cc 100644 --- a/workflow/templates/qhtcp/ExpTemplate/Z_InteractionTemplate.R +++ b/workflow/templates/qhtcp/ExpTemplate/Z_InteractionTemplate.R @@ -1,6 +1,15 @@ -#Based on InteractionTemplate.R which is based on Sean Santose's Interaction_V5 script. -#Adapt SS For Structured Data storage but using command line scripts -###Set up the required libraries, call required plot theme elements and set up the command line arguments +#!/usr/bin/env R +# Based on InteractionTemplate.R which is based on Sean Santos's Interaction_V5 script +# +# Updated 2024 Bryan C Roessler to improve file operations and portability +# NOTE: The script now has 1 REQUIRED argument and 4 OPTIONAL arguments: +# 1. Path to input file +# 2. /output/ directory +# 3. Path to StudyInfo.csv +# 4. Path to SGDgeneList +# 5. Standard deviation value + +# Load libraries library("ggplot2") library("plyr") library("extrafont") @@ -8,73 +17,76 @@ library("gridExtra") library("gplots") library("RColorBrewer") library("stringr") -#library("gdata") -library(plotly) -library(htmlwidgets) +library("gdata") +library("plotly") +library("htmlwidgets") -Args <- commandArgs(TRUE) -input_file <- Args[1] #"!!Results_17_0827_yor1null-rpl12anull misLabeledAsFrom MI 17_0919_yor1-curated.txt" #Args[1] #Arg 1 #"!!ResultsStd_JS 19_1224_HLEG_P53.txt" is the !!results ... .txt -#Tool to find a file and copy it to desired location -destDir= getwd() -#srcFile= file.choose() -#file.copy(srcFile, destDir) -#input_file= tail(strsplit(srcFile,"[/]")[[1]],1) +# Parse arguments +args <- commandArgs(TRUE) +# Set input file (required) +inputFile <- args[1] - -#Path to Output Directory -#W=getwd() #R is F'd up, Can't use, Any legitamate platform could build out dirs from this -outDir <- "ZScores/" #"Args[2] #paste0(W,"/ZScores/") -subDir <- outDir #Args[2] - -if (file.exists(subDir)){ - outputpath <- subDir +# Set output dir +if (length(args) > 2) { + outDir <- args[2] } else { - dir.create(file.path(subDir)) + outDir <- "/ZScores/" } -if (file.exists(paste(subDir,"QC/",sep=""))){ - outputpath_QC <- paste(subDir,"QC/",sep="") +# Set StudyInfo file path +if (length(args) > 3) { + studyInfo <- args[3] } else { - dir.create(file.path(paste(subDir,"QC/",sep=""))) - outputpath_QC <- paste(subDir,"QC/",sep="") + studyInfo <- "../Code/StudyInfo.csv" } -#define the output path (formerly the second argument from Rscript) -outputpath <- outDir -#Set Args[2] the Background contamination noise filter as a function of standard deviation -#std= as.numeric(Args[2]) -#Sean recommends 3 or 5 SD factor. -#Capture Exp_ number,use it to Save Args[2]{std}to Labels field and then Write to Labels to studyInfo.txt for future reference -Labels <- read.csv(file= "../Code/StudyInfo.csv",stringsAsFactors = FALSE) #,sep= ",") -print("Be sure to enter Background noise filter standard deviation i.e., 3 or 5 per Sean") +# Set SGDgeneList file path +if (length(args) > 4) { + SGDgeneList <- args[4] +} else { + SGDgeneList <- "../Code/SGD_features.tab" +} -#User prompt for std multiplier Value -cat("Enter a Standard Deviation value to noise filter \n") -inpChar<- readLines(file("stdin"), n = 1L) -cat(paste("Standard Deviation Value is", inpChar, "\n")) -inpNum= as.numeric(inpChar) -#set std deviation multiplier default if no user entry -if(!is.na(inpNum)){ - std= inpNum -}else{std= 3} +# Set standard deviation +if (length(args) > 5) { + delBGFactor <- args[5] +} else { + # User prompt for std multiplier Value + print("Enter a Standard Deviation value for noise filter") + print("Sean Santos recommends 3 or 5") + delBGFactor <- readLines(file("stdin"), n = 1L) +} +delBGFactor <- as.numeric(delBGFactor) +if(is.na(delBGFactor)){ + delBGFactor <- 3 # Recommended by Sean +} +print(paste("The Standard Deviation Value is: ", delBGFactor)) + +outDir_QC <- paste(outDir,"QC/",sep="") + +if !(file.exists(outDir)){ + dir.create(file.path(outDir)) +} + +if !(file.exists(outDir_QC)){ + dir.create(file.path(outDir_QC)) +} + +# Capture Exp_ number,use it to Save args[2]{std}to Labels field and then +# write to Labels to studyInfo.txt for future reference +Labels <- read.csv(file=studyInfo,stringsAsFactors = FALSE) # sep= "," +expNumber <- as.numeric(sub("^.*?(\\d+)$", "\\1", getwd())) +Labels[expNumber,3] <- delBGFactor +write.csv(Labels,file=studyInfo,row.names = FALSE) -expNumber<- as.numeric(sub("^.*?(\\d+)$", "\\1", getwd())) -Labels[expNumber,3]= as.numeric(std) -Delta_Background_sdFactor <- std -DelBGFactr <- as.numeric(Delta_Background_sdFactor) -#Write Background SD value to studyInfo.txt file -#write.csv(Labels,file=paste("../Code/StudyInfo.csv"),row.names = FALSE) -write.csv(Labels,file=paste("../Code/StudyInfo.csv"),row.names = FALSE) -print('ln 50 write StudyInfo.csv ') -#write.table(Labels,file=paste(outputpath,"StudyInfo.txt"),sep = "\t",row.names = FALSE) - -#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -#++++++BEGIN USER DATA SELECTION SECTION+++++++++++++++++++++++++++++++++++++++++++++++++ +############################################################################### +################### BEGIN USER DATA SELECTION SECTION ######################### +############################################################################### #read in the data -X <- read.delim(input_file,skip=2,as.is=T,row.names=1,strip.white=TRUE) +X <- read.delim(inputFile,skip=2,as.is=T,row.names=1,strip.white=TRUE) X <- X[!(X[[1]]%in%c("","Scan")),] #X <- X[!(X[[1]]%in%c(61:76)),] #Remove dAmp plates which are Scans 61 thru 76 @@ -145,7 +157,6 @@ Xbu= X #Inserted to use SGDgenelist to update orfs and replace empty geneName cells with ORF name (adapted from Sean's Merge script). This is to 'fix' the naming for everything that follows (REMc, Heatmaps ... et.al) rather than do it piece meal later #Sean's Match Script( which was adapted here) was fixed 2022_0608 so as not to write over the RF1&RF2 geneNames which caused a variance with his code results #in the Z_lm_L,K,r&AUC output values. Values correlated well but were off by a multiplier factor. -SGDgeneList= "../Code/SGD_features.tab" genes = data.frame(read.delim(file=SGDgeneList,quote="",header=FALSE,colClasses = c(rep("NULL",3), rep("character", 2), rep("NULL", 11)))) for(i in 1:length(X[,14])){ ii= as.numeric(i) @@ -315,17 +326,17 @@ X$Delta_Backgrd <- X$LstBackgrd - X$X1stBackgrd Raw_l_vs_K_beforeQC <- ggplot(X,aes(l,K,color=as.factor(Conc_Num))) + geom_point(aes(ORF=ORF,Gene=Gene,Delta_Backgrd=Delta_Backgrd),shape=3) + ggtitle("Raw L vs K before QC") + theme_Publication_legend_right() -pdf(paste(outputpath_QC,"Raw_L_vs_K_beforeQC.pdf",sep=""),width = 12,height = 8) +pdf(paste(outDir_QC,"Raw_L_vs_K_beforeQC.pdf",sep=""),width = 12,height = 8) Raw_l_vs_K_beforeQC dev.off() pgg <- ggplotly(Raw_l_vs_K_beforeQC) #pgg -plotly_path <- paste(getwd(),"/",outputpath_QC,"Raw_L_vs_K_beforeQC.html",sep="") +plotly_path <- paste(outDir_QC,"Raw_L_vs_K_beforeQC.html",sep="") saveWidget(pgg, file=plotly_path, selfcontained =TRUE) #set delta background tolerance based on 3 sds from the mean delta background -Delta_Background_Tolerance <- mean(X$Delta_Backgrd)+(DelBGFactr*sd(X$Delta_Backgrd)) +Delta_Background_Tolerance <- mean(X$Delta_Backgrd)+(delBGFactor*sd(X$Delta_Backgrd)) #Delta_Background_Tolerance <- mean(X$Delta_Backgrd)+(3*sd(X$Delta_Backgrd)) print(paste("Delta_Background_Tolerance is",Delta_Background_Tolerance,sep=" ")) @@ -350,12 +361,12 @@ X_Delta_Backgrd_above_Tolerance_L_vs_K <- ggplot(X_Delta_Backgrd_above_Tolerance annotate("text",x=X_Delta_Backgrd_above_Tolerance_L_halfmedian,y=X_Delta_Backgrd_above_Tolerance_K_halfmedian, label = paste("Strains above delta background tolerance = ",X_Delta_Backgrd_above_Tolerance_toRemove)) + theme_Publication_legend_right() -pdf(paste(outputpath_QC,"Raw_L_vs_K_for_strains_above_deltabackgrd_threshold.pdf",sep=""),width = 12,height = 8) +pdf(paste(outDir_QC,"Raw_L_vs_K_for_strains_above_deltabackgrd_threshold.pdf",sep=""),width = 12,height = 8) X_Delta_Backgrd_above_Tolerance_L_vs_K dev.off() pgg <- ggplotly(X_Delta_Backgrd_above_Tolerance_L_vs_K) #pgg -plotly_path <- paste(getwd(),"/",outputpath_QC,"Raw_L_vs_K_for_strains_above_deltabackgrd_threshold.html",sep="") +plotly_path <- paste(outDir_QC,"Raw_L_vs_K_for_strains_above_deltabackgrd_threshold.html",sep="") saveWidget(pgg, file=plotly_path, selfcontained =TRUE) #frequency plot for all data vs. the delta_background @@ -366,7 +377,7 @@ DeltaBackground_Frequency_Plot <- ggplot(X,aes(Delta_Backgrd,color=as.factor(Con DeltaBackground_Bar_Plot <- ggplot(X,aes(Delta_Backgrd,color=as.factor(Conc_Num))) + geom_bar() + ggtitle("Bar plot for Delta Background by Conc All Data") + theme_Publication_legend_right() -pdf(file = paste(outputpath_QC,"Frequency_Delta_Background.pdf",sep=""),width = 12,height = 8) +pdf(file = paste(outDir_QC,"Frequency_Delta_Background.pdf",sep=""),width = 12,height = 8) print(DeltaBackground_Frequency_Plot) print(DeltaBackground_Bar_Plot) dev.off() @@ -432,7 +443,7 @@ Plate_Analysis_Delta_Backgrd_Box_afterQC <- ggplot(X,aes(as.factor(Scan),Delta_B ggtitle("Plate analysis by Drug Conc for Delta_Backgrd after quality control") + theme_Publication() #print the plate analysis data before and after QC -pdf(file=paste(outputpath_QC,"Plate_Analysis.pdf",sep=""),width = 14,height=9) +pdf(file=paste(outDir_QC,"Plate_Analysis.pdf",sep=""),width = 14,height=9) Plate_Analysis_L Plate_Analysis_L_afterQC Plate_Analysis_K @@ -446,7 +457,7 @@ Plate_Analysis_Delta_Backgrd_afterQC dev.off() #print the plate analysis data before and after QC -pdf(file=paste(outputpath_QC,"Plate_Analysis_Boxplots.pdf",sep=""),width = 18,height=9) +pdf(file=paste(outDir_QC,"Plate_Analysis_Boxplots.pdf",sep=""),width = 18,height=9) Plate_Analysis_L_Box Plate_Analysis_L_Box_afterQC Plate_Analysis_K_Box @@ -503,7 +514,7 @@ Plate_Analysis_Delta_Backgrd_Box_afterQC_Z <- ggplot(X_noZero,aes(as.factor(Scan ggtitle("Plate analysis by Drug Conc for Delta_Backgrd after quality control") + theme_Publication() #print the plate analysis data before and after QC -pdf(file=paste(outputpath_QC,"Plate_Analysis_noZeros.pdf",sep=""),width = 14,height=9) +pdf(file=paste(outDir_QC,"Plate_Analysis_noZeros.pdf",sep=""),width = 14,height=9) Plate_Analysis_L_afterQC_Z Plate_Analysis_K_afterQC_Z Plate_Analysis_r_afterQC_Z @@ -512,7 +523,7 @@ Plate_Analysis_Delta_Backgrd_afterQC_Z dev.off() #print the plate analysis data before and after QC -pdf(file=paste(outputpath_QC,"Plate_Analysis_noZeros_Boxplots.pdf",sep=""),width = 18,height=9) +pdf(file=paste(outDir_QC,"Plate_Analysis_noZeros_Boxplots.pdf",sep=""),width = 18,height=9) Plate_Analysis_L_Box_afterQC_Z Plate_Analysis_K_Box_afterQC_Z Plate_Analysis_r_Box_afterQC_Z @@ -556,7 +567,7 @@ X_stats_ALL <- ddply(X, c("Conc_Num","Conc_Num_Factor"), summarise, ) #print(X_stats_ALL_L) -write.csv(X_stats_ALL,file=paste(outputpath,"SummaryStats_ALLSTRAINS.csv"),row.names = FALSE) +write.csv(X_stats_ALL,file=paste(outDir,"SummaryStats_ALLSTRAINS.csv"),row.names = FALSE) #+++++END QC SECTION+++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##### Part 3 - Generate summary statistics and calculate the max theoretical L value @@ -660,7 +671,7 @@ for(s in Background_Strains){ se_AUC = sd_AUC / sqrt(N-1) ) - write.csv(X_stats_BY,file=paste(outputpath,"SummaryStats_BackgroundStrains.csv"),row.names=FALSE) + write.csv(X_stats_BY,file=paste(outDir,"SummaryStats_BackgroundStrains.csv"),row.names=FALSE) #calculate the max theoretical L values #only look for max values when K is within 2SD of the ref strain @@ -708,7 +719,7 @@ for(s in Background_Strains){ ) print(X_stats_BY_L_within_2SD_K) X1_SD_within_2SD_K <- max(X_stats_BY_L_within_2SD_K$sd) - write.csv(X_stats_BY_L_within_2SD_K,file=paste(outputpath_QC,"Max_Observed_L_Vals_for_spots_within_2SD_K.csv",sep=""),row.names=FALSE) + write.csv(X_stats_BY_L_within_2SD_K,file=paste(outDir_QC,"Max_Observed_L_Vals_for_spots_within_2SD_K.csv",sep=""),row.names=FALSE) X_stats_BY_L_outside_2SD_K <- ddply(X_outside_2SD_K, c("Conc_Num","Conc_Num_Factor"), summarise, N = (length(l)), @@ -725,22 +736,22 @@ for(s in Background_Strains){ #X1_SD_outside_2SD_K <- X[X$l %in% X1_SD_within_2SD_K$l,] Outside_2SD_K_L_vs_K <- ggplot(X_outside_2SD_K,aes(l,K,color=as.factor(Conc_Num))) + geom_point(aes(ORF=ORF,Gene=Gene,Delta_Backgrd=Delta_Backgrd),shape=3) + ggtitle("Raw L vs K for strains falling outside 2SD of the K mean at each conc") + theme_Publication_legend_right() - pdf(paste(outputpath_QC,"Raw_L_vs_K_for_strains_2SD_outside_mean_K.pdf",sep=""),width = 10,height = 8) + pdf(paste(outDir_QC,"Raw_L_vs_K_for_strains_2SD_outside_mean_K.pdf",sep=""),width = 10,height = 8) print(Outside_2SD_K_L_vs_K) dev.off() pgg <- ggplotly(Outside_2SD_K_L_vs_K) - plotly_path <- paste(getwd(),"/",outputpath_QC,"RawL_vs_K_for_strains_outside_2SD_K.html",sep="") + plotly_path <- paste(outDir_QC,"RawL_vs_K_for_strains_outside_2SD_K.html",sep="") saveWidget(pgg, file=plotly_path, selfcontained =TRUE) Outside_2SD_K_delta_background_vs_K <- ggplot(X_outside_2SD_K,aes(Delta_Backgrd,K,color=as.factor(Conc_Num))) + geom_point(aes(l=l,ORF=ORF,Gene=Gene),shape=3,position="jitter") + ggtitle("DeltaBackground vs K for strains falling outside 2SD of the K mean at each conc") + theme_Publication_legend_right() - pdf(paste(outputpath_QC,"DeltaBackground_vs_K_for_strains_2SD_outside_mean_K.pdf",sep=""),width = 10,height = 8) + pdf(paste(outDir_QC,"DeltaBackground_vs_K_for_strains_2SD_outside_mean_K.pdf",sep=""),width = 10,height = 8) Outside_2SD_K_delta_background_vs_K dev.off() pgg <- ggplotly(Outside_2SD_K_delta_background_vs_K) #pgg - plotly_path <- paste(getwd(),"/",outputpath_QC,"DeltaBackground_vs_K_for_strains_outside_2SD_K.html",sep="") + plotly_path <- paste(outDir_QC,"DeltaBackground_vs_K_for_strains_outside_2SD_K.html",sep="") saveWidget(pgg, file=plotly_path, selfcontained =TRUE) @@ -1160,7 +1171,7 @@ for(s in Background_Strains){ InteractionScores_RF <- InteractionScores_RF[order(InteractionScores_RF$Z_lm_L,decreasing=TRUE),] InteractionScores_RF <- InteractionScores_RF[order(InteractionScores_RF$NG,decreasing=TRUE),] - write.csv(InteractionScores_RF,paste(outputpath,"RF_ZScores_Interaction.csv",sep=""),row.names=FALSE) + write.csv(InteractionScores_RF,paste(outDir,"RF_ZScores_Interaction.csv",sep=""),row.names=FALSE) for(i in 1:num_genes_RF){ @@ -1234,7 +1245,7 @@ for(s in Background_Strains){ } } print("Pass RF ggplot loop") - write.csv(X_stats_interaction_ALL_RF_final,paste(outputpath,"RF_ZScore_Calculations.csv",sep=""),row.names = FALSE) + write.csv(X_stats_interaction_ALL_RF_final,paste(outDir,"RF_ZScore_Calculations.csv",sep=""),row.names = FALSE) ####### Part 5 - Get Zscores for Gene deletion strains @@ -1564,7 +1575,7 @@ for(s in Background_Strains){ InteractionScores <- InteractionScores[order(InteractionScores$NG,decreasing=TRUE),] df_order_by_OrfRep <- unique(InteractionScores$OrfRep) #X_stats_interaction_ALL <- X_stats_interaction_ALL[order(X_stats_interaction_ALL$NG,decreasing=TRUE),] - write.csv(InteractionScores,paste(outputpath,"ZScores_Interaction.csv",sep=""),row.names=FALSE) + write.csv(InteractionScores,paste(outDir,"ZScores_Interaction.csv",sep=""),row.names=FALSE) InteractionScores_deletion_enhancers_L <- InteractionScores[InteractionScores$Avg_Zscore_L >= 2,] InteractionScores_deletion_enhancers_K <- InteractionScores[InteractionScores$Avg_Zscore_K <= -2,] @@ -1588,16 +1599,16 @@ for(s in Background_Strains){ InteractionScores_deletion_enhancers_lm_Suppressors_AvgZscore_K <- InteractionScores_deletion_enhancers_lm_Suppressors_AvgZscore_K[!is.na(InteractionScores_deletion_enhancers_lm_Suppressors_AvgZscore_K$OrfRep),] InteractionScores_deletion_enhancers_Avg_Zscore_Suppressors_lm_K <- InteractionScores_deletion_enhancers_Avg_Zscore_Suppressors_lm_K[!is.na(InteractionScores_deletion_enhancers_Avg_Zscore_Suppressors_lm_K$OrfRep),] - write.csv(InteractionScores_deletion_enhancers_L,paste(outputpath,"ZScores_Interaction_DeletionEnhancers_L.csv",sep=""),row.names=FALSE) - write.csv(InteractionScores_deletion_enhancers_K,paste(outputpath,"ZScores_Interaction_DeletionEnhancers_K.csv",sep=""),row.names=FALSE) - write.csv(InteractionScores_deletion_suppressors_L,paste(outputpath,"ZScores_Interaction_DeletionSuppressors_L.csv",sep=""),row.names=FALSE) - write.csv(InteractionScores_deletion_suppressors_K,paste(outputpath,"ZScores_Interaction_DeletionSuppressors_K.csv",sep=""),row.names=FALSE) - write.csv(InteractionScores_deletion_enhancers_and_Suppressors_L,paste(outputpath,"ZScores_Interaction_DeletionEnhancers_and_Suppressors_L.csv",sep=""),row.names=FALSE) - write.csv(InteractionScores_deletion_enhancers_and_Suppressors_K,paste(outputpath,"ZScores_Interaction_DeletionEnhancers_and_Suppressors_K.csv",sep=""),row.names=FALSE) - write.csv(InteractionScores_deletion_enhancers_lm_Suppressors_AvgZscore_L,paste(outputpath,"ZScores_Interaction_Suppressors_and_lm_Enhancers_L.csv",sep=""),row.names=FALSE) - write.csv(InteractionScores_deletion_enhancers_Avg_Zscore_Suppressors_lm_L,paste(outputpath,"ZScores_Interaction_Enhancers_and_lm_Suppressors_L.csv",sep=""),row.names=FALSE) - write.csv(InteractionScores_deletion_enhancers_lm_Suppressors_AvgZscore_K,paste(outputpath,"ZScores_Interaction_Suppressors_and_lm_Enhancers_K.csv",sep=""),row.names=FALSE) - write.csv(InteractionScores_deletion_enhancers_Avg_Zscore_Suppressors_lm_K,paste(outputpath,"ZScores_Interaction_Enhancers_and_lm_Suppressors_K.csv",sep=""),row.names=FALSE) + write.csv(InteractionScores_deletion_enhancers_L,paste(outDir,"ZScores_Interaction_DeletionEnhancers_L.csv",sep=""),row.names=FALSE) + write.csv(InteractionScores_deletion_enhancers_K,paste(outDir,"ZScores_Interaction_DeletionEnhancers_K.csv",sep=""),row.names=FALSE) + write.csv(InteractionScores_deletion_suppressors_L,paste(outDir,"ZScores_Interaction_DeletionSuppressors_L.csv",sep=""),row.names=FALSE) + write.csv(InteractionScores_deletion_suppressors_K,paste(outDir,"ZScores_Interaction_DeletionSuppressors_K.csv",sep=""),row.names=FALSE) + write.csv(InteractionScores_deletion_enhancers_and_Suppressors_L,paste(outDir,"ZScores_Interaction_DeletionEnhancers_and_Suppressors_L.csv",sep=""),row.names=FALSE) + write.csv(InteractionScores_deletion_enhancers_and_Suppressors_K,paste(outDir,"ZScores_Interaction_DeletionEnhancers_and_Suppressors_K.csv",sep=""),row.names=FALSE) + write.csv(InteractionScores_deletion_enhancers_lm_Suppressors_AvgZscore_L,paste(outDir,"ZScores_Interaction_Suppressors_and_lm_Enhancers_L.csv",sep=""),row.names=FALSE) + write.csv(InteractionScores_deletion_enhancers_Avg_Zscore_Suppressors_lm_L,paste(outDir,"ZScores_Interaction_Enhancers_and_lm_Suppressors_L.csv",sep=""),row.names=FALSE) + write.csv(InteractionScores_deletion_enhancers_lm_Suppressors_AvgZscore_K,paste(outDir,"ZScores_Interaction_Suppressors_and_lm_Enhancers_K.csv",sep=""),row.names=FALSE) + write.csv(InteractionScores_deletion_enhancers_Avg_Zscore_Suppressors_lm_K,paste(outDir,"ZScores_Interaction_Enhancers_and_lm_Suppressors_K.csv",sep=""),row.names=FALSE) #get enhancers and suppressors for linear regression InteractionScores_deletion_enhancers_L_lm <- InteractionScores[InteractionScores$Z_lm_L >= 2,] @@ -1614,15 +1625,15 @@ for(s in Background_Strains){ InteractionScores_deletion_enhancers_and_Suppressors_L_lm <- InteractionScores_deletion_enhancers_and_Suppressors_L_lm[!is.na(InteractionScores_deletion_enhancers_and_Suppressors_L_lm$OrfRep),] InteractionScores_deletion_enhancers_and_Suppressors_K_lm <- InteractionScores_deletion_enhancers_and_Suppressors_K_lm[!is.na(InteractionScores_deletion_enhancers_and_Suppressors_K_lm$OrfRep),] - write.csv(InteractionScores_deletion_enhancers_L_lm,paste(outputpath,"ZScores_Interaction_DeletionEnhancers_L_lm.csv",sep=""),row.names=FALSE) - write.csv(InteractionScores_deletion_enhancers_K_lm,paste(outputpath,"ZScores_Interaction_DeletionEnhancers_K_lm.csv",sep=""),row.names=FALSE) - write.csv(InteractionScores_deletion_suppressors_L_lm,paste(outputpath,"ZScores_Interaction_DeletionSuppressors_L_lm.csv",sep=""),row.names=FALSE) - write.csv(InteractionScores_deletion_suppressors_K_lm,paste(outputpath,"ZScores_Interaction_DeletionSuppressors_K_lm.csv",sep=""),row.names=FALSE) - write.csv(InteractionScores_deletion_enhancers_and_Suppressors_L_lm,paste(outputpath,"ZScores_Interaction_DeletionEnhancers_and_Suppressors_L_lm.csv",sep=""),row.names=FALSE) - write.csv(InteractionScores_deletion_enhancers_and_Suppressors_K_lm,paste(outputpath,"ZScores_Interaction_DeletionEnhancers_and_Suppressors_K_lm.csv",sep=""),row.names=FALSE) + write.csv(InteractionScores_deletion_enhancers_L_lm,paste(outDir,"ZScores_Interaction_DeletionEnhancers_L_lm.csv",sep=""),row.names=FALSE) + write.csv(InteractionScores_deletion_enhancers_K_lm,paste(outDir,"ZScores_Interaction_DeletionEnhancers_K_lm.csv",sep=""),row.names=FALSE) + write.csv(InteractionScores_deletion_suppressors_L_lm,paste(outDir,"ZScores_Interaction_DeletionSuppressors_L_lm.csv",sep=""),row.names=FALSE) + write.csv(InteractionScores_deletion_suppressors_K_lm,paste(outDir,"ZScores_Interaction_DeletionSuppressors_K_lm.csv",sep=""),row.names=FALSE) + write.csv(InteractionScores_deletion_enhancers_and_Suppressors_L_lm,paste(outDir,"ZScores_Interaction_DeletionEnhancers_and_Suppressors_L_lm.csv",sep=""),row.names=FALSE) + write.csv(InteractionScores_deletion_enhancers_and_Suppressors_K_lm,paste(outDir,"ZScores_Interaction_DeletionEnhancers_and_Suppressors_K_lm.csv",sep=""),row.names=FALSE) - write.csv(Labels,file=paste("../Code/StudyInfo.csv"),row.names = FALSE) + write.csv(Labels,file=paste(StudyInfo),row.names = FALSE) print('ln 1570 write StudyInfo.csv after ') #write.table(Labels,file=paste("../Code/StudyInfo.txt"),sep="\t",row.names = FALSE) @@ -1697,7 +1708,7 @@ for(s in Background_Strains){ } } print("Pass Int ggplot loop") - write.csv(X_stats_interaction_ALL_final,paste(outputpath,"ZScore_Calculations.csv",sep=""),row.names = FALSE) + write.csv(X_stats_interaction_ALL_final,paste(outDir,"ZScore_Calculations.csv",sep=""),row.names = FALSE) @@ -1706,7 +1717,7 @@ for(s in Background_Strains){ Blank <- ggplot(X2_RF) + geom_blank() - pdf(paste(outputpath,"InteractionPlots.pdf",sep=""),width = 16, height = 16, onefile = TRUE) + pdf(paste(outDir,"InteractionPlots.pdf",sep=""),width = 16, height = 16, onefile = TRUE) X_stats_X2_RF <- ddply(X2_RF, c("Conc_Num","Conc_Num_Factor"), summarise, mean_L = mean(l,na.rm=TRUE), @@ -1863,7 +1874,7 @@ for(s in Background_Strains){ - pdf(paste(outputpath,"RF_InteractionPlots.pdf",sep=""),width = 16, height = 16, onefile = TRUE) + pdf(paste(outDir,"RF_InteractionPlots.pdf",sep=""),width = 16, height = 16, onefile = TRUE) X_stats_X2_RF <- ddply(X2_RF, c("Conc_Num","Conc_Num_Factor"), summarise, mean_L = mean(l,na.rm=TRUE), @@ -2155,7 +2166,7 @@ for(s in Background_Strains){ theme_Publication() - pdf(paste(outputpath,"RankPlots.pdf",sep=""),width = 18, height = 12, onefile = TRUE) + pdf(paste(outDir,"RankPlots.pdf",sep=""),width = 18, height = 12, onefile = TRUE) grid.arrange(Rank_L_1SD,Rank_L_2SD,Rank_L_3SD,Rank_K_1SD,Rank_K_2SD,Rank_K_3SD,ncol=3,nrow=2) grid.arrange(Rank_L_1SD_notext,Rank_L_2SD_notext,Rank_L_3SD_notext,Rank_K_1SD_notext,Rank_K_2SD_notext,Rank_K_3SD_notext,ncol=3,nrow=2) @@ -2262,7 +2273,7 @@ for(s in Background_Strains){ geom_hline(yintercept=c(-3,3)) + geom_point(size=0.1,shape=3) + theme_Publication() - pdf(paste(outputpath,"RankPlots_lm.pdf",sep=""),width = 18, height = 12, onefile = TRUE) + pdf(paste(outDir,"RankPlots_lm.pdf",sep=""),width = 18, height = 12, onefile = TRUE) grid.arrange(Rank_L_1SD_lm,Rank_L_2SD_lm,Rank_L_3SD_lm,Rank_K_1SD_lm,Rank_K_2SD_lm,Rank_K_3SD_lm,ncol=3,nrow=2) grid.arrange(Rank_L_1SD_notext_lm,Rank_L_2SD_notext_lm,Rank_L_3SD_notext_lm,Rank_K_1SD_notext_lm,Rank_K_2SD_notext_lm,Rank_K_3SD_notext_lm,ncol=3,nrow=2) @@ -2297,7 +2308,7 @@ for(s in Background_Strains){ get_lm_AUC <- lm(X_NArm$Z_lm_AUC~X_NArm$Avg_Zscore_AUC) AUC_lm <- summary(get_lm_AUC) - pdf(paste(outputpath,"Avg_Zscore_vs_lm_NA_rm.pdf",sep=""),width = 16, height = 12, onefile = TRUE) + pdf(paste(outDir,"Avg_Zscore_vs_lm_NA_rm.pdf",sep=""),width = 16, height = 12, onefile = TRUE) print(ggplot(X_NArm,aes(Avg_Zscore_L,Z_lm_L)) + geom_point(aes(color=Overlap),shape=3) + geom_smooth(method = "lm",color=1) + ggtitle("Avg Zscore vs lm L") + @@ -2329,7 +2340,7 @@ for(s in Background_Strains){ pgg <- ggplotly(lm_v_Zscore_L) #pgg - plotly_path <- paste(getwd(),"/",outputpath,"Avg_Zscore_vs_lm_NA_rm.html",sep="") + plotly_path <- paste(outDir,"Avg_Zscore_vs_lm_NA_rm.html",sep="") saveWidget(pgg, file=plotly_path, selfcontained =TRUE) X_NArm$L_Rank <- rank(X_NArm$Avg_Zscore_L) @@ -2357,7 +2368,7 @@ for(s in Background_Strains){ num_genes_NArm2 <- (dim(X_NArm)[1])/2 - pdf(paste(outputpath,"Avg_Zscore_vs_lm_ranked_NA_rm.pdf",sep=""),width = 16, height = 12, onefile = TRUE) + pdf(paste(outDir,"Avg_Zscore_vs_lm_ranked_NA_rm.pdf",sep=""),width = 16, height = 12, onefile = TRUE) print(ggplot(X_NArm,aes(L_Rank,L_Rank_lm)) + geom_point(aes(color=Overlap),shape=3) + geom_smooth(method = "lm",color=1) + ggtitle("Rank Avg Zscore vs lm L") + @@ -2481,7 +2492,7 @@ for(s in Background_Strains){ theme_Publication() - pdf(paste(outputpath,"RankPlots_naRM.pdf",sep=""),width = 18, height = 12, onefile = TRUE) + pdf(paste(outDir,"RankPlots_naRM.pdf",sep=""),width = 18, height = 12, onefile = TRUE) grid.arrange(Rank_L_1SD,Rank_L_2SD,Rank_L_3SD,Rank_K_1SD,Rank_K_2SD,Rank_K_3SD,ncol=3,nrow=2) grid.arrange(Rank_L_1SD_notext,Rank_L_2SD_notext,Rank_L_3SD_notext,Rank_K_1SD_notext,Rank_K_2SD_notext,Rank_K_3SD_notext,ncol=3,nrow=2) @@ -2588,7 +2599,7 @@ for(s in Background_Strains){ geom_hline(yintercept=c(-3,3)) + geom_point(size=0.1,shape=3) + theme_Publication() - pdf(paste(outputpath,"RankPlots_lm_naRM.pdf",sep=""),width = 18, height = 12, onefile = TRUE) + pdf(paste(outDir,"RankPlots_lm_naRM.pdf",sep=""),width = 18, height = 12, onefile = TRUE) grid.arrange(Rank_L_1SD_lm,Rank_L_2SD_lm,Rank_L_3SD_lm,Rank_K_1SD_lm,Rank_K_2SD_lm,Rank_K_3SD_lm,ncol=3,nrow=2) grid.arrange(Rank_L_1SD_notext_lm,Rank_L_2SD_notext_lm,Rank_L_3SD_notext_lm,Rank_K_1SD_notext_lm,Rank_K_2SD_notext_lm,Rank_K_3SD_notext_lm,ncol=3,nrow=2) @@ -2619,7 +2630,7 @@ get_lm_6 <- lm(X_NArm$Z_lm_AUC~X_NArm$Z_lm_r) L_lm_6 <- summary(get_lm_6) -pdf(file=paste(outputpath,"Correlation_CPPs.pdf",sep=""),width = 10, height = 7, onefile = TRUE) +pdf(file=paste(outDir,"Correlation_CPPs.pdf",sep=""),width = 10, height = 7, onefile = TRUE) ggplot(X_NArm,aes(Z_lm_L,Z_lm_K)) + geom_point(shape=3,color="gray70") + geom_smooth(method="lm",color="tomato3") + @@ -2707,24 +2718,7 @@ ggplot(X_NArm,aes(Z_lm_r,Z_lm_AUC)) + geom_point(shape=3,color="gray70") + theme_Publication_legend_right() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.text.x = element_text(size=16),axis.title.x = element_text(size=18), axis.text.y = element_text(size=16),axis.title.y = element_text(size=18)) - - - dev.off() - #write.csv(Labels,file=paste("../Code/Parameters.csv"),row.names = FALSE) timestamp() - -#BoneYard*********************************************** -#I'm thinking this parameter needs to be save somewhere "permanent' for the record so outputs can be reproduced. -#take this out of the Arguments. In Matlab I could for future in .mat file. Maybe I could save the SD Args[2] as part of the StudyInfo.txt. -#Corruptable but better than nothing. -#if(is.na(Args[2])){ -# std=3 -#}else { -# std= Arg[2] -#Delta_Background_sdFactor <- 2 #Args[3] -#DelBGFactr <- as.numeric(Delta_Background_sdFactor) -#} -