I only learned about Gibson when I started my postdoc, and it completely changed how I approached science. In some experiments with Ethan when I was in the lab, I was blown away when we realized that you don’t even need Gibson mix to piece a plasmid back together; this is something we were exploring to try to figure out if we could come up with an easier & more economical library generation workflow. I was disappointed but equally blown away when I realized numerous people had repeatedly “discovered” this fact in the literature already; the most memorable of the names given to it was IVA, or In-Vitro Assembly. Ethan had tried some experiments, and had said it worked roughly as well as with Gibson. Of course, I can’t recall exactly what his experiment was at this point (Although probably a 1-piece, DNA recircularization reaction, since this was in the context of inverse PCR-based library building, after all). So the take away I had was that it was a possible avenue for molecular cloning in the future.

We’ve done a fair amount of molecular cloning in the lab already, creating ~ 60 constructs in the first 4 months since Sarah joined. I forgot exactly the circumstances, but something was right where it made sense to try some cloning where we didn’t add in Gibson mix. I was still able to get a number of intended constructs on that first try, so I stuck to not adding Gibson mix for a few more panels of constructs. I’ve been trying to keep very organized with my molecular cloning pipelines and inventories, which included keeping track of how often each set of mol cloning reactions yielded correctly pieced-together constructs. I’ve taken this data, and broken it down based on two variables: whether it was a 1- or 2-part DNA combination (I hardly ever try more than 2 in a single reaction, for simplicities’ sake, and also because properly combined cloning intermediates may still be useful down the line, anyway), and whether Gibson mix was added or not. Here’s the current results:

So this was extremely informative. Some points
1) I’m willing to screen at least 4 colonies for a construct I really want. Thus, I’m counting a success rate > 0.25 as being a “successful” attempt at cloning a construct. In the above plot, that means any area above the dotted red line. Thus, 1-part DNA recircularizations have pretty decent success rates, since the area of the colored curve above the red dotted like >> the area below it. Sure, Gibson mix helps, but it’s not a night-and-day difference.
2) 2-part DNA combinations are a completely different story,. Lack of Gibson means that I have just as many failed attempts at cloning something as successful attempts. Those are not great odds. Adding Gibson mix makes a big difference here, since it definitely pushes things in favor of a good outcome. Thus, I will ALWAYS be adding GIbson mix before attempting any 2-part DNA combinations.

Other notes: I’m using home-grown NEB 10-beta cells, which give me pretty decent transformation rates (high-efficiency 1-part recircularization reactions can definitely yield many hundreds of colonies on the plate from a successful attempt), so there have been relatively few plates where I literally have ZERO colonies, where I’m more likely to have a few colonies that are just hard-to-remove residual template DNA).

So as I bring in trainees into the lab, I’ll want them to learn how to do some (at the very least) basic data analyses. I realize they may not know exactly where to start, so as I go about making my own basic analysis scripts for analyzing data relevant to me, I’ll post them here and make sure they’re reasonably well commented to explain what is going on.

OK, the specific backstory here. Back in UW, we had next-day IDT orders. It became very clear that was not going to be the case at CWRU, especially after talking to the IDT rep (who did not seem to really care about having a more thriving business here). So, I priced out my options, and Thermo ended up handedly winning the oligo price battle ($0.12 a nucleotide. which is a slight improvement over the $0.15 we seemed to be paying at UW *Shrug*). Thermo also does no shipping cost (another bonus), with the downside being that they only deliver on Tuesdays and Thursdays. I wanted to figure out how long it takes to receive primers are ordering them, so I’ve been keeping track of when I made each order, and how long it took to arrive. Now that I have a decent number of data-points, I decided to start analyzing it to figure out if there were any patterns emerging.

Here’s a link to the data. Here’s a link to the R Markdown file. You’ll want both the data and the R Markdown file in the same directory. I’m also copy-pasting the code below, for ease:

Primer_wait_analysis

This is the first chunk, where I’m setting up my workspace by 1) Starting the workspace fresh 2) Importing the packages I’ll need 3) Importing the data I’ll need

# I always like to start by clearing the memory of existing variables / dataframes / etc.

rm(list = ls())

#Next, let's import the packages we'll need.
#1) Readxl, so that I can import the excel spreadsheet with my data
#2) Tidyverse, for doing the data wrangling (dyplr) and then for plotting the data (ggplot)

library(readxl)
library(tidyverse)

## ── Attaching packages ─────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ ggplot2 3.2.1     ✓ purrr   0.3.3
## ✓ tibble  2.1.3     ✓ dplyr   0.8.3
## ✓ tidyr   1.0.0     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.4.0

## ── Conflicts ────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

# Lastly, use readxl to import the data I want

primer_data <- read_excel("How_long_it_takes_to_receive_items_after_ordering.xlsx")

The idea is to make a bar-graph (or equivalent) showing how many days it takes for the primers to arrive depending on what day I order the primers. I first have to take the raw data and do some wrangling to get it in the format I need for plotting.

# I care about how long it took me to receive primers in a *normal* week, and not about what I encountered during the holidays. Thus, I'm going to filter out the holiday datapoints.

primer_data_normal <- primer_data %>% filter(holiday == "no")

# Next, I want to group data-points based on the day of the week

primer_data_grouped <- primer_data_normal %>% group_by(day_of_week) %>% summarize(average_duration = mean(days_after_ordering), standard_deviation = sd(days_after_ordering), n = n())

# Now let's set the "day_of_week" factor to actually follow the days of the week

primer_data_grouped$day_of_week <- factor(primer_data_grouped$day_of_week, levels = c("M","T","W","R","F","Sa","Su"))

# Since I'll want standard error rather than standard deviation, let's get the standard error
primer_data_grouped$standard_error <- primer_data_grouped$standard_deviation / sqrt(primer_data_grouped$n)

Now that the data is ready, time to plot it in ggplot2.

Primer_plot <- ggplot() + geom_point(data = primer_data_grouped, aes(x = day_of_week, y = average_duration)) +
  geom_errorbar(data = primer_data_grouped, 
                aes(x = day_of_week, ymin = average_duration - standard_error, 
                    ymax = average_duration + standard_error), width = 0.5) + 
  geom_text(data = primer_data_grouped, aes(x = day_of_week, y = 0.2, label = paste("n=",primer_data_grouped$n))) + geom_hline(yintercept = 0) + scale_y_continuous(limits = c(0,7), expand = c(0,0.01)) +
  theme_bw() + xlab("Day primer was ordered") + ylab("Days until primer arrived") + theme(panel.grid.major.x = element_blank())

ggsave(file = "Primer_plot.pdf", Primer_plot, height = 4, width = 6)
Primer_plot

Hmm, probably should have made that figure less tall. Oh well.

So the n values are still rather small, but it looks like I’ll get my primers soonest if I order on Monday or maybe Tuesday. In contrast, ordering on Thursday, Friday, or Saturday give me the longest wait (though well, some of that are the weekend days, which don’t matter as much). Thus, if I have any projects coming up where I have to design primer, it’s probably worth me taking the time to do that on a Sunday or Monday night instead of late in the workweek.

UPDATE 2/22/2020: Uhhh, so the plot above was my first draft attempt, and I have since honed in on the right representation of hte data:

And that’s because average numbers are only somewhat meaningful, and the distribution of frequencies is much more relevant / accurate. Here’s a link to the updated script for generating the plot with the data file at this other link.

Installing Enrich2 / using a Conda environment

I have a couple of new Macs in the lab that need Enrich2 installed. The goal today is to go through the steps of making a Conda environment specifically for running Enrich2 (and installing Enrich2), developed by Alan Rubin.

The Enrich2 repository

The Enrich2 docs

Alan instructions for making a Conda environment specifically for Enrih2

Alan’s instrucitons for doing this is pretty good, if I remember correctly. But, Alan is a seasoned programmer / computational scientist, while people like me are novices and far less familiar with these steps. Furthermore, depending on the specifics of your computer / system, you may get different errors in the installation process. Thus, here’s my interpretation of this process for the benefit of others like me.

0) This supposes that you have already installed an updated version of Anaconda on your Mac. Go back and do this now, if you haven’t done so already.

1) Download the Anaconda Enrich2 environment file and put it somewhere you can access using Terminal.

2A) Update to the newest version of Anaconda, just in case. Hit “y” if prompted.

conda update conda

2B) OK, so now was installing the right version of Pandas it needs. To do this, I first went into Anaconda Navigator and made a new environent called python2, that uses Python 2. Then in terminal, I called:

conda activate python2

2C) This activated python2, so that now the prompt didn’t say “base”, but now it said “(python2)”. Cool. Envinroment activated. I then installed pandas 0.19.2 by typing in:

pip install pandas==0.19.2

2D) OK, this seemed to work too. I got an error when trying to re-run the “enrich2_env.yml file”, so I removed the “=0.19” part form the .yml file and ran “conda env create -f enrich2_env.yml”. This actually seemed to work giving me a list of packages being extracted. Now to actually test it out.

conda activate enrich2

3) Cool, that worked, and the terminal prompt now says “(enrich2)”. Now we’re in business. Next is actually installing Enrich2 now that we’re in the enviroment. Go to the Enrich2 repository, download the file, unzip it, and then move to its directory in terminal. Then run:

python setup.py install

That seemed to work since it didn’t throw any errors.

4) Next is actually trying to run the application. Type in the following:

enrich_gui

Oh Jesus christ. It literally crashes the Finder due to the following errors:

CGSTrackingRegionSetIsEnabled returned CG error 268435459
CGSTrackingRegionSetIsEnabled returned CG error 268435459
CGSTrackingRegionSetIsEnabled returned CG error 268435459
HIToolbox: received notification of WindowServer event port death.
port matched the WindowServer port created in BindCGSToRunLoop

Looks like this may have been a problem with the MacOS operating system. Updating my OS to Catalina and then trying again.

5) OK, MacOS has been updated to Catalina. Now let’s try running Enrich again. (You’ll likely want to run this command after you’ve navigated to the directory with your raw files to simplify the file locating process).

enrich_gui

Awesome. It worked!

PS. I followed these instructions for my second Mac and it worked like a charm. I even went straight for the Cataline update early on and didn’t run into the error in Step4.

TL;DR: Statistics is everywhere, and simulating bottlenecks that happen during routine lab procedures such as dilutions of cells can potentially help you increase reproducibility, and at the least, help you better conceptualize what is happening with each step of an experiment.

I’m still working on getting a cell counter for the lab. In the meantime, we’ve been using an old school hemacytometer to count cells before an experiment. Sarah had used a hemacytometer more recently than me, and knew to dilute the cells from a T75 flask 10-fold to get them into a countable range for the hemacytometer. She said she had performed the dilution by putting 10 ul cells in 90 ul media (and then putting 10 ul of the dilution into the hemacytometer). But as she said this, she asked whether it was OK to perform the dilution as described; a grad student in her previous lab had taught her to do it that way, but a postdoc there said it was a bad idea. My immediate response was that if the cells are sufficiently mixed, then it should be fine. And while that was my gut reaction, I realized that it was something I could simulate and answer myself using available data. Would the accuracy of the count be increased if we diluted 100 ul of cells into 900ul of media, or 1ml of cells into 9ml of media?

Here are the methods (skip if you don’t want to dive in and want to save yourself a paragraph of reading): To me, it would seem the answer to whether the dilution matters depends on how the cells are dispersed in the media / how variable the count is when the same volume is sampled numerous times. Sarah’s standard practice is to count four squares of the hemacytometer, so I had four replicate counts for each volume pipetted. She had repeated this process three times by the time I had performed the analysis, giving me a reasonable dataset I could roll with. I got the mean and standard deviations for each of the three instances, all corresponding to a volume of 0.1 ul. They were all quite similar, so I created a hypothetical normal distribution from the average mean and standard deviation. Next was seeing how different ways of performing the same dilution impacted the accuracy of individual readings. I recreated the 10 ul cells by sampling from this distribution 100 times, 100 ul cells by sampling 1,000 times, and 1 ml by sampling 10,000 times, and taking the mean. I repeated this process 5,000 times for each condition, and looked at how wide each distribution was.

I then turned the counts into concentration (cells / ml):

Instead of stopping there, I thought about the number of cells I was actually trying to plate, which was 250,000. The number the distributions were converging to was ~ 27.3 (black line), so I used that as the “truth”, and saw how many “true” cells would be plated if I had determined the volume needed to be plated based on each of the repeat concentrations calculated by each of the conditions of dilutions. The resulting plot looked like this:

So as you can tell based on the plot, there are slight differences in cells plated depending on imprecision propagated by the manner in which the same 10-fold dilution as performed: while all distributions are centered around 250k, the 10 ul dilution distribution was quite wide, while the 1 ml in 9 ml dilution resulted in cell counts very close to 250k each time. To phrase it another way, ~15% of the time, a 10 ul in 90 ul dilution would cause the “wrong” number of cells to be plated (less than 24k, or more than 26k). In contrast, due to the increased precision, a 100 ul in 900 ul dilution would never result in the “wrong” number of cells being plated. So speaking solely about the dilution, the way the dilution was being performed could have some light impacts on the accuracy of how many cells would be actually plated.

I was going to call this exercise complete, but I ran this analysis by Anna, and she mentioned that I wasn’t REALLY recreating the entire process; sure I had recreated the dilution step, but we would have also counted cells from the dilution in the hemacytometer to actually get the cell counts in real life. Thus, I modified the code such that each dilution step was followed by a random sampling of four counts (using the coefficient of variation determined from the initial hemacytometer readings), and taking the mean of those counts; this represented how we would have ACTUALLY followed up each dilution in real life. The results were VERY different:

In effect, the imprecision imparted by the hemacytometer counts seemed to almost completely drown out the imprecision caused by the suboptimal dilution step. This was pretty mind-blowing for me; especially considering that I would have totally missed this effect had I not run this post by Anna. Now fully modeling the counting process, a 10 ul in 90 ul dilution would cause the “wrong” number of cells (less than 24k, or more than 26k) to be plated ~ 42.5% of the time, and a 100 ul in 900 ul dilution would still cause a “wrong” cell number to be plated ~ 42.2 % of the time; almost identical! Thus, while a 100 ul in 900 ul dilution does impart some slightly increased accuracy, it’s quite minor / negligible over a 10 ul in 90 ul dilution. So while in a sense this wasn’t the initial question asked, it’s still effectively the real answer.

At the end of the day, I think the more impactful aspect of this exercise is the idea that even routine aspects of wet-lab work are deeply rooted in stats (in this case, propagation of errors caused by poor sampling), and that the power of modern computational simulations can be used to optimize these procedures. There’s something truly empowering to having a new tool / capability that gives you new perspectives on procedures you’ve done a bunch of times, and allows you to fully rationalize it rather than relying on advice given to you by others.

Here’s the code if you want to try running it yourself.
Acknowledgements: Thanks to Sarah for bringing this question to my attention. Also, BIG THANKS to Anna for pointing out where I was being myopic in my analysis, which got me to a qualitatively different (and more real-life relevant) answer. It really is worth having smart people look over your work before you finalize it!

UPDATE 5/18/2020: Well, for whatever reason, (at least my version of) Pymol stopped turning the camera (turn command) or molecule (rotate command) in script form, so the below no longer works. Buuutttt. I just told used the drop-down menu to do the movie > program > camera loop > y roll > 4 (or whatever) seconds. And then went to export movie in the file menu and made my movie that way. *shrug*

I’m sure there are many ways to do this, but this is the way I’ve been doing it most recently:

1) Set up your pymol session. That means importing your structure, turning the background to white, and enabling any other setting to make it the desired quality.

2) I use a custom python script to make a series of commands that make the structure turn slightly, ray-trace the structure (for a high quality image), and export the image to a PNG file. The script that I’ve linked to here makes pymol export 360 images, making the structure spin around completely.

Note: Depending on your settings, this process can take 10 minutes or 3 hours. An easy parameter to change would be the resolution (default setting in the script is 2000). Obviously the settings that are toggled during step 1 can drastically change how long it takes as well.

3) Use ffmpeg (at least, on a mac) to turn the 360 images into a video.

$ ffmpeg -framerate 20 -pattern_type glob -i '*.png' -c:v libx264 -preset slow  -profile:v high -level:v 4.0 -pix_fmt yuv420p -vf pad="width=ceil(iw/2)*2:height=ceil(ih/2)*2" -crf 22 -codec:a aac Output.mp4

I used homebrew to install ffmpeg last time I had to do it (brew install ffmpeg). I followed these instructions last time I had to install homebrew.

4) Drag into your powerpoint presentation, and voila!