Showing posts with label science. Show all posts
Showing posts with label science. Show all posts

Friday, 28 February 2014

Intro to GitHub (for scientists)

The problem

I talked to a young woman yesterday who is a bio-engineering postdoc at Stanford. She has some code that she'd like to 'upload' to GitHub. She admitted that, well, actually, she hadn't managed to get any of her code onto her GitHub account yet, and she looked so overwhelmed and dejected that I felt bad. I know that frustration.
So here is yet another blog post to try and help with the learning curve that is git. I'm writing this for the scientist who has written some code and wants to share it on GitHub. Most scientists write code in order to accomplish a particular task, and thus are not familiar with professional programming practices, including writing documentation, unit testing, SCM (source control management) and version control. There are great benefits to learning these techniques, and their usefulness is becoming more and more apparent to academics as research relies more and more on computer programs. In fact, the table of contents for the journal Nature Methods just arrived in my inbox with a lead editorial on reproducible research. I quote:

Nature Methods strongly encourages researchers to take advantage of the opportunity that code repositories, such as GitHub, provide to improve a software tool before submission. Even if others do not examine and test the code, the act of preparing the code and necessary documentation for deposit and use by others will help avoid delays in publication.
In short, if you are coding and publishing work derived from your code, the process of uploading your code for collaboration will help bring it up to a publishable quality. And yes, it's even more important than making the figures pretty.

GitHub == collaboration

Firstly, if you just want to upload some code, you need to take step back, a deep breath and 'think different'. There is no 'upload' button on GitHub for a good reason. It is not built for uploading code and leaving it to rot on a server, but for fostering collaboration between programmers. Thus, in order to share your code on GitHub, you first need to get it ready for collaboration. To do that, you need to set up version control. This is somewhat more complicated than finding the 'track changes' setting in Microsoft Word, but it is also far more useful.
As you try to think different, keep in mind that GitHub was built with a particular set of workflows in mind. Those workflows have to do with managing a code-base that is constantly being updated by multiple contributors. GitHub tries to ease the difficulties of this type of collaboration by bringing several things together.
  • version control -- similar to 'track changes' in a word processor, but for whole projects
  • cloud-based storage -- simultaneously backup and share your code
  • user accounts -- keep it private, share with a team or make it public as your needs change
  • social -- there's messaging, so you can talk, argue, document, discuss.... collaborate!
But version control is the main thing that sets GitHub apart from any number of social sharing platforms, and makes it so powerful for people who code. To use GitHub effectively you therefore need to understand the basics of version control with git. It's hairy and scary at first, but no worse than... well, ok, it is worse than a lot of things, but sometimes a learning curve is the price you pay for really capable tools. So, put a set of bookmarks in your browser, make a cheatsheet and keep it handy. A local cheatsheet with a good searchable title does wonders. If you're thinking 'Oh, I'll just use it once' or 'I'll remember', well, git is for professionals. Are you a professional?

Git == version control

Version control is to 'track changes' for a text document, as Superman is to Tarzan; as Microsoft Word is to TextEdit, as New York City is to Detroit. Firstly, with version control, you are tracking changes to an entire project, not a single document. Changes are tracked line by line with comments and attribution through time as the project grows, changes, and splits; as subroutines develop into full projects of their own; as new owners take control of the code base. It is flexible and thorough and reliable. You should learn to use it.
Use case: You've got some code you developed for your thesis that you want to upload to GitHub. Maybe someone else will find it useful. There are GUI front ends to git, which may help with many tasks, but git was designed to be run from the command line, and this simple use case it not too difficult to master at that level, so let's just go for it.

Download and install git

Walk through the steps at github set-up. Today, they suggest that you download their native app, but note that this only manages part of the workflow. The steps you have to do to get git working still involve the terminal. Be brave. The steps are:
  1. get a GitHub account
  2. download and install git
  3. link your local system to your GitHub account
    • tell git your name, email, GitHub login information
    • set up security keys so that GitHub knows you are you

Prepare your codebase

There are a couple of adjustments you probably want to make before releasing your code in the wild. GitHub recommends that all code comes with a License, a Readme and a .gitignore file.
  • License If your code comes with a license, it's easier for collaborators to re-use the code and build on it. Specifically, it makes clear what they are allowed to do. It may not be important to you, but your code will be easier to share if you make it clear what the rules are.
  • Readme
    You have to explain your work at some point. If you do it in a file named README, GitHub will automatically put it on the front page of the repository. This is very helpful for anyone trying to understand what you did and why. The README can be just a text file. If it is in markdown, perhaps with additional flavoring GitHub will render it with headings and styles, which is much nicer for the reader and not difficult for the writer.
  • .gitignore
    Some files don't need to be tracked. For instance, some old Mac directories contain .DS_Store files with directory display information for the Finder app. That doesn't need to be part of the repository. So here's my .gitignore for an old MatLab project:
      $cat .gitignore
      .DS_Store 
    
    Pretty simple. You might also add *.log or tmp/ to the .gitignore file, depending on your context. Basically, any files that are automatically generated or updated on compile should not be tracked.

Make a local repository

The project you want to get onto GitHub probably consists of a directory or directory tree containing a series of text files, and possibly some image files or data files. In order for git to track changes to this project, you have to put these files into a repository. This is simple to do once you've prepared it for sharing.
Find your command line. On a Mac, you can use the Terminal app. Navigate to the directory holding your project. If you have never ever used the command line before, this might be challenging. If you want to dive in, you can certainly do so with three little commands: ls, cd, pwd. You can look at the man pages for these commands by typing, for example $man ls, or you could try a crash course in using the command line.
Once you get to your project directory, type:
$git init   
You should see a message from git, something like:
Initialized empty git repository in /Users/suz/programming/octave/OrX/.git/
This initialises an empty repository, which looks like a file named .git. You can check that it's there by typing
$ls -a
Git gives some feedback about what it has done, but I often find it useful to check with
$git status
after each command to see what has happened.
Now we can get the project into the repository. To do this, first type
$git add .   
This prepares git to add your files to the repository, a process known as 'staging'.
Note: The . tells git to stage all the files in the directory tree to the repository. This is very handy when we add files because they don't all have to be specified by name. On the other hand, it isn't ideal, because there are often binary files or .log or even image files that update automatically. You won't want to keep track of those changes. Fortunately, git will automatically look for the .gitignore file we already prepared to get the list of exceptions.

Ready for some commitment? Type:
$git commit -m 'initial commit' 
You should get a full list of the files that git is committing to the repository. Check the status when it's done and you should give a reassuring message:
# On branch master
nothing to commit (working directory clean)

Success! now your project is actually in the repository and git can track any changes to the files. The repository will also keep all the messages that you put with each commit. Always use the -m flag and take the time to add a meaningful message.

Congratulations! Your code is now in a git repository, under version control. You are ready to collaborate.

Share your work


Make a repository on GitHub

  1. Log into your GitHub account
  2. Create a new repository on GitHub
    On your profile page, in the 'Repositories' tab is a list of repositories that you've contributed to. At the upper right should be a bright green 'New' button.
  3. Follow the directions in the form, adding an informative description so others know what treasure they have found.
Congratulations! You have a GitHub repository to share your code from!

Link the repositories

In git terminology, the current state of your code is the head of a branch. By default, the initial commit is called the 'master' branch. You can make other local branches, and probably should to try out new features. You can also make remote branches. At this point, your new GitHub repository is essentially an empty remote branch. By custom, this branch will be referred to as 'origin'. To point git to it, type (on one line), with your appropriate username and project title:

$git remote add origin https://github.com/username/project.git

This command translates roughly as "Dear git; Please be so kind as to add a connection to a remote repository. I will be referring to the remote repository as 'origin' in our future correspondence. It can be found at https://...... Thank you for your kind assistance in this matter. Sincerely, yours truly, esq."

Upload your code

Ok, ok, there is no upload on GitHub, but it is payoff time. Once you have a local repository linked to a remote repository, you can just push the code from one to the other.

$ git push -u origin master

Translation: "Dear git; Please push the recent changes I committed to my local repository, known as master, into the remote repository known as origin. Also, please be aware that I'd like this to be the default remote repository, sometimes referred to as 'upstream'. Thank you again for your kind assistance. I am forever in your debt. Sincerely, Yours truly, esq. and, etc."

Success!!! You have now successfully pushed your code to GitHub.

Or at least I hope you were successful. If not, if you've tried to follow this post and the directions at GitHub and you still feel lost, there is more help out there. Many universities are running Software Carpentry bootcamps to help students and faculty develop more professional programming skills. The skills taught aim to improve software collaboration and impart the skills needed to carry out reproducible research. Two key tools they teach are version control with git and collaboration via GitHub.
Live long and collaborate!

Wednesday, 29 January 2014

UV reactive bead spectra

UV reactive bead spectra

Background

I was interested in designing some educational experiments using UV-reactive beads to teach about the presence/absence and intensity of UV light under different atmospheric conditions. The beads turn from clear or colorless to various bright colors when exposed to UV light. Unfortunately, these beads are just too sensitive – they even react strongly to the stray 400 nm light coming through glass windows. And reach full color intensity within a minute, even at 51° N in February! It's fun to watch, but not very useful for detecting physiological UV conditions for vitamin D synthesis.

I don't have a lot of scientific information abou these, so the range of UV light needed for the color change is not well characterized, and the resultant absorbance spectra of the beads (related to their aparent colors) is not readily available. I decided to study them with a reflectance spectrophotometer to see what I could learn.

It may not be useful for the educational project I had in mind, but I'll write it up here anyway, along with the R-code used to analyze the data….

Experimental

I ordered UV Reactive beads from UV gear, UK.

I used a table in the garden on a sunny day for the experiments. The sky was open directly above and direct sunlight came from the south. There were a few leafless tree branches in the way, and some buildings nearby, so it was not full sky exposure, but at least there was full direct sunlight. The table was covered with a small white blanket to provide a consistent background. The beads equibilbrated 10-15 min in direct sunlight before spectral measurements were made.

I used an Ocean Optics model USB2000+ spectrometer connected to a 1 meter UV-Vis fiber optic cable (part number QP400-1-UV-VIS from Ocean Optics) with a cosine correction filter at the end (Ocean Optics part CC-3-UV-S) to minimize stray light. The spectrometer was connected by USB to a MacBook Pro (OS X mountain lion 10.8.1) running SpectraSuite software for data collection. Spectra were measured with 30 averages in reflectance mode and saved as tab-delimited text files. The reflectance measurements required a 'reference spectrum' (light spectrum) of the white blanket and a second, 'dark' reference, which was obtained by blocking the end of the cosine correction filter. To collect the reflectance spectra, I pointed the fiber optic light pipe at individual beads from a distance of < 0.5 cm.

Spectral analysis can be done in the SpectraSuite software, however, many basic functions had not (yet?) been implemented in the OS X version of the software. Also, I wanted to practice my R skills, so I decided to load it into R and see what I could learn. R has multiple spectroscopy packages, and in fact a new one has come out since I started this project. I looked into using hyperspec, but the data structure seemed overly complicated for the simple analysis I had in mind. So what follows is my simple R analysis of reflectance data on sunlight-exposed UV beads.

Load in the Data

# change working directory
try(setwd("/Users/suz/Documents/vitD Schools/UV bead exp 19022013/spectrometer data/"))

# data is tab-separated, with a header. The end of file wasn't recognized
# by R automatically, but a read function can specify the number of rows
# to read.
spec.read <- function(spec.name) {
    read.table(spec.name, sep = "\t", skip = 17, nrows = 2048)
}

# the function can then read in all the files ending in '.txt', put them
# in a list.
spec.files <- list.files(pattern = ".txt")
df.list <- lapply(spec.files, spec.read)

# convert the list to a matrix: the first column is the wavelengths and
# the other columns are experimental data -- one reflectance measurement
# at each wavelength.
spec.mat <- matrix(df.list[[1]][, 1], nrow = 2048, ncol = (1 + length(spec.files)))
spec.mat[, 2:11] <- sapply(df.list, "[[", 2)

matplot(spec.mat[, 1], spec.mat[, 2:11], type = "l", lty = 1, lwd = 3, xlab = "wavelength (nm)", 
    ylab = "reflectance", main = "UV color changing beads, after UV exposure")
text(500, 25000, "Raw Data")

Looks like we need to do some clean-up!

Clean up the data

Baselines and edges

At the edges of the spectral range, the reflectance data is dominated by noise, so it isn't useful. The baselines for the different spectra also need aligning, and we'll scale them to the same intensity for comparison. The intensity range observed in the data depends somewhat on the angle of the probe, even with a cos filter in place.

# define terms for the processing:
nm.min <- 400
nm.max <- 800  # edges of the displayed spectrum
base.min <- 720
base.max <- 900  # define baseline correction range
peak.range.min <- 420
peak.range.max <- 680  # where to find peaks for scaling.

# remove ends of the data range that consist of noise
spec.mat <- spec.mat[(spec.mat[, 1] > nm.min) & (spec.mat[, 1] < nm.max), ]

# normalize baselines, set baseline range = 0.
spec.base <- colMeans(spec.mat[(spec.mat[, 1] > base.min) & (spec.mat[, 1] < 
    base.max), ])
spec.base[1] <- 0  # don't shift the wavelengths
spec.mat <- scale(spec.mat, center = spec.base, scale = FALSE)

Choose colors for the plot by relating the file names to R's built-in color names.

bead.col <- sapply(strsplit(spec.files, " bead"), "[[", 1)
# replace un-recognized colors with r-recognized colors (see 'colors()')
bead.col <- gsub("darker pink", "magenta", bead.col)
bead.col <- gsub("dk ", "dark", bead.col)
bead.col <- gsub("lt ", "light", bead.col)
bead.col <- gsub("lighter ", "light", bead.col)

# plot corrected data
matplot(spec.mat[, 1], spec.mat[, 2:11], type = "l", lty = 1, lwd = 3, col = bead.col, 
    xlab = "wavelength (nm)", ylab = "reflectance", main = "UV color changing beads, after UV exposure")
text(500, 10, "Baseline Corrected")

From this plot, we can see that the lighter colored beads have smaller peaks than the darker beads. The lighter color probably represents less dye in the beads. There seems to be a lot of variation in the peak intensity of some beads, particularly the yellow beads and the dark blue beads. Based on the width of the peaks, the purple and magenta beads appear to have mixtures of dyes for both pink and blue colors. The dark blue beads appear to be either a mixture of all the dye colors or a mixture of pink and blue, but much more dye is used than for the paler pink or blue beads. The yellow bead spectra is oddly shaped on the short wavelength side, probably due to instrument cutoffs around 400 nm.

Scaling and smoothing

Now I'll scale the data to the same range. It turns out that the R command 'scale' is perfect for this.

# scale the peaks based on the min reflected intensity
peak.range <- which((spec.mat[, 1] > peak.range.min) & (spec.mat[, 1] < peak.range.max))
spec.min <- apply(spec.mat[peak.range, ], 2, min)
spec.min[1] <- 1  # don't scale the wavelengths
spec.mat <- scale(spec.mat, center = FALSE, scale = abs(spec.min))

The spectra are also jittery due to noise. This can be removed by filtering. This filters over a range of 10 points.

data.mat <- spec.mat[, 2:11]
dataf.mat <- apply(data.mat, 2, filter, rep(1, 10))
dataf.mat <- dataf.mat/10
specf.mat <- matrix(c(spec.mat[, 1], dataf.mat), nrow = dim(spec.mat)[1], ncol = dim(spec.mat)[2], 
    byrow = FALSE)
matplot(specf.mat[, 1], specf.mat[, 2:11], type = "l", lty = 1, lwd = 3, col = bead.col, 
    xlab = "wavelength (nm)", ylab = "reflectance", main = "UV color changing beads, after UV exposure")
text(500, 0.2, "Scaled and Smoothed")

I could attempt to prettify the spectra further by using actual colors from pictures of the beads. It looks like the Bioconductor project has a package 'EBImage' that should be just what I want, but it looks like I need to update to R version 3.0.1 in order to run it. So I guess I get some spot colors from JImage.

bead.yellow <- rgb(185, 155, 85, maxColorValue = 255)
bead.orange <- rgb(208, 155, 82, maxColorValue = 255)
bead.purple <- rgb(110, 25, 125, maxColorValue = 255)
bead.pink <- rgb(190, 100, 120, maxColorValue = 255)
bead.dkblu <- rgb(22, 18, 120, maxColorValue = 255)
bead.ltblu <- rgb(130, 135, 157, maxColorValue = 255)
bead.dkpink <- rgb(200, 25, 140, maxColorValue = 255)

bead.col2 <- c(bead.dkpink, bead.dkblu, bead.dkblu, bead.pink, bead.ltblu, bead.orange, 
    bead.pink, bead.purple, bead.yellow, bead.yellow)

matplot(specf.mat[, 1], specf.mat[, 2:11], type = "l", lty = 1, lwd = 3, col = bead.col2, 
    xlab = "wavelength (nm)", ylab = "reflectance", main = "UV color changing beads, after UV exposure")
text(500, 0.2, "Colors from Photo")

Analysis

I can carry this anlaysis further by quantifying the wavelength of the peak absorbance and the peak width (usually Full Width at Half Maximum – FWHM) for each spectrum. This could be useful in further analyses, reports, or as a feature in a machine learning approach.

To extract this information, I could try to fit a series of gaussians to the peak, representing the fraction of pink, blue or yellow dye present, but the quality of the data doesn't really justify this, particularly as I don't have a good shape for the yellow absorbance peak or adequate reference data for each of the dyes. As a quick and dirty method, I'll take the median position of the data values that are at 95% of peak. That should give something in the center of the peak. Since the peaks have all been scaled, that corresponds to the center of the data values < -0.95.

Likewise, the usual peak width (FWHM) woud be the range of values < -0.5, however, the poor baseline at shorter wavelengths makes this range more or less unusable. For this reason, we can take a more well-behaved approaximate peak width measurement as the range of reflectance values < -0.8.

# approximate peak wavelength
indcs <- apply(specf.mat[, 2:11], 2, function(x) median(which(x <= (-0.95))))
pks <- specf.mat[indcs, 1]

# approximate peak width
lowidx <- apply(specf.mat[, 2:11], 2, function(x) min(which(x <= (-0.8))))
highidx <- apply(specf.mat[, 2:11], 2, function(x) max(which(x <= (-0.8))))

pkwidth80 <- (specf.mat[240, 1] - specf.mat[140, 1])/100 * (highidx - lowidx)

features <- data.frame(pks, pkwidth80, bead.col, bead.col2)
features <- features[order(pks), ]
features
##      pks pkwidth80  bead.col bead.col2
## 10 436.4     75.21    yellow   #B99B55
## 6  449.8     88.86    orange   #D09B52
## 9  451.7     98.07    yellow   #B99B55
## 7  529.3    111.35      pink   #BE6478
## 4  533.6     99.18 lightpink   #BE6478
## 1  549.8    110.61   magenta   #C8198C
## 8  568.7    133.10    purple   #6E197D
## 3  581.1    179.93  darkblue   #161278
## 2  585.4    136.79  darkblue   #161278
## 5  603.0     76.32 lightblue   #82879D

# save your work!
write.csv(features, "UVbead data features.csv", quote = TRUE)

One side effect of this is that we now have a small table of metrics that describe the relatively large original data set reasonably accurately and could be used to classify new data. In current data analytics parlance, this is known as 'feature extraction'. In traditional science, these characterizing features could be combined with others from different studies (crystallography, electrochemistry, UV-vis, IR, Raman, NMR, …) to help predict the effects of a chemical change or a different chemical environment on the behaviour of the molecule. Such studies are traditionally used to help direct synthetic chemists toward better products.

Conclusions

From this analysis, we can see that the blue beads are absorbing at longer wavelengths than the yellow, and the pink and purple beads have absorbances in between. The darker colors have broader absorbance peaks than the lighter colors, with the darkest blue having a range that appears to cover the purple, yellow and blue regions. These darker beads probably contain combinations of the dyes used for the different colors, and not just more of the dye used for the light blue beads.

The reflectance spectra are not as observant as our eyes. For instance, the yellow and orange beads are readily distinguished by eye, but not so clearly in the observed spectra or the extracted features. This may be due to the 400nm cutoff of the light pipe, which distorts the peak shapes for the yellow and orange beads. Ideally, we could observe the changes in the absorption spectra over the whole range 280-800 nm as a function of time. The current reflectance spectrometer setup, however, is only capable of capturing the 400-750 nm range. This means that we do not have access to the interesting behavior of these dyes at shorter wavelengths.

Chemically, I expect that absorption of the UV light in the 300-360 nm range causes a reversible conformation change in the dye molecules, as has long been known for the azo-benzenes and their derivatives. After the conformation change, the absorption maximum of the dye is shifted to a much longer wavelength. If the UV dyes in these beads are very closely related to each other, is likely that the yellow beads absorb at relatively short UV wavelengths and the blue beads at longer wavelengths before the transition, however, without measuring the UV absorption spectra we cannot know this. It is quite possible that their UV absorption spectra are nearly indistinguishable and the effects of their chemical differences are only apparent in the visible spectra. Without access to better hardware or more information about the molecules involved there is no way to know.

Next steps

We do have access to one thing that varies in the appropriate way. The experiments shown here were taken at relatively low UV light levels (February in England). During the summer, the angle of variaton of the sun is much larger. At dawn, it is at the horizon, while at noon, it is about 60° higher. Since light scattering in the atmosphere is a strong function of wavelength, the relative intensities of light at 300, 360, and 400 nm will vary with the angle of the sun. If the beads have different absorption peaks in the UV range, the time dependence of their color changes should vary by the time of day.

The simplest way to measure this is not with a spectrometer, but probably by following color changes in video images.