Subsetting data with string matching in R

Getting a subset from a larger pool of data is a major strength in R. This subsetting can be based on numbers or strings; the latter being more challenging. String matching (using grep, grepl, sub) and subsetting (using, well, subset in R) are two separate feature and one can combine them to great effect. Details of what grep, grepl etc can do are at
http://stat.ethz.ch/R-manual/R-devel/library/base/html/grep.html . This site
http://atomicules.co.uk/2010/10/09/invert-a-subset-selection-when-using-grepl-in-r.html
was of great help when I wanted (had) to use grepl in subset but I also wanted to use the invert feature available only in grep and not in grepl. Turns out grepl returns a vectors of T and F which can be easily inverted with a simple !grepl(...)

Area under a curve

R has a neat way to calculate area under a curve (whether a function/equation or a trend line). The example below describes the two ways in which this can be achieved. rollmean requires library(zoo).

Making and filling grids with R

I had a table with 64 rows and 130 columns and I needed to map a third dimension as a simple X on the 2D plot based on some criteria. I got this done in R using a code similar to the one below.

Contact Maps using VMD

Contact maps are a quick way to identify residues of a protein (or monomer) that interact with residue from another protein, ligand or monomer. The script below helps obtain a distance list:

set seg1 [atomselect top "segname SG03 and name CA"]

set seg2 [atomselect top "segname SG04 and name CA"]

set file [open "Contact_map.dat" w]

set list1 [$seg1 get index]

set list2 [$seg2 get index]

foreach atom1 $list1 {

        foreach atom2 $list2 {

                set index1 [atomselect top "index $atom1"]

                set index2 [atomselect top "index $atom2"]

                set resid1 [[atomselect top "index $atom1"] get resid]

                set resid2 [[atomselect top "index $atom2"] get resid]

                set resnm1 [[atomselect top "index $atom1"] get resname]

                set resnm2 [[atomselect top "index $atom2"] get resname]

     puts $file "$resnm1 $resid1 $resnm2 $resid2 [veclength [vecsub [measure center $index1] [measure center $index2]]]"

                $index1 delete

                $index2 delete

        }

puts $file " "   # I include this line to make it easier to make contour plots using gnuplot.

}

close $file

Aligning hexameric proteins along an axis

This tcl script aligns a hexamer such that its central axis is aligned along the y-axis.


set all [atomselect top protein]
$all moveby [vecscale [measure center $all] -1]

set res10 [atomselect top "resid 10"]
set res40 [atomselect top "resid 40"]

set cenres10 [measure center $res10]
set cenres40 [measure center $res40]

set vector [vecsub $cenres10 $cenres40]

$all move [transvecinv $vector]
$all move [transaxis z 90]

foreach i {A B C D E F} {
        set chain [atomselect top "chain $i and protein and resid 4 to 105"]
        $chain writepdb chain$i.pdb
}

This script is also applicable to trimers, tetramers and other aggregates that have a well defined central axis.

Curve fitting in R

This is a simple example of a linear fit where a perfect fit is distorted with the help of a gaussian function. This data is then fit to a linear model. The most important factor determining the success of a fit is the model equation and the second most important factor is suitable starting values of the fit variables (a and b in this case). Further, if you assign the output to a variable called fit, values a and b can be extracted using:

fit < - nls
coef(fit)[1]
coef(fit)[2]

Saving data frames in R

If 'df' is your data frame,
save(df,file="yourfile.rdf",ascii=TRUE)
save the data frame in to a .rdf file which is readable. Although the file as lot more data than you might want, loading the file back into R irons out any wrinkles.
load("yourfile.rdf")
http://stackoverflow.com/questions/8345759/how-to-save-a-data-frame-in-r

Edit Jan 2014
If however you really want to avoid additional lines in the saved file, try the following:
write.table(df,filename,row.names=FALSE,col.names=FALSE)
This should produce a clean file with your data in a neat column.

Speaker and Mic settings on Ubuntu Linux

Although Skype and Google video get installed fairly easily on Ubuntu, bad sound settings could be a show stopper. Alsamixer controls the sound and mic and can be force rebooted using:

alsa force-reload
or
alsa reload

This usually fixes most problems. Other useful links:
http://askubuntu.com/questions/164518/hard-resetting-alsa-configuration
http://chriscarey.com/wordpress/2008/07/01/reset-alsa-sound-on-linux/
http://ubuntuforums.org/showthread.php?t=1068612
http://forums.bodhilinux.com/index.php?/topic/3800-solved-how-does-a-user-restart-alsa-without-rebooting/

Mapping 3d data on a 2d color map

Gnuplot has a very easy and straightforward way to convert a x-y-z data into a 3d plot or a 2d color map. Let's write the data in the form;
1 1     0.0
1 2     0.4451
1 3     0.4249

2 1     0.13142
2 2     0.0
2 3     0.983

3 1     xxx
3 2     xxx
3 3

A space between chunks of data is important. Using commands below you can get a 2d plot;

set size square
set pm3d map    # Remove the word 'map' to get a 3d plot
set xrange [1:261]  # Change the range to suit your needs
set yrange [1:261] 
splot 'all.dat'

Reducing the file size of PDF documents

gs -sDEVICE=pdfwrite -dCompatibilityLevel=1.4 -dPDFSETTINGS=/screen -dNOPAUSE -dQUIET -dBATCH -sOutputFile=output.pdf input.pdf

I recently made a nice pdf figure using R but the file was 18Mb heavy. With 7 such figures my manuscript was easily crossing the 100Mb mark which felt rather unnecessary. The command above scales down the18Mb file to a tenth of its size. Since these figures are based on vectors, scaling down the file size has little to no impact on the quality of the figure when zoomed in. gs is turning out to be far more potent than I initially thought.

Sweave, R and Latex

If you use R for your statistical analysis and graph plotting needs and then use those figures in your Latex documents, you will find this interesting. Sweave is a package in R which lets takes in an .rnw file (which contains usual text from a tex document with R commands embedded in it), compiles it (in R) and gives out a tex file. This tex file can in turn be compiled by pdflatex (if you have figures) or plain latex to get a properly formatted pdf/ps/dvi file with relevant R output. There are many excellent examples on the interwebs and I strongly recommend using Sweave right away.  The beauty of this approach is that on can save a lot of time by not spending on choosing the appropriate dimensions of a a png or an eps output. Sweave does it for you.

Edit: If you want to insert a R-generated figure then add this to your tex file:

\begin{figure}
\begin{center}
<>=
out3 <- br="" i="" lm="" x="" y="">plot(x, y)
@
\end{center}
\end{figure}


If you want to insert text from R but properly tabbed and with all the ampersands added, then add these lines to your tex file:

<>=
library(xtable)
print(xtable(foo, caption = "ANOVA Table", label = "tab:one",
         digits = c(0, 0, 2, 0, 2, 3, 3)), table.placement = "tbp",
         caption.placement = "top")
@

Manipulating dihedral angles in VMD

VMD's scripting interface is quite powerful. It has a trans bond command which gives the transformation matrix required to achieve desired dihedral rotation. If you want to edit the dihedral 1-2-3-4 (i.e. rotate around bond 2-3), then you first select the set of atoms on which the transformation is to be applied; 1,2,5,6,7,8,9 in this case. To this selection you apply the transformation;

$moveselection move [trans bond [lindex [$a1 get {x y z}] 0] [lindex [$a2 get {x y z}] 0] 10 deg]

trans bond expects only the atom indices. The value 10 indicates the magnitude and sign of change to the dihedral 1-2-3-4. deg says the change should be in degrees and not radians. To get a nice picture of rotation around a bond use this script:

set movesel [atomselect top "index 1 2 5 6 7 8 9"]
for {set i 0} {$i <360 } {incr i 10}
{
$movesel move [trans bond [lindex [$a1 get {x y z}] 0] [lindex [$a2 get {x y z}] 0] 10 deg]
render TachyonInternal [format "rotate%03d.tga" $i]
}

Joining multiple pdf files in Linux

gs -dBATCH -dNOPAUSE -q -sDEVICE=pdfwrite -sOutputFile=finished.pdf file[1-3].pdf

This is a great way to join many pdf files into one file. Mac's Preview has another way to do the same but in a much easier drag-and-drop fashion. The line below in my .bashrc file helps automate the concatenation process.

alias join='gs -dBATCH -dNOPAUSE -q -sDEVICE=pdfwrite -sOutputFile=awkwardName.pdf'

time and script

These are great commands in Linux while debugging and timing codes.

Combining figures into one PDF file

convert -compress JPEG *.jpg output.pdf
convert -compress Zip *.jpg output.pdf

Both the commands join a bunch of jpeg files into one pdf file. This came in handy recently. Dropping the compress flag and its options gave a smaller file. Might look into it later.

convert *.jpg output.pdf

Bibtex being bibtex

Do you notice any difference between this:


@article{Dashnau2007Computational,
    abstract = {.....}
    author = {Dashnau, J. L. and Vanderkooi, J. M.},
    journal = {J. Food Sci.},
    month = jan,
    number = {1},
    pages = {R001--R010},
    title = {{Computational Approaches to Investigate How Biological Macromolecules Can Be Protected in Extreme Conditions}},
    volume = {72},
    year = {2007}
}

and this:


@article{Dashnau2007Computational,
    abstract = {.....}
    author = {Dashnau, J. L. and Vanderkooi, J. M.},
    journal = {J. Food Sci.}
    month = jan,
    number = {1},
    pages = {R001--R010},
    title = {{Computational Approaches to Investigate How Biological Macromolecules Can Be Protected in Extreme Conditions}},
    volume = {72},
    year = {2007}
}

Well, I couldn't either. Notice the comma at the end of the journal entry? Well without the comma Bibtex ignores the ensuing metadata creating problems when compiling your document. Something to be careful about.

Batch cropping images using Mogrify

The command 'mogrify' can crop multiple images simultaneously. Sample case is below.
mogrify -crop WidthxHeight+x+y        images*.png
mogrify -crop 2796x983+429+327        images*.png
The values of width, height, x and y can be obtained by opening an image in Gimp and selecting the rectangle of interest.

Firefox/Thunderbird already running but not responding

Firefox and Thunderbird can sometimes give this error when you restart a linux box or when you try to make multiple log ins into a shared file system. The result is that Firefox/Thunderbird will not open. The causes for this are not entirely clear to me but I found a solution. Simply delete the .parentlock file in both cases and restart the service. Locations for this file are as follows:
Firefox: ~/.mozilla/firefox/xxxxx.default/.parentlock
Thunderbird: ~/.thunderbird/xxxxx.default/.parentlock
Since both the folders and files start with a dot, they are not visible with a simple
ls -lrt
You have to use
ls -lrta

Emacs Primer

I decided just a little while ago that it was time I learned how to use a text editor other than vi. Below are my notes as I keep learning more about Emacs.

Open a new file: emacs newfile.txt &

Starting writing in the file.

Save: ctrl x s
Search Forward: ctrl s
Search Reverse: ctrl r

Exit emacs: ctrl x c

Select Region:      ctrl space and arrows
Mark Activate:      ctrl space
Mark Deactivate: ctrl space or ctrl g
Select All: ctrl x (wait) h

Delete Region: ctrl w
Copy Region:   Esc  w
Paste buffer: ctrl y

Delete till end of sentence: ctrl k

Delete line: ctrl S backspace

Delete next word: Esc d

Undo last action: ctrl _ (underscore)

Jump to end of file: Esc >
Jump to start of file: Esc <
Go to line 10: Esc x 10   

List all buffers: ctrl x (leave ctrl) b

Split screen:                                   ctrl x (leave ctrl) 2
Close all but current buffer:     ctrl x (leave ctrl) 1
Close only the current buffer:  ctrl x (leave ctrl) 0

Repeat 9 times: ctrl u (leave ctrl) 9
Type - 9 times: ctrl u (leave ctrl) 9 -

If you start a command and then decide you don't want to complete it, type C-g to cancel it.

http://doors.stanford.edu/~sr/computing/emacs.html
http://www.nbcs.rutgers.edu/newdocs/unx01101/unx01101.php3
http://dept.cs.williams.edu/~kim/cs136/s04/emacs.html

Edit: Dec 28, 2012
I was writing a manuscript in Emacs and the bibtex reference key had a beta symbol which meant that I had to type this into Emacs;   Kheterpal2000Aβ and not Kheterpal2000$\beta$. To do this I copied the following two lines into ~/.emacs and solved my problem. Thanks to this site.

(global-set-key (kbd " a") "α") ; F9 followed by a
(global-set-key (kbd " b") "β")

Java script to get Long/Lat in Google Maps

javascript:alert(window.gApplication.getMap().getCenter());

Latex Template

Reminding myself of every detail in the assigning of packages and page type in Latex seems like a bit of a chore. So I copied the following from a working .tex file to reduce redundancy.

\documentclass[12pt,doublespace]{article}

\usepackage[top=3cm,bottom=3cm,left=2.5cm,right=2.5cm]{geometry}

\usepackage{rotating}

\usepackage{apacite}

\begin{document}

\title{}

\author{

Name

\small Affiliation

\and

New Name

\small Affiliation

}

\maketitle

\begin{abstract}

...

\end{abstract}

\section{}

Main body

\end{document}

pdflatex vs pdftex

I have a document with lot of .eps figures embedded in it. These figures were of high quality but also large size. To cut down on the size I converted them to .png which were a tenth of the .eps format. Now using
latex file.tex (with eps figures)
worked well. But to use .png figures I thought I had to use pdftex.
pdftex file.tex (with png figures)
Turns out I should have used pdflatex instead. I ended up losing quite some time on this. The switch to pdflatex did the job and gave me a file which was a tenth of the original size. More on this from here:

"pdftex expects the file to be plain TeX. pdflatex expects it to be LaTeX. These
are two different dialects of TeX. You do not expect a C compiler to work with
Java, do you?

Note that in practice pdftex and pdflatex are often the same command, which
analyzes how it was called to switch to the proper input language. On some
systems gcc behaves in the same way, compiling C if it is called as gcc or
Fortran if it is called as gfortran."

Bottom line, figures can be heavily compressed by converting them to png format w/o a large loss in clarity; and pdflatex can handle png figures. Even if there is a loss of clarity, this is usually just fine for supporting information; maybe not so much for main article.

Enough said...

Handling references using BibTex

Scientific journals with arbitrary requirements for reference styles can be a source of great heart-ache. When writing articles using Latex-BibTex, ability to handle reference styles can be useful. These articles are a great starting point to learn more about this:

http://amath.colorado.edu/documentation/LaTeX/reference/faq/bibstyles.html#styles

http://chenfuture.wordpress.com/2007/09/24/diy-your-bibtex-style-file/

http://www.math.washington.edu/tex-archive/macros/latex/contrib/custom-bib/merlin.pdf

Controlling axis properties in R

http://www.programmingr.com/content/controlling-margins-and-axes-oma-and-mgp/

The mgp option
In addition to changing the margin size of your charts, you may also want to change the way axes and labels are spatially arranged. One method of doing so is the mgp parameter option. The mgp setting is defined by a three item vector wherein the first value represents the distance of the axis labels or titles from the axes, the second value is the distance of the tick mark labels from the axes, and the third is the distance of the tick mark symbols from the axes. As with the oma option discussed above, the distances are given in line widths. The defaults for the mgp setting are c(3, 1, 0). The examples below illustrate the effects of changing the various mgp values. Note: the mgp.axis() function in the Hmisc package can be used to change these settings for each axis individually.

ImageJ Resources

http://rsbweb.nih.gov/ij/
http://www.fmhs.auckland.ac.nz/sms/biru/facilities/analysis_resources.aspx

Learning ImageJ

ImageJ is a fantastic light-weight tool to analyze and manipulate images. I needed to use it for quickly analyzing microscopy images. Here I will list features of ImageJ that I found particularly useful for my purposes:
- I took a bunch of images at various gains to find out the point where bleaching starts occurring. So I have a series of images in which the pixel intensity varies from 0 to 255.

ImageJ allows me to load the images as a stack and calculate the maximum pixel value in each image by selecting the whole image and hitting 'm' (for measure). Plotting the maximum value gives this:

So, bleaching occurs somewhere between a gain of 450 and 475.

Print page-range using lpr

I use the extremely light-weight pdf viewer, xpdf. For some unknown reason, I cannot print pages from xpdf. A nice workaround I found was with dear old lpr which can print a range of pages thus;
lpr -P HP_LaserJet_6MP file.pdf -o page-ranges=2-5
lpstat -a # Gives the status (and names) of all the connected printers

Batch processing images in Gimp

Over the past few weeks I've had to apply a specific set of changes to a large number of images. Doing this repeatedly took a great toll on my soul and Gimp came to my soul's rescue. Take a read at this and this and help yourself.

Memory convservation in R

Depending on how large a data set you are handling, R can take up large swaths of your RAM and slow down other processes. It is therefore a good idea to remove unused variables.
ls()  # Gives the list of variables/arrays currently in use.
rm(list=ls()) # Removes all the variables and frees up memory.
Useful links:
 http://www.matthewckeller.com/html/memory.html

Benchmarking Simulations

It is always a good idea to benchmark your simulations and select the optimal number of cores for your simulations. This exercise greatly reduces wastage of resources and time. The script below helps with the setup. You need files conf.in and template.pbs

for i in 24 48 96 192 240 576 960 1200; do
        mkdir scaling$i
        cp conf.in scaling$i
        sed 's/xxx/'$i'/g' template.pbs > scaling${i}/sub.pbs
        echo "cd scaling$i";
        echo "aprun -n $i $namd conf.in 2>&1 1> conf.out &";  
        echo "cd \$PBS_O_WORKDIR";
done

exit;



----- Template.pbs ------
#PBS -q debug
#PBS -l mppwidth=xxx
#PBS -l walltime=00:30:00
#PBS -N scaling
#PBS -S /bin/bash
#PBS -j oe

namd=/usr/common/usg/namd/2.8/bin/namd2
cd $PBS_O_WORKDIR
aprun -n xxx $namd conf.in 2>&1 1> conf.out

Making 2D heat-maps of protein secondary structure

I needed to make 2D heat-maps of protein structure with time on x-axis and protein residues on y-axis. This helps to see the time evolution of structure quite nicely. This is easily accomplished by loading the trajectory in VMD and sourcing the file ss_traj.tcl and going through your entire trajectory once. This creates a buffer with all the structure information. ss_traj.tcl creates a ss_traj.dat file which has the structure information. This ss_traj.dat file is then used by ss_plot.c, after compilation of course, to generate a .ps file with the desired output. All this is explained quite well here. I made a personal copy of the files here:
ss_traj.tcl
ss_plot.c
Sample ss_traj.dat file.

Once you have the output in .ps format you can do interesting things to change the file. The 7th line has the scaling information. The width and height of the output were not to my satisfaction; so I changed the scaling factors to adjust the display. Also the executable could crash when running due to memory requirements. Simply change the dimensions of the array in ss_plot.c according to your needs.

Reading/processing multiple files in R

j<-1;
for (i in seq(300,390,15)) {
df<-read.table(paste("sasa",i,".txt",sep=``));
# Files are named "sasa300.txt" for ex.
tab[j,]<-cbind(i,mean(df[,2]),sd(df[,2]));
print(tab[j,]);
j<-j+1;
}
> tab
     [,1]     [,2]      [,3]
[1,]  300 3408.514  96.48097
[2,]  315 3422.388 119.14861
[3,]  330 3391.590  81.57252
[4,]  345 3424.176  84.84211
[5,]  360 3394.205 100.37193
[6,]  375 3462.469 122.60486
[7,]  390 3438.328  92.70519

Quick and dirty way to plot error bars in R

If you have data that looks as follows, a one line command can draw error bars to the plot.
-x---T-----y-------------error-----
p5 315 2.513402 0.5711628                                                      
p5 330 2.173868 0.4711594                                                      
p5 345 2.132954 0.634045                                                       
p5 360 2.092686 0.4851832                                                      
p5 375 2.966971 0.9165328                                                      
p5 390 2.214886 0.4288291

r<-read.table("data.txt")
plot(r[,2],r[,3],ylim=c(0,4))
lines(r[,2],r[,3])
for(i in 1:length(r[,1])) {segments(r[i,2],r[i,3]+r[i,4],r[i,2],r[i,3]-r[i,4])}

I am drawing a simple line from (x,y+err) to (x,y-err). The result is not as fancy as other efforts but it is enough for my purposes.

Using psfgen to delete unwanted atoms in simulations

psfgen's user manual can be found here.

Adding labels to figures using ImageMagick

This is a 1200x1200 pixel blank file to which I added labels as an exercise in using the annotate flag in ImageMagick's convert.
 
convert blank.gif -pointsize 50 \
-gravity south -append -annotate +000+000 '0s' \
-gravity south -append -annotate +100+100 '1s' \
...
-gravity south -append -annotate +900+900 '9s' \
-gravity north -append -annotate +000+000 '0n' \
-gravity north -append -annotate +100+100 '1n' \
...
-gravity north -append -annotate +900+900 '9n' \
-gravity northwest -append -annotate +000+000 '0nw' \
-gravity northwest -append -annotate +100+100 '1nw' \
...
-gravity northwest -append -annotate +900+900 '9nw' \
-gravity northeast -append -annotate +000+000 '0ne' \
-gravity northeast -append -annotate +100+100 '1ne' \
...
-gravity northeast -append -annotate +900+900 '9ne' \
out.png

Joining images using Latex

Do not forget to use pdflatex (not pdftex) to compile tex files which contain non-eps figures.
 
\documentclass{article}
\usepackage{rotating}
\usepackage{graphicx,subfig}

% Default margins are too wide all the way around. I reset them here
\setlength{\topmargin}{-.5in}
\setlength{\textheight}{9in}
\setlength{\oddsidemargin}{.125in}
\setlength{\textwidth}{6.25in}

\begin{document}


\begin{sidewaysfigure}
\begin{tabular}{cc...c}
\includegraphics[scale=0.10, bb= 0 0 669 834]{view1_mol0} &
\includegraphics[scale=0.10, bb= 0 0 669 834]{view1_mol10} \\

\includegraphics[scale=0.10, bb= 0 0 669 834]{view0_mol0} &
\includegraphics[scale=0.10, bb= 0 0 669 834]{view0_mol10} \\

\includegraphics[scale=0.10, bb= 0 0 669 834]{view2_mol0} &
\includegraphics[scale=0.10, bb= 0 0 669 834]{view2_mol10} \\

\end{tabular}
\end{sidewaysfigure}

\end{document}

Joining images using ImageMagick's convert

convert +append list_of_figures out_1.png
convert +append list_of_figures out_2.png
convert +append list_of_figures out_3.png
convert -append out_1.png out_2.png out_3.png final_out.png

+append puts images side by side while -append adds them top down.
convert -border 2x2 in.png out.png
This puts a light border around images which comes in handy when creating a grid of images from a bunch of individual images which have no border.

Empty plots in R

When plotting a grid of graphs in R, I needed leave some cells empty. For this the following commands come handy;

plot(0,xaxt='n',yaxt='n',bty='n',pch='',ylab='',xlab='')

plot.new()

grid.newpage()

grid()

Subsetting in R

Following are some example of subsetting in R I have used recently.


> head(crystal)
    V1     V2      V3      V4
1 1ACJ PHE120 290.980 289.101
2 1ACJ PHE153 307.673 283.814
3 1ACJ PHE155 296.153  89.152
4 1ACJ PHE186 311.606 292.186
5 1ACJ PHE187 298.920 270.613
6 1ACJ PHE197  66.678 323.351

>subset(crystal, crystal$V2=="PHE330" & grepl("^2A",crystal$V1))

http://www.ats.ucla.edu/stat/r/faq/subset_R.htm
http://stackoverflow.com/questions/2125231/subsetting-in-r-using-or-condition-with-strings

Charmm2lammps

Charmm2lammps.pl is a useful little code which takes in forcefield [top/par], psf and pdb files in Charmm format and converts them to Lammps data and sample input file. This is extremely useful for running simulations in the much better written MD code, LAMMPS. While trying to simulate a polymer chain I came across a very interesting feature in charmm2lammps; it really expects the force constant in the forcefield parameter file to be non-zero. This works just fine when you never want to specify angle/dihedral parameters to free chains. But when you want to control the parameters in between runs, it is a good idea to define all the angle and dihedrals and then decide to turn them off/on in LAMMPS by setting the force constant to 0/nonzero.

Compiling C programs on Ubuntu

On Ubuntu, it turns out, this does not work
gcc -lm code.c
but this does
gcc code.c -lm

Running job on clusters continuously

Here is a way to run jobs on cluster in a continuous fashion. This example addresses NAMD jobs but should also work for others.

Labeling figures

http://www.imagemagick.org/Usage/annotating/
This website offers a  nice way to add labels to images that "convert" can handle. I thought it was a nice idea to have my name on figures I generate.

Public folder on Dropbox

http://dl.dropbox.com/u/65895350/

Extracting coordinates of protein from a simulation

set dir [pwd]
mol new     ${dir}/ache.psf       type {psf}
mol addfile ${dir}/ache.pdb       type {pdb} waitfor all
mol addfile ${dir}/acheminnpt.dcd type {dcd} waitfor all

set molnum 0
set skip 1
set backframe0 [atomselect $molnum backbone frame 0]

for {set i 0} {$i < [molinfo $molnum get numframes]} {incr i $skip} {
set pro  [atomselect $molnum protein  frame $i]

# Change 'protein' to "resid N to Z" to suit your needs.
set back [atomselect $molnum backbone frame $i]
set trans_mat [measure fit $back $backframe0]
$pro move $trans_mat
set k [expr $i/$skip]
$pro writedcd ${dir}/dcd/project-pdb-$k.dcd
}

Now generate a protein-only psf file. At this stage you should have n number of dcd files. To join them all in to one dcd file,

mol new    ${dir}/ache_protein.psf
set pro_frames [expr [molinfo $molnum get numframes]/$skip ]
for {set j 0} {$j < $pro_frames} {incr j 1} {
mol addfile ${dir}/dcd/project-pdb-$j.dcd
}
set top [molinfo top]
animate write dcd {${dir}/dcd/project_protein.dcd} beg 0 end [expr $pro_frames - 1] skip 1 [molinfo top]

The last step is where the individual files are stitched together.

Installing Gimp on Ubuntu 11.1

Last month I installed Gimp on Ubuntu 11.1 which dropped Gimp since they deemed it too complex for new linux converts. Since this was last month I don't remember all the steps I took but I did copy the commands to a file and they might come handy in future. Here they are...

   69  sudo add-apt-repository ppa:matthaeus123/mrw-gimp-svn
   71  sudo apt-get install gimp gimp-data gimp-plugin-registry gimp-data-extras
   73  sudo add-apt-repository ppa:lucid-bleed/ppa
   74  sudo add-apt-repository ppa:matthaeus123/mrw-gimp-svn
   79  mkdir gimp
   80  cd gimp/
   81  wget ftp://ftp.gimp.org/pub/gimp/v2.7/gimp-2.7.2.tar.bz2
   82  tar xjf gimp-2.7.2.tar.bz2
   83  export PATH=/opt/gimp-2.7/bin:$PATH
   84  export PKG_CONFIG_PATH=/opt/gimp-2.7/lib/pkgconfig
   85  export LD_LIBRARY_PATH=/opt/gimp-2.7/lib
   86  cd gimp-2.7.2/
   87  ./configure —prefix=/opt/gimp-2.7
   89  ./configure
   90  cd ..
   92  sudo apt-get install git
   93  git clone git://git.gnome.org/babl
   94  cd gimp-2.7.2/
   95  ./configure --prefix=/opt/gimp-2.7
   96  cd ..
   97  ./autogen.sh --prefix=/opt/gimp-2.7
   99  cd babl/
  100  ./autogen.sh --prefix=/opt/gimp-2.7
  101  make -j5
  102  cd ../gimp-2.7.2/
  103  ./configure --prefix=/opt/gimp-2.7
  104  cd ../babl/
  105  cd ..
  106  git clone git://git.gnome.org/babl
  107  ./autogen.sh --prefix=/opt/gimp-2.7
  108  sudo apt-get install libtool
  109  ./autogen.sh --prefix=/opt/gimp-2.7
  110  cd babl/
  111  ./autogen.sh --prefix=/opt/gimp-2.7
  112  make -j5
  113  sudo make install
  114  cd ../gimp-2.7.2/
  115  ./configure --prefix=/opt/gimp-2.7
  116  cd ..
  117  wget ftp://ftp.gimp.org/pub/gegl/0.1/gegl-0.1.6.tar.bz2
  118  tar -xjf gegl-0.1.6.tar.bz2
  119  sudo apt-get install libjpeg62-dev libopenexr-dev librsvg2-dev
  120  cd gegl-0.1.6/
  121  ./configure --prefix=/opt/gimp-2.7
  122  sudo make install
  123  cd ..
  124  ls
  125  cd gimp-2.7.2/
  126  ./configure --prefix=/opt/gimp-2.7
  127  sudo apt-get install libtiff4-dev python-dev python-gtk2-dev
  128  ./configure --prefix=/opt/gimp-2.7
  129  which gimp
  130  sudo apt-get install libexif-dev liblcms1-dev
  131  ./configure --prefix=/opt/gimp-2.7
  132  make -j5
  133  make install -j5
  134  sudo make install -j5
  135  ls /opt/gimp-2.7/bin/gimp-2.7
  136  /opt/gimp-2.7/bin/gimp-2.7
  137  history > ~/gimpinstallcommands.dat

Making movies with VMD/PyMol

These two scripts help greatly in making and handling movies. I need scripts instead of the GUI because most of my data is on machines on the other coast and it gets rather cumbersome to pull a gigantic trajectory file just to see a quick snapshot of the system. I could make PDB (or better dcd) files of parts of trajectory I am interested in, but I will need to make movies at some point.

view_change_render.tcl
trajectory_movie.tcl
http://www.myofilament.org/documents/howto/modeling/VMD_Tcl.htm

PyMol uses get_view and "set view" commands to save and retrieve views of molecules. VMD handles things a little differently. It needs 4 parameters per view. They are [rotate/center/scale/global]_matrix. Saving and retrieving are done by processes below. They have been taken from view_change_render.tcl
I added the last part get_vp to see what the view point actually looks like (similar to get_view in PyMol).

proc save_vp {view_num} {
  global viewpoints
  if [info exists viewpoints($view_num)] {unset viewpoints($view_num)}
  # get the current matricies
  foreach mol [molinfo list] {
    set viewpoints($view_num,$mol,0) [molinfo $mol get rotate_matrix]
    set viewpoints($view_num,$mol,1) [molinfo $mol get center_matrix]
    set viewpoints($view_num,$mol,2) [molinfo $mol get scale_matrix]
    set viewpoints($view_num,$mol,3) [molinfo $mol get global_matrix]
  }
} 

 proc retrieve_vp {view_num} {
  global viewpoints
  foreach mol [molinfo list] {
    if [info exists viewpoints($view_num,$mol,0)] {
      molinfo $mol set rotate_matrix   $viewpoints($view_num,$mol,0)
      molinfo $mol set center_matrix   $viewpoints($view_num,$mol,1)
      molinfo $mol set scale_matrix   $viewpoints($view_num,$mol,2)
      molinfo $mol set global_matrix   $viewpoints($view_num,$mol,3)
    } else {
      puts "View $view_num was not saved"}
  }
}

 proc get_vp {view_num} {
  global viewpoints
  foreach mol [molinfo list] {
    if [info exists viewpoints($view_num,$mol,0)] {
puts "global viewpoints"
puts "molinfo $mol set rotate_matrix     \"$viewpoints($view_num,$mol,0)\""
puts "molinfo $mol set center_matrix     \"$viewpoints($view_num,$mol,1)\""
puts "molinfo $mol set scale_matrix      \"$viewpoints($view_num,$mol,2)\""
puts "molinfo $mol set global_matrix     \"$viewpoints($view_num,$mol,3)\""
    } else {
      puts "View $view_num was not saved"}
  }
}

Once you save this output (as a text file) it becomes very easy to retrieve the view even 

after closing down VMD. Simply copy and paste the output of get_vp.

Loops in Tcl/Tk

If loop

if { test1 } {
body1 
} elseif { test2 } {
body2 
} else {
bodyn
}

For loop

for { set i 1 } { $i <= 10 } { incr i } { # This will append all multiples of all numbers for the number n # into a variable called table set table "$table $i x $n = [expr $i \* $n]\n" }

 #Make the list
set list [list "Item 1" {Item 2} Last]
set text ""
foreach item $list {
 #Append the every item to $text
 set text "$text\n$item"
}

Foreach loop

#Make the list 
set list [list "Item 1" {Item 2} Last] 
set text "" 

foreach item $list { 
#Append the every item to $text 
 set text "$text\n$item"
 }

Find non-zero size files in Linux

http://www.unix.com/unix-dummies-questions-answers/25514-how-find-files-not-empty.html

You can use the find command to find non-zero length files:

find path -type f ! -size 0

In most UNIX implementations, the -size expression can also be used to search for file sizes of exactly N bytes (-size Nc), greater-than N bytes (-size +Nc), and less-than N bytes (-size -Nc).

The confusing thing is that the numeric following -size is, by default, 512-byte blocks not a byte count. The numeric must be followed by a 'c' for that. The following command will find files less than 2KB:

find . -type f -size -2048c -print

Conversely, for files greater than 2KB:

find . -type f -size +2048c -print

One note, the 512-byte block count does not directly translate into bytes. It's a long story. You can display a file's block usage with, under Solaris, ls -s. Stick to byte counts.

Oh, yes, one more thing, you can search for a specific range by using multiple -size expressions as long as they can all be satisfied by a single file. For example, locate files larger than 2KB and less than 8KB (inclusive):

find . -type f -size +2047c -size -8193c

Making VMD detect improperly named atoms

mol reanalyze does the job.
More here; http://www.ks.uiuc.edu/Research/vmd/mailing_list/vmd-l/19657.html

List of webpages that have come handy

http://ubuntugeeknerd.blogspot.com/2010/01/you-have-not-chosen-to-trust-thawte.html

For a better solution, try this:
sudo ln -s /etc/ssl/certs/Thawte_Premium_Server_CA.pem /usr/lib/ICAClient/keystore/cacerts/
Useful when you get this error;
" Citrix SSL Error 61: You have not chosen to trust “Thawte Server CA”".

NAMD sample scripts

http://www.ks.uiuc.edu/Research/namd/2.7/ug/node64.html

This file demonstrates the analysis of a DCD trajectory file using NAMD. The file pair.pdb contains the definition of pair interaction groups; NAMD will compute the interaction energy and force between these groups for each frame in the DCD file. It is assumed that coordinate frames were written every 1000 timesteps. See Sec. 12.1 for more about pair interaction calculations.

# initial config
coordinates     alanin.pdb
temperature     0

# output params
outputname      /tmp/alanin-analyze
binaryoutput    no

# integrator params
timestep        1.0

# force field params
structure       alanin.psf
parameters      alanin.params
exclude         scaled1-4
1-4scaling      1.0
switching       on
switchdist      8.0
cutoff          12.0
pairlistdist    13.5
stepspercycle   20
 
# Atoms in group 1 have a 1 in the B column; group 2 has a 2.
pairInteraction  on
pairInteractionFile pair.pdb
pairInteractionCol B
pairInteractionGroup1 1
pairInteractionGroup2 2

# First frame saved was frame 1000.
set ts 1000

coorfile open dcd /tmp/alanin.dcd

# Read all frames until nonzero is returned.
while { ![coorfile read] } {
  # Set firstTimestep so our energy output has the correct TS.
  firstTimestep $ts
  # Compute energies and forces, but don't try to move the atoms.
  run 0
  incr ts 1000
}
coorfile close

Drawing arrows in Gimp

Gimp (Gnu image manipulation program) is an open source image editor/manipulator that I use regularly to edit images and add text. I often need to add arrows to images and for some strange reason Gimp does not have an arrow button by default. But thanks to
http://registry.gimp.org/node/20269 there is now an easy way to accomplish that. For my purposes, I found these values ideal;

I could have learned a lot by simply reading the asterisk which says that negative values indicate fractions. For example a value of -2 indicates that the length of wings is 1/2 of the length of arrow. A value of -8 for brush thickness indicates that the brush is 1/8th of the length of arrow. For very long arrows one should select very low negative values (-100 or -200) to get reasonable looking arrows.

PLINK

PLINK is a free, open-source whole genome association analysis toolset, designed to perform a range of basic, large-scale analyses in a computationally efficient manner.

The focus of PLINK is purely on analysis of genotype/phenotype data, so there is no support for steps prior to this (e.g. study design and planning, generating genotype or CNV calls from raw data). Through integration with gPLINK and Haploview, there is some support for the subsequent visualization, annotation and storage of results.

PLINK (one syllable) is being developed by Shaun Purcell at the Center for Human Genetic Research (CHGR), Massachusetts General Hospital (MGH), and the Broad Institute of Harvard & MIT, with the support of others.


Basic usage for SNP lookup function
plink --lookup rs1475515
which outputs to the LOG file the following information

Bookmarks

My favorite online bookmarking tool is delicious.com. It is a pretty neat, no-nonsense website currently run by Chad Hurley and Steve Chen. I typically assign a tag to a saved link and looking at what other people have saved under that tag can be quite educational. Here's how to do it:
delicious.com/tag/name-of-tag

Solvating molecules with VMD

To perform MD simulations with NAMD with Charmm ff, I typically use Charmm to build the initial psf/pdb of the protein and then use VMD to solvate it. Now VMD does not like any psf file generated by Charmm. It likes to its own version of psf file. For this one has to first read in a pdb file and read out a psf file. For this the following script is helpful:


package require psfgen
topology ./top_all27_prot_na.inp
segment PROT { pdb xp5helixonXaxis.out.pdb }
coordpdb xp5helixonXaxis.out.pdb PROT
guesscoord
writepsf vmd.psf

package require solvate
solvate vmd.psf xp5helixonXaxis.out.pdb  -o solvate  -s WT -t 15

Useful resource:
http://www.ks.uiuc.edu/~villa/dna/

Fitting Protein Structures in VMD

# This script can be used to first load a reference structure and
# fit the second structure/trajectory to the reference.

mol new reference.pdb

mol new enzfibsol.psf
mol addfile md2.dcd type dcd first 0 last -1 step 1 waitfor all
mol addfile md3.dcd type dcd first 0 last -1 step 1 waitfor all
mol addfile md4.dcd type dcd first 0 last -1 step 1 waitfor all
mol addfile md5.dcd type dcd first 0 last -1 step 1 waitfor all

set backbone [atomselect 0 backbone]

for {set i 0} {$i < [molinfo 1 get numframes]} {incr i 1} {

set bk [atomselect 1 backbone frame $i]
set all [atomselect 1 all frame $i]
$all move [measure fit $bk $backbone]

}

Image manipulation in Matlab

a = double(imread('figures_from_microscope/500nm_particles_1200.tif'));
colormap('gray'), imagesc(a);

b = bpass(a,0,5,20);
#    bpass(img,lnoise, lobject, threshold)
colormap('gray'), image(b);

pk = pkfnd(b,25,3);
#     pkfnd(im, th, sz)

fid = fopen('exp.txt', 'w');
fprintf(fid, '%6.2f %12.8f\n', pk);
fclose(fid);

pk= pkfnd(b(200:400,200:400), 30,10);

# Below is a way to display multiple figures
figure(1); colormap('gray'), imagesc(a);
figure(2); colormap('gray'), image(b);
figure(3); colormap('gray'), image(b(200:400,200,400));
figure(4); colormap('gray'), image(a(200:400,200:400));

# This will put dots of 'r'ed on the particles of interest.
plot(pk(:,1),pk(:,2),'r.')
Looking at a small section of an image and checking for accuracy of the image manipulation is necessary. This is important since we make many assumptions on what particles to look at and what to ignore [this is what values in bpass and pkfnd do].

If var is a variable you are using in Matlab, it can be saved into a file using:
save('filename.txt', 'var', '-ascii')

These are the relevant links:
http://physics.georgetown.edu/matlab/tutorial.html
http://www.physics.emory.edu/~weeks/idl/tracking.html
http://www.physics.emory.edu/~weeks/idl/index.html

Getting header names in R and playing with csv files

names(myTable)

a<-read.csv("file.csv", header=T, sep=",")
names(a)
motility<-subset(a, Putative.orf.function=="Motility & Attachment")
write.csv(motolity, file="motility.csv")

Edit Feb 15th 2012

grep("Motility", a$Putative.orf.function)

gives the list of rows whose Putative.orf.function column contains the string "Motility".

Using this, one can get the entire row information.

df<-read.csv("file.csv", header=T, sep=",")

names(df)

 > names(df)
 [1] "Strain.Name"           "ID"                    "Location"             
 [4] "PA.ORF"                "Gene"                  "Gene.Abbrev."         
 [7] "Putative.orf.function" "Position.in.ORF"       "Frame"                
[10] "Genome.Position"       "Transposon.Direction"  "Transposon"           
[13] "F.Primer.Name"         "F.Primer.Seq"          "F.Primer.Position"    
[16] "F.Primer.Tm"           "Insert...F"            "R.Primer.Name"        
[19] "R.Primer.Seq"          "R.Primer.Position"     "R.Primer.Tm"          
[22] "R...Insert"            "WT.PCR.Len"            "Confirmed."           
[25] "How.confirmed."        "Notes"                
> df$Strain.Name[grep("Attachment",df$Putative.orf.function)]
  [1] PW10299 PW10300 PW1728  PW1729  PW1730  PW1731  PW1749  PW1750  PW1751 
 [10] PW1752  PW1753  PW1754  PW1755  PW1756  PW1757  PW2791  PW2792  PW2794 
 [19] PW2795  PW2941  PW2942  PW2943  PW2944  PW2945  PW2946  PW2947  PW2948 
...

[118] PW8681  PW8682  PW9465  PW9466  PW9467  PW9468  PW9469  PW9470  PW9471 
[127] PW9472  PW9473  PW9474 

Useful commands in Pymol

load file.pdb
create protein, all and segid PROT
create water, all and resname WAT
create res40, protein and resid 40
# To set display to orthoscopic
cmd.set('''orthoscopic''','''1''',quiet=0)
# To reset to default view
reset
hide everything, all
hide everything, protein
show sticks, protein
show surface, protein
set transparency=0.5
bg_color white
# To save images
ray 1200,1200
save file.png
distance (selection1), (selection2)
distance name = (selection1), (selection1) [,cutoff [,mode] ]
fit (selection), (target-selection)
# Colors everything by chain id
util.cbc

Resources:
http://pymol.sourceforge.net/newman/ref/S1000comref.html
http://www.pymolwiki.org/index.php/TOPTOC
http://www.pymolwiki.org
http://www.pymol.org/support
http://www.pymolwiki.org/images/7/77/PymolRef.pdf
http://people.mbi.ucla.edu/sawaya/tutorials/Graphics/pymol.html

The last link is particularly useful in the post-beginner stage of using Pymol.

load xxx.pdb
hide everything, all
create protein, all and not name h*
show sticks, protein
distance resi 22 and name CD and chain A, resi 22 and name CD and chain B
distance resi 22 and name CD and chain B, resi 22 and name CD and chain C
distance resi 22 and name CD and chain C, resi 22 and name CD and chain D
distance resi 22 and name CD and chain D, resi 22 and name CD and chain E
distance resi 28 and name NZ and chain D, resi 23 and name OD1 and chain C
distance resi 28 and name NZ and chain D, resi 23 and name OD2 and chain C
distance resi 28 and name NZ and chain C, resi 23 and name OD1 and chain B
distance resi 28 and name NZ and chain C, resi 23 and name OD2 and chain B
distance resi 28 and name NZ and chain B, resi 23 and name OD1 and chain A
distance resi 28 and name NZ and chain B, resi 23 and name OD2 and chain A
hide labels, all

March 2nd 2012

create protein, all and not resn DT+DG+NAG+HOH+OG6
create rna, all and resn DG+DT
create nag, all and resn NAG
create og6, all and resn OG6
hide everything, all
show cartoon, rna
show lines, rna
create closeprotein, (rna around 10) and not resn DG+DT
show sticks, closeprotein
distance /rna//D/DT`13/O4, /closeprotein//H/ARG`75/NH1
distance /rna//D/DT`13/O4, /closeprotein//H/ARG`75/NH2
distance /rna//D/DG`2/OP2, /closeprotein//H/ARG`77A/NH1
distance /rna//D/DG`2/OP2, /closeprotein//H/ARG`77A/NH2

March 5th, 2012

 # simple example: label residue 22's atoms with their names
label i. 22, name
 
# Label residue #44's alpha carbon with it's residue name, number and B-factor.
label n. CA and i. 44, "(%s, %s, %s") % (resn, resi, b)

This link has a lot of info on labeling. http://pymolwiki.org/index.php/Label

Printers

Here are a few commands to work with printers on Suse Linux.

# To check the status of all printers
lpstat -a
# To check the status of a printer
lpstat -p printer
# To set a default printer
lpoptions -d printer-name
lpq
# To print a file
lp -d printer file
# To remove a print job
lprm -Pprinter job-id

sudo -u pkc  /usr/sbin/cupsenable HP_LaserJet_6MP

To merge pdf files in Mac

Here's a neat little trick to combine pdf files on a Mac
http://macintoshhowto.com/leopard/how-to-merge-pdf-files-with-preview-in-leopard.html

VMD script to pick a box within a box

set cellxmin [lindex [lindex [measure minmax $all] 0] 0]
set cellxmax [lindex [lindex [measure minmax $all] 1] 0]
set cellymin [lindex [lindex [measure minmax $all] 0] 1]
set cellymax [lindex [lindex [measure minmax $all] 1] 1]
set cellzmin [lindex [lindex [measure minmax $all] 0] 2]
set cellzmax [lindex [lindex [measure minmax $all] 1] 2]

for {set i 0} {$i < 100} {incr i} {
# This lets you pick a random box of size 10x10x10 in the simulation box
set xmin [expr $cellxmin + ( ($cellxmax-$cellxmin-$boxsize) * rand())]
set ymin [expr $cellymin + ( ($cellymax-$cellymin-$boxsize) * rand())]
set zmin [expr $cellzmin + ( ($cellzmax-$cellzmin-$boxsize) * rand())]
set xmax [expr $xmin + $boxsize]
set ymax [expr $ymin + $boxsize]
set zmax [expr $zmin + $boxsize]

set selection [atomselect top "(x>$xmin and x<$xmax and y>$ymin and y<$ymax and z>$zmin and z<$zmax)"]
# This gives the number of atom in that selection
set n [$selection num];
#This gives the number of non-bonded interactions in the 10x10x10 box selected.
puts $file [expr $n * ($n -1)/2 ]
}

Useful perl code

#!/usr/bin/perl

$infile  = "./test.pdb";
$outa = "./test.pdb.A";
$outb = "./test.pdb.B";
open (FH,$infile) || die;
open (OUTA, ">$outa") || die;
open (OUTB, ">$outb") || die;

while ($buf = )
{
    $x = substr($buf,16,1);
    if ($x eq 'A') { print OUTA $buf; }
    elsif ($x eq 'B') { print OUTB $buf; }
    else { print OUTA $buf; print OUTB $buf; }
}

close FH;
close OUTA;
close OUTB;

# This code reads in a file and if the 17th column in each line matches A/B, that line is put in OUTA/OUTB.This is useful in separating multiple coordinates for particular residues in 3D coordinate file of a protein. Sometimes studying all possible conformations is important especially if it involves an aromatic ring flip or flips in residues with amine ends.

Delete lines matching a string on VIM

:g/pattern/d

:v/pattern/d

Go to col on vim/vi

You can quickly jump to a specified column in Vim by simply typing the column number followed by a pipe.

80| # go to the 80th column

Thanks to DailyVim

Min-Max on R plots

d <- density(rnorm(1000))
plot(d)
d$x[which.max(d$y)]
abline(v=d$x[which.max(d$y)])

http://r.789695.n4.nabble.com/Find-x-value-of-density-plots-td4165926.html

Great source for 'sed' examples

http://sed.sourceforge.net/sed1line.txt

Bash commands in R

x<-as.matrix(read.table(pipe("awk -f cut.awk F120x1")))

NAMD inputs from VMD

set all [atomselect top all]
measure minmax $all
set cellxmin [lindex [lindex [measure minmax $all] 0] 0]
set cellxmax [lindex [lindex [measure minmax $all] 1] 0]
set cellymin [lindex [lindex [measure minmax $all] 0] 1]
set cellymax [lindex [lindex [measure minmax $all] 1] 1]
set cellzmin [lindex [lindex [measure minmax $all] 0] 2]
set cellzmax [lindex [lindex [measure minmax $all] 1] 2]

puts "cellOrigin              [format "%7.3f %7.3f %7.3f" [lindex [measure center $all] 0] [lindex [measure center $all] 1] [lindex [measure center $all] 2]]"
puts "cellBasisVector1        [format "%7.3f %7.3f %7.3f"      [vecsub $cellxmax $cellxmin] 0 0]"
puts "cellBasisVector2        [format "%7.3f %7.3f %7.3f"  0   [vecsub $cellymax $cellymin] 0]"
puts "cellBasisVector3        [format "%7.3f %7.3f %7.3f"  0 0 [vecsub $cellzmax $cellzmin] ]"



Pretty primitive, but it works.

Common VMD Tcl/Tk commands

mol modselect 0 0 resname BGLC and noh
mol modselect 0 1 resname BGLC and noh
mol modselect 0 2 resname BGLC and noh
mol modselect 0 3 resname BGLC and noh

mol smoothrep 0 0 5
mol smoothrep 1 0 5
mol smoothrep 2 0 5
mol smoothrep 3 0 5


mol modstyle 0 0 Licorice 0.300000 10.000000 10.000000

mol addrep 0
mol modselect 1 0 protein
mol modstyle  1 0 Licorice 0.300000 10.000000 10.000000
mol modstyle  1 0 vdw
mol smoothrep 1 1 5

mol off 0
mol off 1
mol off 2
mol off 3

mol showrep 0 0 0
mol showrep 0 0 1

color Display Background white

Calculating phi/psi values in VMD

First source the following script:
--- start ---
proc all_dihed_angle { a1 a2 a3 a4 } {
  # Delete all existing Dihdral labels so that the one we add has index 0.
  label delete Dihedrals

  # Use the top molecule
  set molid [molinfo top]

  # Add the dihedral monitor
  label add Dihedrals $molid/$a1 $molid/$a2 $molid/$a3 $molid/$a4

  # Get all values
  return [label graph Dihedrals 0]
}
--- end ---

Then run this script:

mol new enzfibsol.psf
mol addfile md7.dcd first 0 last 10 waitfor all
set phifile [open "cellulose_dihed_phi.tcl" w]
set psifile [open "cellulose_dihed_psi.tcl" w]

set oxygen [atomselect top "resname BGLC and resid 436 and name O1"]
set ox [$oxygen get index]
set h1 [expr {$ox-1}]
set c1 [expr {$ox-2}]
set c4 [expr {$ox-9}]
set h4 [expr {$ox-8}]

puts $phifile [all_dihed_angle $h1 $c1 $ox $c4]
puts $psifile [all_dihed_angle        $c1 $ox $c4 $h4]

close $phifile
close $psifile
quit

Resize images on a console

sips --resampleWidth 2400 IMG_0062.JPG --out test.jpg
(double dash) resampleWidth and (double dash) out

This Mac command does a fantastic job of resizing large images in to more manageable file sizes.

Edit: A width of 2400 is too large. Use 1600 or lower instead.

Wildcards in Tcl/Tk

http://phaseit.net/claird/comp.lang.tcl/fmm.html

Very useful link on wildcards in Tcl Tk. This line solved my problem;
mol new [glob fitted*.pdb rev*.pdb]

A little background would help greatly. I had a set of directories each of which had either a fittedxxx.pdb or revxxx.pdb file. I wanted a single tcl script that would go into these directories and read the fitted/rev file and perform calculations. Since the name of the file in each directory was different, I needed a way to select whatever was available. glob does exactly that. 

Install packages in R

R CMD INSTALL mypkg -l /my/own/R-packages/

install.packages("mypkg", lib="/my/own/R-packages/")

ssh-keygen

To reset IP and known_host keys, use the following command:

ssh-keygen -R moldyn.ornl.gov -f /home/pkc/.ssh/known_hosts
RSA host key for moldyn.ornl.gov has changed and you have requested strict checking.

That #$%@ beep

1. xset b off
2. set b off
3. Turn off beep on Gmixer (or whatever audio controller you run)

These commands turn off the annoying beep that accompanies [tab], [auto-fill], and errors in the linux environment.

png figures in a loop in R

This sample script creates 100 png figures with the iteration number included in the name.

for(i in 1:100) {

a<-rnorm(100)

name<-paste("test",i,".png",sep="")

png(name)

plot(a/max(a),type="l")

dev.off()

}

Running multiple serial jobs on clusters

This script describes a way to run multiple serial jobs on clusters.

#PBS -N serial.analysis
#PBS -j oe
#PBS -l walltime=4:00:00,nodes=22:ppn=16,gres=widow2

cd $PBS_O_WORKDIR
date
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$PBS_O_WORKDIR

sort $PBS_NODEFILE | uniq > ./TEMP.NODEFILE

n=`wc -l ./TEMP.NODEFILE |awk '{print $1}'`
for((i=1;i<=$n;i++)); do j=`expr $i - 1`; sed -n ''$i','$i'p' ./TEMP.NODEFILE > host.$j; done

$MPI -np 1 --hostfile host.0 $NAMD ./w1_md3_resid_100.namd > ./w1_md3_resid_100.namd.out &
$MPI -np 1 --hostfile host.0 $NAMD ./w1_md3_resid_101.namd > ./w1_md3_resid_101.namd.out &
....
$MPI -np 1 --hostfile host.0 $NAMD ./w1_md3_resid_174.namd > ./w1_md3_resid_174.namd.out &
$MPI -np 1 --hostfile host.1 $NAMD ./w1_md3_resid_175.namd > ./w1_md3_resid_175.namd.out &
....
$MPI -np 1 --hostfile host.1 $NAMD ./w1_md3_resid_246.namd > ./w1_md3_resid_246.namd.out &
$MPI -np 1 --hostfile host.2 $NAMD ./w1_md3_resid_247.namd > ./w1_md3_resid_247.namd.out &
$MPI -np 1 --hostfile host.2 $NAMD ./w1_md3_resid_251.namd > ./w1_md3_resid_251.namd.out &
..

and so on.

Conditional plots in R

This link says it all.

 Al - 19 Jan 2007 15:33 GMT

Hi!
Given a 1000x2 matrix x, i.e. 1000 2d-points. How can I directly (!)
plot/visualize only those which fulfill a certain condition, say
x[,1]>0 & x[,2]>0?
Thanks,
Al

Marc Schwartz - 19 Jan 2007 21:24 GMT

> Hi!
> Given a 1000x2 matrix x, i.e. 1000 2d-points. How can I directly (!)
> plot/visualize only those which fulfill a certain condition, say
> x[,1]>0 & x[,2]>0?
> Thanks,
> Al

First, R specific programming queries are best posted to the r-help
list. More information at:

http://www.r-project.org/

under the Mailing Lists link.

Second, the easiest way to do this is to subset the data while using
plot.formula(), which is a plot method that will be dispatched if the
first argument to plot() is a formula, such as 'y ~ x'. See ?plot and
?plot.formula for more information.

As an example, let's create a data frame with 2 columns and 1000 rows,
each column comprised of 1000 random normal deviates:

DF <- data.frame(x = rnorm(1000), y = rnorm(1000))

> str(DF)

'data.frame': 1000 obs. of 2 variables:
$ x: num -0.0265 0.1554 -0.1050 -0.9697 -0.3430 ...
$ y: num 1.386 -1.356 -1.170 0.426 0.204 ...

Now, let's plot the data meeting the criteria you indicated above:

plot(y ~ x, data = DF, subset = (x > 0) & (y > 0))

This says to generate a scatter plot, using 'DF' as the input data set,
and specifying the logical subset where x and y are both > 0. Modify
the subset argument as you require.

Note that in this usage, the 'data' argument must be a data frame. So
you could coerce your matrix (if needed) by using:

DF <- as.data.frame(matrix)

HTH,

Marc Schwartz

And this is the way (ex.) to get the length of the array that satisfies the conditions.

sum((x > 0) & (y > 0), na.rm=TRUE)

Custom replacements in VI

Thanks to Steve for this pearl of VI wisdom.
For search and replace operation in VI, parenthesis remembers the item. This can be used for custom replacements later. Try the command below on a file with this data:

29 5 2010
30 6 2011

:1,$s/[0-9]* \([0-9]*\) \([0-9]*\)/\2 \1/
VI can remember 9 such fields.

Potential Energy Scans with NWChem

While performing PES on dihedrals of relatively large molecules using NWchem, one can get stuck during the minimization process. Possible ways to overcome those problems are:
1. Try a smaller step size
2. Use the "tight" function in the driver
3. Redo everything with a smaller molecule that can mimic the dihedral in question

Useful "vi" tips

1. If you want to merge rows ( for exp. row 1 to 20 ) do:
:1,20 j
2. To replace string with a new line
:%s/ TI/^MTI
To get ^M use following
Hold down the Ctrl, and press v m
DO NOT USE Shift+6 and Shift+M
3. To move a string to a new line w/o hitting Enter, do this;
:s/string/\r&/g

Flush VMD

While running VMD analysis scripts over large number of trajectories, it is a good idea to free up memory after every step like so:
for {set j 1} {$j < 101} {incr j}
{
mol addfile /.$j.dcd first 0 last -1 waitfor all
set ca [atomselect top "name CA"]
set cb [atomselect top "name CB"]
....
$ca delete
$cb delete
animate delete all
}
The last three lines will make life much better.

Using Bit Bucket

Following my friend's recommendation I decided to check out this site called Bit Bucket which is a free code sharing/hosting website where one can essentially store their code and any changes they make to that file. Following are the steps to get things up and running:
1. Install Mercurial (if you don't already have it)
2. Open an account with https://bitbucket.org
3. Create a repository on bitbucket.org. Let's call it bit_repo.
4. Create a /home/user/.hgrc file with the following lines on your desktop.
[ui]
; editor used to enter commit logs, etc. Most text editors will work.
editor = vi
username = firstname lastname
5. cd in to a directory on your desktop which contains codes that you want to backup.
5.1. Type> hg clone https://bitbucket_user_id@bitbucket.org/bitbucket_user_id/bit_repo
destination directory: bit_repo
http authorization required
...
Then cd bit_repo
6. To add a specific file to Bitbucket, type
hg add file1.c
hg commit -m "I commit to add this file to bitbucket" [The msg in quotes can by anything]
hg push https://bitbucket_user_id@bitbucket.org/bitbucket_user_id/bit_repo
7. To remove a file type
hg remove file2.c
hg commit -m "I commit to remove this file"
hg push https://bitbucket_user_id@bitbucket.org/bitbucket_user_id/bit_repo
8. To copy files/repositories from bitbucket, type
hg clone https://bitbucket_user_id@bitbucket.org/bitbucket_user_id/bit_repo
9. Here is a Bitbucket 101 http://confluence.atlassian.com/display/BITBUCKET/Bitbucket+101
10. To keep things safe I would just softlink files to their sources on my desktop.