Data Analysis and Visualization Environment (DAVE) is our primary 3D visualization tool. DAVE help and some DAVE images are viewable online.
This program takes one or more outputs of a simulation (e.g. output of run1 or run2) For each z slice (time point) in each output: blurs it with a psf, pulls out the in-focus slice. Then calculates total concentration CHANGE (relative to time 0) as a function of distance from the calcium entry site (assumed to be the brightest point over all time after blurring). This is written to a file with the same root as the input simulation, but with a .xy extension. This can be viewed with xmgrace, or read into matlab (readxmgr) and visualized with acurves().
Assumes simulation was at 100 nm resolution. Uses a fixed psf with 300x300x150 nm voxels. Several scale factors are applied, so comparisons of magnitudes between files are not valid (but comparisons between time pts from the same file are fine).
analyze_simulation.pl [options] sim.i2i [additional input files]
options: -zrange=#:#:# startz stopz dz, inclusive, zero indexed. default is all z slices (i.e., all time pts). -verbose print out extra into to stderr as it goes along -keepfiles keep all intermediate files. may take lots of space!
examples: analyze_simulation.pl -verbose -zrange=0:10:1 run1_cafl.i2i > & sim.out
matlab: data = readxmgr('run1_cafl.xy'); acurves(data);
analyze_simulation.pl -verbose run1_cafl.i2i run1_ca.i2i > & sim.out
This program reads in a 2D diffusion tracks and calculates a D for each file (track_\d+.out) by fitting a line to the MeanSquaredDisplacement vs time:
MSD(t) = 4Dt (it would be MSD(t) = 4Dt + 2*sigma^2 if there is noise, see ref below).
so the slope of the line divided by 4 should yield D (even in the case with noise). note: Lag intervals overlap one another unless -no_overlap is specified. meanD, varD, and standard error of the meanD is printed to stderr. note: if -no_overlap is not specified, the standard error of the mean (seD) may not be correct since the points aren't independent.
-minlag=# if not specified, minlag = 1 -maxlag=# if not specified, maxlag = 1/2 the total time in the first track file analyzed -duration=# rather than analyze each track in its entirety, split it into (non-overlapping) sections of the specified duration (number of time points). -files=file1:...:filen list of filenames to analyze (instead of all track_#.out files)
-filerange=#:# only look at files track_#1.out to track_#2.out (inclusive)
-logfit also do the D calculation (via a line fit) in the log-log space (ln (MSD(t) = ln(4D) + ln(t)) -no_overlap
-print_Ds print all the Ds calculated, one per line (tracknum D), to stdout. -print_wDs print all the weighted Ds calculated, one per line (tracknum D), to stdout.
-print_graph prints to stdout MSD vs time -notracknum the filename being read does not have tracknumber in it, just autonumber them. -dummy just print some info to stderr -verbose which file is being analyzed printed to stderr
example: calc_diffusion3.pl -print_graph files=track_1.xy | cull -c -s -1 -2 |xmgrs
# take a 1,000 pt track and analyze it as 10 separate 100 length tracks to see what those Ds look like:
calc_diffusion3.pl -print_Ds -print_wDs -files=track_1.xy calc_diffusion3.pl -duration=100 -print_Ds -print_wDs -files=track_1.xy
calc_diffusion.pl same program but only operates on one file. calc_diffusion2.pl same program but without the -duration option SingleParticleTracking_DiffusionCoefficients_Saxton.pdf p. 1746 for the algorithm SubdiffusionOfSingleParticleTracks.pdf - p. 4 gives the formula if there is noise in the estimate.
This program combines data produced by cross_sections_of_voronoi (-vfile) with an analysis of the same image within fiji (via the stereology_lml macro, -sfile). Results printed (on each line there will be sfile data then vfile data) to stdout.
compare_volume_measurements.pl [options] -vfile=filename1.txt -sfile=filename2.txt
options: -print_volumes rather than print entire combined lines, just print volume info in columns. the columns are: region: the region id, this is the same as the pixel value. sarea: the cross-sectional region area as reported by the fiji stereology_lml.txt macro. svolume: the volume estimated by stereology_lml.txt swvolume: the weighted volume estimated by stereology_lml.txt (not really a true volume). true_volume: the true volume of the region (as created by cross_sections_of_voronoi). area: area as reported by cross_sections_of_voronoi volume_from_area another way to estimate volume from area, w/out stereology: ((1/.68)*area)^(3/2) -vprefix=# foreach field printed out, change its name by adding the specified prefix -sprefix=# foreach field printed out, change its name by adding the specified prefix this may be used so that field names become unique (otherwise, eg, there will be two "area" fields and two "volume" fields on each line). -verbose stuff to stderr
examples: cross_sections_of_voronoi -P voronoi.i2i > filename1.txt Fiji Edit->options->conversions - make sure "scale when converting" is NOT checked.
Plugins-> open i2i -> voronoi.i2i, then Image->type->32, Image->type->16, to get rid of the table lookup Plugins->macros->install... install stereology_lml.txt Plugins->macros->stereology_lml.txt -> voronoi_fiji.txt
compare_volume_measurements.pl -vfile=filename1.txt -sfile=voronoi_fiji.txt
# look at volume estimated via stereology (x axis) vs. true volume (y axis): compare_volume_measurements.pl -print_volumes -vfile=filename1.txt -sfile=voronoi_fiji.txt | cull -c -s -3 -5 | xmgrace -free -pipe
# look at volume estimated from area (x axis) vs. true volume (y axis): compare_volume_measurements.pl -print_volumes -vfile=filename1.txt -sfile=voronoi_fiji.txt | cull -c -s -7 -5 | xmgrace -free -pipe
# compare true volume (x axis) to both the stereology estimated volume and the volume from area measurements (y axis): compare_volume_measurements.pl -print_volumes -vfile=filename1.txt -sfile=voronoi_fiji.txt | cull -c -s -5 -3 -7 > tmp.nxy
xmgrace -nxy tmp.nxy -free
see also: /home/lml/krypton/Corvera/Raz/stereology/README, run_voronoi.pl cross_sections_of_voronoi /home/lml/krypton/packages/fiji/Fiji.app/macros/stereology_lml.txt
Reads in a file produced by makebeads2 (bead.info) and one produced by print_pixel_coords (mask.info). A bead's info is kept only if at least one of its voxels is under the mask. Kept beads are printed to stdout with their original values (bead id, position, etc.). Thus, if you actually mask the bead image (using maski2i) which was produced by makebeads2, this will keep the associated bead info for those beads which have at least one pixel left in the image after masking.
Note: this can only be run on data produced by versions of makebeads2 newer than 2/12/09.
filter_beads.pl bead.info mask.info > filtered_bead.info
-keepcenter=# : rather than keep a bead if any of its voxels is under the mask, only keep the bead
if its center is under the mask. # is the scale factor to convert from the original resolution the bead was created at (-S option to makebeads2) and the resolution of the mask. For example: makebeads2 -S 0.005 0.005 0.005 ... > beads.info If you then use a mask at .1 um you would specify -keepcenter=20
makebeads2 -v -d 1280 1280 60 -S 0.005 0.005 0.005 -n 300 -i 1 beads300.i2i > bead.info maskimage beads300.i2i mask.i2i beads_masked.i2i # now want to filter bead.info to leave only beads at least partially under mask.i2i (i.e., beads with at least one pixel # left in beads_masked.i2i) print_pixel_coords -2 mask.i2i > mask.info filter_beads.pl bead.info mask.info > filtered_bead.info
# now filtered_bead.info "matches" what is left in beads_masked.i2i
An example where resolution is changed: # voxel coords printed to bead.info by makebeads2 will reflect the new (chres'ed) coords appropriately:
makebeads2 -v -d 1280 1280 60 -S 0.005 0.005 0.005 -c 0.025 0.025 0.01 -n 300 -i 1 beads.i2i > bead.info
# now actually change the res of the bead image to match:
chres -S 5 5 2 beads.i2i beads_ch.i2i
# apply the mask to the lower res image. mask must already be at that same resolution:
maskimage beads_ch.i2i mask_ch.i2i beads_ch_masked.i2i
# get out coords of all pixels in this low res mask:
print_pixel_coords -2 mask_ch.i2i > mask_ch.info
# now just keep the bead info for beads which are under the mask at this low resolution:
filter_beads.pl bead.info mask_ch.info > filtered_bead.info
# now filtered_bead.info "matches" what is left in beads_ch_masked.i2i (I hope!)
See also: makebeads2, print_pixel_coords
Reads specified columns of data from an xpp output file. Sends output stdout.
note: time is always the first variable output (i.e., the first column).
graph_xpp2.pl [options] output.dat > name1_name2.dat
-allstates: extract all columns of data having to do with channel state.
This is all the basic data. Other data is derived from this.
-openstates: all open states. cannot be specified with -closedstates. -nopen: sum of all the open states. cannot be used with -closedstates. -closedstates: all closed states. cannot be specified with -openstates. -verbose
xpp -silent bk.ode -> output.dat graph_xpp.pl -vars=time:name1:name5:name7 output.dat > xpp_subset.out xmgrace -nxy xpp_subset.out
NOT DONE YET.
Calculates steady state oxygen levels as a function of distance from a blood vessel (see ~/krypton/Corvera/Olga). Writes to (rstar,Pstar) columns to stdout; this is the Krogh-Erlang equation. Equation 17 of KroghOxygenDiffusionModel1.pdf
see also: oxygen_krogh.c
Generate a random path (ie, diffusion).
random_path.pl [options] > track_1.out
options: -D=# diffusion coefficient, default = 1 -length=# pathlength (number of time points). default = 1000 -xy only print x,y rather than t,x,y
-method=oversample|gaussian oversample takes dt size substeps to calculate each 1 time unit step. gaussian just samples into an equivalent Gaussian distribution once for each 1 time unit step. default: oversample. gaussian option not tested yet... -dt=# substep time interval, default = .001 (only used if -method=oversample). -seed=# so random numbers are the same from run to run (only if -method=oversample)
note: using the default random number generator which is probably not great.
see also: run_random_path.pl calls random_path.pl to generate lots of paths. calc_diffusion3.pl can analyze these paths ~/krypton/diffusion/mytracks/
Reads in file produced by makebeads to check to see if sizes and amplitudes are Gaussian and Poission (-G and -P options). Produces _tmpsize and _tmpamp. You may want to use filterbeads.pl prior to using readbeads.pl. You can read these into xmgrace and then create histograms if so desired.
-skipzero skip values which are zero -outer put outer diameter (rather than "size") into _tmpsize. If makebeads used the -g option, then size always equals 1. outer should be the gaussian mean + 1/2 shell width. -rpts=# also prints out _coords.rpts file of coordinates of center of beads. This can be read in by addlines. Can also be read by pt_distances. # is a scale factor to divide by to get the new coordinates (so you can compare these to a lower resolution image). Set # to 1 if you don't want to scale the values. Note: coords are truncated an integer value (new = integer(old/# + 1)). skipzero command doesn't affect this option. note: if more than 1 point has the same new coordinates, then the w value for each duplicate point is 1 higher than the previous value (e..g, if -w=-10, then the first (5,6,3) has a w = -10, next has a value -9, next -8, etc. note: this also causes coords to be converted from 0-indexed (makebeads2 output) to 1-indexed (play assumes this for rpts files). -dontrescalez: when using -rpts=#, don't rescale the z coord -twod ignore z dimension when printing rpts file (still must specify -rpts option). -w=# in rpts file, set w to specified number, default = -1000 -amp use the amplitude of the bead for the w value in rpts file output
See also: filterbeads.pl, makebeads2, objs2dist121.pl, print_pixel_coords, filter_beads.pl
Generate lots of random paths files (ie, diffusion) by calling random_path.pl multiple times..
options: -seed=# so random numbers are the same from run to run -D=# diffusion coefficient, default = 1 -length=# pathlength (number of time points). default = 1000 -xy only print x,y rather than t,x,y -numpaths=# default = 10 paths -rootname=string default: track_ (so files will be called track_#.out) -verbose stuff to stderr
note: using the default random number generator which is probably not great.
see also: random_path.pl calc_diffusion3.pl
This script calls cross_sections_of_voronoi many times, and produces an image and output text file from each run.
options: -num=# number of times to run. default = 10 -D=# -D option to cross_sections_of_voronoi, average diameter, default = 150 -dummy don't do anything, just echo to stderr what would be done
compare the true volumes of voronoi regions (as reported by cross_sections_of_voronoi, e.g. voronoi_D170_0001_g0_G250_volume.txt) to those calculated within fiji using stereology_lml.txt (e.g., voronoi_D170_0001.txt)
Graphs volume (x) vs predicted volume (y) (or prints to stdout if -stdout is specified)
note: if looking at weighted volumes (eg, in weighted/ ), the individuals stereologic volumes
really are meaningless, and you should really just pull out the wmean values from each run.
voronoi_vs_fiji.pl [options] *volume.txt
-vdir=/path/to/directory location of directory with the output from cross_sections_of_voronoi (eg. voronoi_D170_0001.txt) default is the current directory
-debug info about all fields on each line printed to stderr
-stdout everything to stdout
-print_avg for each file, just graph mean and std error of the mean (gotten by combining all the data for each file, so may include multiple stereological runs) rather than all data pts note: this will average together the multiple fiji stereology runs which are often in the same *volume.txt file. I think that is ok (but may make the sem look better than it would be from just one run).
-use_wmean just use the mean or wmean stored in the *volume.txt file and compare these values to ??? ... not done yet
examples: cd /storage/big1/lml/Raz/voronoi run_voronoi2.pl -D=170 -num=20 # calls cross_sections_of_voronoi run_voronoi2.pl -D=220 -num=20 # calls cross_sections_of_voronoi Fiji # then run run_stereolgy_voronoi.ijm cd results/
voronoi_vs_fiji.pl -vdir=.. *D170*volume.txt # -> voronoi_vs_fij_D170.xmgr voronoi_vs_fiji.pl -vdir=.. *D220*volume.txt # -> voronoi_vs_fij_D220.xmgr
see also: /storage/big1/lml/Raz/voronoi/ - this is where I've done what is shown in the examples above /storage/big1/lml/Raz/voronoi/README, true_vs_predicted.xmgr /home/lml/krypton/Corvera/Raz/stereology/README, run_voronoi.pl, run_voronoi2.pl, compare_volume_measurements.pl, plot_results.pl cross_sections_of_voronoi /home/lml/krypton/packages/fiji/Fiji.app/macros/stereology_lml.txt /home/lml/krypton/packages/fiji/run_stereology.txt analyze_cross_sections_of_voronoi.pl
incorrect command line. This program creates an elliptical shell of brightness the specified brightness (-i option) It is used to create test images for Fischer/hmc.
Usage: ball [options] ball.i2i options: -d xdim ydim zdim: dimensions of image, default = (50,100,50) -r x y z: x,y, and z radius in voxels, default = (20.000000,20.000000,20.000000) -s samplerate: number of subsamples in x and z directions per pixel, default = 1 -i intensity: intensity of the ball's shell, default = 100.000000 -b intensity: intensity background voxels, default = 0.000000 -f image: read in the specified image, add the ball image to that (just replaces voxels along the ball's shell with new values). -t thickness: thickness of shell, default = 1.000000
Note: you can call this multiple times to get a shell within a shell
code in /home/lml/krypton/facil/ball.c
incorrect command line. Reads in 4 floating point images (or 4 short int images) as produced by simulation, calculates the difference image deltaca.i2i (as a floating pt image) between the observed free calcium and what would be expected in equilibrium conditions, i.e., [ca]-[predicted_ca], where [predicted_ca]
is kd*[CaFl]/[Fl]. Output concentration in nM.
All input concentrations should be in nM Each z slice is analyzed separately (i.e., z is assumed to be time).
Usage: calc_ca3 [options] ca.i2i fl.i2i cafl.i2i cafixed.i2i deltaca.i2i options: -k #: kd in uM, default = 1.125000
-h: print cumulative histogram for each z slice. options for printing graphs of free calcium (not concentration) vs. distance from ca entry location:
the units for the graph are nM*um^3, but 1nM = .6ions/um^3 , so multiply be .6 to get ions
(and 1pA = 3,000 ions/msec if you want to convert to pA-msec). -c: print delta calcium (delta from first slice) distance graph for each z slice. -n: normalize delta calcium distance graph for each z slice (so sum of amplitudes = 1), (otherwise when graph spacing is changed, -g option, the height of the graph will change). -a: print absolute calcium, not delta calcium. when this is applied to a uniform image it should produce a 4pi*r*r distribution. options for printing graphs of total calcium, i.e. ca+cafl+cafixed (not concentration): -C: print delta total calcium distance graph for each z slice. -N: normalize delta total calcium distance graph for each z slice (so sum of amplitudes = 1). -A: print absolute total calcium, not delta total calcium
-p # #: pixel size in microns in the radial (x) and lengthwise (y) directions, default = 0.10 and 0.10 -P # # #: pixel size in microns in x,y,z directions (rectangular not cylindrical coords), default = ( 0.10, 0.10, 0.10) -r: data is rectangular, not cylindrical coords (-P sets this too) -e x y: 0-indexed coords of calcium entry location (default: x= 0, y = ydim/2)
-g #: graph spacing for -c,-C,-n,-N,-a,-A options. defauls is pixel size. -s #: subsample each voxel this number of times in the x (radial) and y (longitudinal) directions, default = 3. Subsampling only applies to histogram and graph generation, not image generation. -z low hi dz: only look at z slices low to hi (0-indexed, inclusive), stepping by dz. The zdim of deltaca.i2i will not change, but skipped slices will be 0.
float2int, in matlab: readxmgr, acurves (see ~lml/krypton/JSinger for an example)
source code in ~lml/krypton/Walsh directory.
: in main, before MPI_Init()
: in main, after MPI_Init() : /home/lml/krypton/bin/cellsim -h help: unable to parse command: -h
This program simulates Ca and Na in a cylindrical cell with cylindrical symmetry. Note: It is recommended that the stderr output be examined for lines with the words: warning, Warning, error, or Error on them.
cellsim [options] >& output_datafile
Options controlling the duration of the simulation:
-time #: total time to simulate (in sec), default = 0.000000 -time_rate #: slow down time, ie take smaller dt steps during simulation. this may be necessary if you are creating a large ca gradient, for example, with -ca_ring or -ca_src options (.5 means 1/2 the dt). The number of iterations needed for convergence is roughly: I= (T*diff)/(dr*dr), where T is time, dr is annulus thickness and diff is diffusion constant. This ignores ion influx,etc. If you set diffusion rates very small, the iter # calculated will be too small if other stuff is happening, use -time_rate to compensate.
Options controlling the geometry of the model:
-radius #: radius of the cell in microns, default = 6.000000 -rings #: radius of the cell in rings, default = 60 Ring thickness = cell radius/ # of rings. -length #: length of the cell in disks, default = 300 Each disk will be the same thickness as a ring. -bc Rmax Lmin Lmax: boundary conditions at maximum radius and minimum and maximum length coordinate (i.e., at the cell ends). These must be the letters I (for infinite) or S (for solid, i.e. no ions go through the membrane). The I option is not always accurate. default: S -mask image.i2i: only voxels in image.i2i > 0 are part of the cell. Use this to specify an arbitrary shape for the cell. Not extensively tested. -spacing file: file specifies the spacing along each axis (i.e, size of each voxel). spacing is in nm, one entry per voxel, comment lines ignored (#). The cell direction (X,Y, or Z) must be the first letter on the line. Other characters on the line will be ignored. The number of entries must be < 2000 each direction. Do not specify -rings,-radius,and -length if you use this option. NOTE: the simulation time step goes as (length)^2, where length is the length of the shortest side of all the voxels specified. Any number of entries per line. First R, then D (length). Example: # R direction (10 pixels in this direction, high res at low R) X 25 50 75 100 200 300 400 500 600 1000 Y (lengthwise) direction (9 pixels in D, very high res in the middle) 100 200 300 400 500 540 550 560 1000 Use cyl2cart with -spacing option to convert output images to a uniform spacing. (I haven't added that option yet) Not extensively tested.
Options specifying substances (ions, buffers, compounds) present in the model (10 max):
-ion name free_concentration D valence: concentration in mMolar, D in cm^2/sec, valence in absolute value. valence is only used if -ion_src is also specified. concentration is the initial free concentration. -buffer name total_concentration D: concentration in mMolar, D in cm^2/sec. concentration is the initial total concentration (sum of its free and all bound forms). -kin bufname ionname boundname onrate offrate D: kinetics for interaction between bufname and ionname. Boundname is the name to give the bound pair. bufname and ionname must have been specified (with -ion -buffer or -kin options) prior to this option. onrate:in 1/(mMolar*seconds), offrate: in 1/seconds, D: in cm^2/sec -equil_error relerr iter: maximum relative equilibrium error allowed during the initial calculation of concentrations from the user specified free ion concentrations and total buffer concentrations. rel_error = ABS((true_kd - calculated_kd)/true_kd). Default = 1e-05 This should be small enough unless you have huge buffer concentrations (>100 mM?). iter (an int) is maximum number of iterations