Modeling bacterial growth

I do a lot of molecular cloning, which means a lot of transformations of chemically competent e.coli. Using 50 uL of purchased competent bacteria would cost about $10 per transformation, which would be an AWFUL waste of money, especially with this being a highly recurring expense in the lab. I had never made my own competent cells before, so I had to figure this out shortly after starting my lab. It took a couple of days of dedicated effort, but it ended up being quite simple (I’ll link to my protocol a bit later on). Though my frozen stocks ended up working fine, I became quite used to creating fresh cells every time I need to do a transformation. The critical step here is taking a saturated overnight starter culture, and diluting it so you can harvest a larger volume of log-phase bacteria some short time later. A range of ODs [optical density here defined as absorbance at 600 nm] work, though I like to use bacteria at an OD around 0.2. I had gotten pretty good at being able to eyeball when a culture was ready for harvesting (for LB in a 250 mL flask, I found this was right when I started seeing turbidity), but I figured there was a better way to know when it’s worth sampling and harvesting.

I started keeping good notes about 1) the starting density of my prep culture (OD of the overnight culture divided by the dilution factor), 2) the amount of time I left the prep culture growing, and 3) the final OD the prep culture. I converted everything into cell density which is a bit more intuitive than OD (I found 1 OD[A600] of my bacteria roughly corresponded to 5e5 bacteria per mL), and worked in those units from there on out. Knowing bacteria exhibit exponential growth, I log base-10 transformed the counts. Much like the increasing number of COVID-19 deaths experienced by the US from early March through early April, exponential growth becomes linear in log-transformed space. I figured I could thus estimate the growth of my prep culture of competent cells by making a multi-variate linear model, where the final density of the bacteria was dependent on the starting bacterial density and how long I left it growing. I figured the lag-phase from taking the saturated culture and sticking it into cold-LB would end up being a constant in the model. Here’s my dataset, and here’s my R Markdown analysis script. My linear model seemed to perform pretty well, as you can see in the below plot. As of writing this, the Pearson’s r was 0.98.

The aforementioned analysis script has a final chunk that allows you to input the starting OD of your starter culture, and assuming a 1000-fold dilution, tells you how long you likely need to wait to hit the right OD of your prep culture. Then again. I don’t think anyone really wants to enter this info into a computer every time they want to set up a culture, so I made a handy little “look-up plot”, shown below, where a lab member could just look at their starter culture OD on the x-axis, choose the dilution they want to do (staying 2x within 1000-fold since I don’t know if smaller dilutions can affect bacterial competency), and figure out when they need to be back to harvest (or at least stick the culture on ice). I’ve now printed this plot out and left it by my bacterial shaker-incubator.

I’m still much more of a wet-lab scientist than a computational one. That said, god damn do I still think the moderate amount of computational work I can do is still empowering.

Gibson / IVA success rates

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:

Note: This is a *stacked* smoothed histogram. Essentially, the only real way to look at this data is consider the width of a given color across the range of the x-axis, relative to its thickness in other portions.

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).

Plasmid Lineages

Recombinant DNA work is integral to what we’re doing here, so I’ve become extremely organized with keeping track of the constructs we are building. This includes having a record of how sequences from two constructs were stitched together to create a new construct. Here’s a network map showing how one or more different plasmid sequences were combined to create each new construct.

[The series of letters and numbers prefixed with G (for Gibson) are unique identifiers I started giving new constructs when it became clear partway through my postdoc that I was going to need a better way of tracking everything I was building. Those prefixed with A are constructs obtained through addgene. Those prefixed with R are important constructs I had built before this tracking system, where I had to start giving them identifiers retroactively.]

HEK293Ts with melanin

I think synthetic biology is really cool, and I like playing around with recombinant DNA elements so I can see how well they work in my own hands. If they work OK, then I just let that knowledge stew in the back of my brain until I can eventually figure out a use for it. Reading this paper by Martin Fussenegger made me realize just how easy it is to make cultured cells express melanin. Here was my first foray in creating melanin in HEK cells by overexpressing tyrosinase/

Cells pelleted in the tubes on the left are expressing tyrosinase. The cells pelleted in the tubes on the right are not.

Doesn’t quite work well enough to use as a general reporter (it’s really hard to tell in a cell monolayer, and only becomes noticeable as colonies of cells or in a pellet, like above), but still kind of fun to see. Let’s see if I find an eventual use for this in some future work.

Primer Inventory Google Sheet To Benchling

Here’s a Python script that converts the MatreyekLab primer google sheet into a csv file that is easily imported into benchling.

1) Go to the lab primer inventory google sheet -> “https://docs.google.com/spreadsheets/d/15RDWrPxZXN34KhymHYkKeOglDYzfkhz5-x2PMkPo0/edit?usp=sharing”

2) Go to file -> download -> Microsoft Excel (.xlsx)

3) Then take above file (MatreyekLab_Primer_Inventory.xlsx) and put it in the same directory as the Google_sheet_to_benchling.py file.

4) Open terminal, go to the right directory, and then enter:

Python3 Google_sheet_to_benchling.py

5) It should make a new file called “Matreyeklab_primers_benchling.csv”. The text in this file can be copy-pasted into benchling and imported into the “Primer” folder.

Uploading the list to Benchling

6) Next, Log onto Benchling, go into the “MatreyekLab” project and into the “O_Primers” folder. Make a new folder named with the date (eg. “20200130” for January 30th, 2020). Once in the folder, select “Import Oligos”, and copy-paste the contents of the csv file (Opened using a basic text editor like “TextEdit” or “Sublime Text 3”. 

7) You’ll have to scroll down all the way to the bottom, but then you can hit insert. Benchling is nice enough to get back to work as it uploads. 

Using the new primer list

8) Once it does finishes uploading, you can go to whatever plasmid map you want to annotate with our current primers. Go to the right-hand side, two icons down to “Primers”. Hit attach existing, add the new folder as the new location, and hit “find binding sites”. Select all of the primers (top check box), and then hit the “Attach Selected Primers” button in the top right.

9) Now click on hte sequence map tab and Voila!, you can see the plasmid map now annotated. Find the primer you want (sequencing or otherwise) and go do some science.

Basic data analysis in R Studio

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.

Writing a simple Python script (for Biologists)

Being able to write short scripts that help you do what you need to do is a very empowering ability for modern-day biologists like me.

In the process of thinking about planning a new plasmid that encodes two proteins fused together by a flexible linker, I realized that I should just write a short script that, when run, gives me a random set of glycine and serine residues composing that linker (riffing off of the ol’ GGS-repeat linkers that seem to prevail in synthetic molecular biology). I’m also in the spirit of documenting the things I’ve learned so that future trainees could have them as reference for their own learning. So here goes my attempt at explaining how I’m approaching this:

First is thinking about the scope of the script you’re trying to put together. This should be an extremely simple one, where I’ll have lists of all of the glycine and serine codons, and a series of random numbers will determine which of those codons to use. Maybe I’ll also include in this script a user-feedback feature that will make it so you can tell the script how many codons it should be randomnly stringing together. Since it’s so simple, it shouldn’t require loading too many custom libraries / packages for performing more advanced procedures.

In the spirit of good practice, I’ll try to do this in iPython as well, though I’ll make a simple python script at the end for quick running from the command line. OK, here goes:

Note: I posted the final code at my Github page (link at the bottom). But for people wanting to just follow along here, here is the code in its final form at the outset (so you can see how I went about building it):

from random import randint

codons = ["GGT","GGC","GGA","GGG","TCT","TCC","TCA","TCG","AGT","AGC"]

length = input("What amino acid length flexible linker would you like a nucleotide sequence for?")

linker_sequence = ""
length = int(length)
for x in range(0,length):
    new_codon = codons[randint(0,9)]
    linker_sequence = linker_sequence + new_codon
    x = x + 1

print(linker_sequence)

1) Had to remind myself of this, but firs you type in “jupyter notebook” to open the interactive web-browser interface (assuming you’ve already installed it).

2) Go to the right directory, and then make a new iPython3 notebook.

3) Let’s start simple and make the random codon generator. First, let’s make a list of the codons we want to include.

codons = ["GGT","GGC","GGA","GGG","TCT","TCC","TCA","TCG","AGT","AGC"]

4) Next, let’s figure out how to choose a random index in the list, so that a random codon is chosen. I don’t remember how to do this off-hand, so I had to google it until I got to this page.

5) OK, so i guess I do have to load a package. I’m now writing that at the top of the script.

from random import randint   
randint(0,9)

It gave the correct output so I saw that it worked and commented out the randint call.

# randint(0,9) #This line won't run

6) Next is having the script randomly pull out a glycine or serine codon. This is simply done by now including:

codons[randint(0,9)]

7) Cool. All that has worked so far. Let’s now have the script take in a user-input for the number of residues so we can repeat this process and spit out a linker sequence of desired length. I don’t remember how to ask for user input in a python script, so I had to google this as well.

8) Allright, so I need to use the use the “raw imput” function and lead it to a variable. Well, tried that and it said “raw_input” wasn’t defined. Google that, and found this link saying that advice as deprecated, and instead it was just “input()”.

9) Thus, I typed in:

length = input("What amino acid length flexible linker would you like a nucleotide sequence for?")

And this asked me for a number like I had hoped. Great.

10) To make this process iterable, I went for my trusty “for” loop. Actually, I tried to make a “while” loop first, but realized I didn’t know how to make it work off the top of my head. So I just ended up making a “for” loop and putting a x = x + 1 statement at the end to effectively turn into a “while” loop. Personally, I think the ends justify the means, and going with what you know works well is a valid option for most basic scripts, when efficiency isn’t a huge priority.

length = 3  #Giving an arbitrary number for testing the script
for x in range(0,length):
    new_codon = codons[randint(0,9)]
    x = x + 1 

11) I then had to be able to keep track of the codons that were pulled during each iteration of the loop. I thus created an empty variable called “linker_sequence” and just added the string for the new codon at the back of “linker_sequence” during each iteration.

linker_sequence = ""  # A blank variable to keep track of things
length = int(length)  # To convert the input text into a number
for x in range(0,length):   
    new_codon = codons[randint(0,9)]  
    linker_sequence = linker_sequence + new_codon  
    x = x + 1

12) Lastly is putting a final print function so it returns the desired string of nucleotides to the user.

print(linker_sequence)

Nice, I think that does it for the script. Super easy and simple!

13) Finally, let’s test the script. I typed in 5 amino acids, and it gave me “TCTAGTTCAGGCTCT” as the output string. I google searched “Transeq” to get to the ebi for a simple codon translator, and translation of the nucleotide sequence above gave “SSSGS” as the protein sequence. So great, it worked! Sure, a little serine heavy, but that’s random chance for you. Should still be perfectly fine as a linker, regardless. Now to finish planning this plasmid…

Note: I’ve posted both the iPython notebook file (.ipynb) and a simple python script (.py) on my Github page. Feel free to use them!