(20 min read) Creating a CRISPR-Cas3 Phage Vector to Knock Out the Essential Genes of E.
For a while, I thought phage therapy was only good for curing bacterial infections that have gained antibiotic resistance. Although that…
Creating a CRISPR-Cas3 Phage Vector to Knock Out the Essential Genes of E. coli CFT073
For a while, I thought phage therapy was only good for curing bacterial infections that have gained antibiotic resistance. Although that on its own sounds really good, I realized that there is so much more to that and that phage therapy has so many other possibilities.
Phage therapy has been researched to transport other gene materials/medicines besides their own, become vectors for cancer treatments, detect bacteria in the body and other substances (like water or urine, for example), etc.
I was specifically interested in the phage therapy X gene editing possibility, because, to me, it sounded like gene editing could actually make phage therapy a reality. I thought it meant that it would make the FDA approval of phage therapy much easier and much faster.
I’m actually late to the game. So many large companies are working with phages and gene editing them for different purposes. But one company that stood out to me the most was Locus Biosciences, led by Paul Garofolo.
Locus is working to gene-edit phages to help overcome FDA approval barriers like treatment safety/consistency, phage host range, phage stability, and so on. But they are also using gene editing on phages in a different way.
Locus is gene editing phages to carry their own CRISPR systems. Those phages then act as vectors that are able to deliver those CRISPR systems to bacteria or other hosts. Once the CRISPR systems are inside the host, it’s then able to modify or suppress the genes of their host. If you don’t know what CRISPR is, read here.
When I read about Locus creating CRISPR-phage vectors, I thought that this could definitely help solve the barriers of consistency and safety between phage treatments. So, I decided to try and build one of my own
[embed]A video of me explaining my CRISPR-Cas3 phage vector system.
What Did I Do?
I built my own CRISPR-Cas3 phage vector to deliver a CRISPR-Cas3 system to the E. coli strain CFT073. Once inside E. coli CFT073, the CRISPR system will knock out its essential genes and speed up the bacteria elimination process.
I did this by selecting my target genes, building the CRISPR system, and editing the gene into the chosen phage to create a phage vector. Then, I ran three simulations to test the efficacy of my vector when trying to knock out entire populations of E. coli CFT073.
The entire project was built and tested online using biocomputing software and writing my own simulation code.
Why Did I Do It?
I learned that one of the problems with implementing phage therapy is due to a matter of safety and biological variability between phage species, and even bacteriophage to bacteriophage. Gene-editing a phage treatment batch to all carry the same kind of CRISPR-Cas3 system makes the entire treatment more consistent and reliable.
My ultimate end goal with all the projects and research I do is to try and find a way to get phage therapy FDA-approved, or at least one step closer. This project does just that because I’m exploring something already being done by companies related to phage therapy FDA approval, but putting my own spin on it.
Early-Stage Project Setup
Choosing Bacterial Strain E. coli CFT073
E. coli is known for developing antibiotic resistance really quickly, and at a very strong level. When focusing on developing different phage therapy treatments, developing a treatment that is E. coli-specific will have more impact and can be better adopted and implemented for other bacterial strains.
I wanted to choose E. coli for those reasons because if I am able to develop something that works for E. coli, a dangerous and persistent pathogen, it will have more real-world impact and be versatile for other situations.
I chose the specific strain CFT073 because I had looked at the strains Locus Biosciences was working with and did some of my own research, and CFT073 isn’t a “popular” strain to work with, so I thought I’d try to work with something more niche and less researched.
Electron microscopy of E. coli CFT073. Credit
Deciding to Work With CRISPR-Cas3
CRISPR-Cas3 is used less compared to Cas9 or Cas12a, but Locus was working with Cas3, and I wanted to work with it as well. But, because it was used less, I was also wondering why a large company like Locus was investing so much time into developing phage treatments based on this system.
It all made sense when I was able to compare Cas9 to Cas3. Cas9 uses one protein. It’s like a pair of scissors, which is great to use for small, precise cuts. Cas3, on the other hand, is made up of six different Cas proteins. This system is like a shredder, which is able to take out larger portions of DNA efficiently, and multiple portions in one run.
Because, for this project, we are trying to knock out the essential genes of E. coli, we aren’t worried about preserving the bacteria’s genomic integrity. If something gets cut that wasn’t intended to, it’s not a big deal because we aren’t trying to keep the host alive. So, it’s better to use Cas3, which can cut larger portions of DNA, and multiple portions of it (so multiple locations), even if we are running the risk of it not being as precise. If we wanted to make specific cuts (at one location), where being precise was very important, then we would use Cas9.
Process of CRISPR-Cas3 gene editing. Credit
Choosing the Five Target Genes Essential to E. coli CFT073
When doing research on what some general essential genes are for E. coli, I found they had these main functions:
- Maintaining genomic integrity and structure
- Repairing damaged DNA
- Functions for transcription and translation of DNA to proteins
So, I focused on choosing genes that had these common functions. But, I also chose one of my genes to inhibit population control and prevent the E. coli from multiplying. The five target genes I chose are:
- ftsZ: Involved in bacterial cell division. Without this gene, the bacteria cannot properly divide and replicate themselves. This prevents it from expanding its population.
- gyrA: DNA is structured in a coil shape, and this gene is what helps the DNA keep that structure. If the gene is knocked out, it disrupts DNA supercoiling, causing the structure to unravel and creating fatal DNA damage.
- murA: This gene is responsible for cell wall synthesis. Some bacteria, like CFT073, have a cell wall, which is a thick substance that surrounds the bacteria, protecting them, as well as helping them maintain their outer structure. Without the cell wall, the structure begins to fall apart, leading to bacterial rupture.
- recA: recA is one of the “SOS-response” genes, using homologous recombination, that helps to immediately repair compromised DNA material, as well as to help maintain genomic integrity within the CFT073 genome. Without it, the genome is more susceptible to damage and can’t be repaired easily, which can cause fatal damage to the bacteria.
- rpoB: Involved in the transcription and translation process of DNA to proteins. Transcription is DNA to RNA, and translation is RNA to proteins. Those proteins are then essential to carry out different tasks within the host and maintain stability and longevity. rpoB is necessary in the transcription process, from DNA to RNA. Without it, RNA can’t be synthesized from DNA and can’t go on to make the essential proteins needed for the bacteria to function, which can compromise the bacteria and endanger it.
Building the CRISPR-Cas3 Complex
Finding the gRNAs for my CRISPR-Cas3 System
The gRNAs are what’s going to guide the Cas3 complex to the location for gene deletion. Without them, the Cas system has no idea where the genes are for deletion and won’t be able to perform its job.
A gRNA is a portion of the DNA you want to cut, which is transcribed into RNA (turned into RNA from DNA, like in the transcription and translation process). For Cas3, a gRNA is roughly 32 base pairs (bp) long. For Cas9, it’s roughly 20 bp.
One of the things that determines if a gRNA will be effective is if it’s compatible with a protospacer adjacent motif (PAM) sequence. A PAM sequence is a short DNA sequence that appears immediately next to the target DNA sequence that a CRISPR system wants to cut. Cas proteins need the PAM site in order to recognize it and effectively bind to the target DNA. For Cas3, the optimal PAM sequence is GAA (or TTC on the reverse strand). A gRNA is extracted 32 bp upstream from where the PAM sequence is found in the target DNA.
There are many software types online that allow you to paste a copy of your target gene and select your target host, and then they’re able to find the most compatible gRNAs based on what CRISPR type you’re using. But, Cas3 isn’t a commonly used Cas protein, so a lot of the current software out there doesn’t support Cas3 as an option for gRNA search.
Because of this roadblock, I was having a lot of trouble moving forward with my project. I was looking for different software to use to find my gRNAs, but nothing was turning up. So, eventually, I decided to create my own software to find my own gRNAs compatible with Cas3.
I used AI-powered simulation/code-building software to help me build software that could accurately find the gRNAs I needed for my CRISPR-Cas3 system. Getting into the code for the gRNA software can get technical, so you can read more about it here.
The chosen gRNAs are:
- ftsZ gRNA Sequence: CGCGCCACGCTATACAAACATTTACGGGAGT PAM: GAA Location: 1092–1124 Mismatches: 0 GC Content: 0.53 Score: 0.9215
- gyrA gRNA Sequence: GCGCAAATGCCTGTTTCAACTGCCGTTGCG PAM: GAA Location: 203–235 Mismatches: 0 GC Content: 0.60 Score: 0.8945
- murA gRNA Sequence: TGCAATTTGCAGTGCGCTGCGATGAACTCA PAM: TTC Location: 187–219 Mismatches: 0 GC Content: 0.50 Score: 0.9387
- recA gRNA Sequence: TGGGAAAACCCACGTCACCAGCAGGCTCTT PAM: TTC Location: 238–270 Mismatches: 0 GC Content: 0.57 Score: 0.9102
- rpoB gRNA Sequence: TGCAATTTGCAGTGCGCTGCGATGAACTCA PAM: TTC Location: 187–219 Mismatches: 0 GC Content: 0.50 Score: 0.9387
Getting the Cas3 Genes for the Complex
I now had to collect the Cas3 genes for the Cas3 proteins needed to create the complex. CRISPR was originally discovered in bacteria as a defence mechanism against viruses like phages, so normally, a lot of bacteria are already equipped with the entire complex. If the strain already had its own CRISPR-Cas3 complex, I would only have to transport the gRNAs. But, E. coli CFT073 didn’t have its own CRISPR-Cas3 system, so I had to create the entire construct and have the phage transport the whole thing.
I took the Cas3 genes from the E. coli K12 strain and then used them to put together the complex. As I have said before, Cas3 is made up of six different Cas proteins, which are (in order for complex design):
- Cas3
- Cas8e
- Cas11e
- Cas7
- Cas5
- Cas6e
I got the genes for each of the Cas proteins and started building the construct by putting in the corresponding genes in that order.
Gathering the Other Components for a Functional CRISPR Construct
There are other components needed for the full construct, like a promoter, terminator, and spacer sequence.
A promoter “activates” the entire CRISPR-Cas3 system to start looking for the matching gRNA sequences. I used the Anderson promoter BBa_J23119: TTGACGGCTAGCTCAGTCCTAGGTACAGTGCTAGC.
This promoter is widely used in synthetic biology and works well in E. coli hosts. This promoter is considered a constitutive promoter (it’s always on), meaning that we don’t need to add another component that triggers it to turn on. Using a constitutive promoter is simpler because once the phage infects the bacteria and delivers the CRISPR complex, the promoter can activate it right away because it’s always on.
A terminator tells the RNA to stop transcription to prevent run-on and to keep genes organized. I’m using the terminator BBa_B0015: CCGCTGAGCAATAACTAGCATAACCCCTTGGGGCCTCTAAACGGGTCTTGAGGGGTTTTTTG.
This terminator has been tested and validated in real-world projects and is also highly compatible with E. coli hosts.
The last main component we need is a repeat. A repeat is also classified as a CRISPR Repeat RNA (crRNA). I am using a natural CRISPR repeat that is generally used for published Cas3 work: GTGTTCCCCGCGCCAGCGGGGATAAACCG.
A repeat has a few different purposes. It acts as a barrier between each of the gRNA sequences in the construct to prevent any from getting evaluated or used together. Cas6e is the protein responsible for cutting the entire CRISPR-Cas3 array into different portions, and it’s the protein that is able to identify and read the spacer sequence.
The entire complex is constructed in the order:
- Promoter
- Cas genes (in order shown before)
- Repeat, then gRNA (x5, for each of the target genes)
- Terminator
To see what the entire CRISPR-Cas3 complex looks like put together, you can see it here.
Creating the Phage Vector
Selecting the Killian Phage as the Vector
The phage I chose to be my phage vector was the Killian phage. It’s a less-researched phage compared to others like T4 or T7 phages (I was considering T4 phage as my vector). I ended up choosing the Killian phage because it has a large genome and is specific to E. coli CFT073. The T4 phage, while also really good for carrying large insertions and acting as a vector, was not specific to strain CFT073, but rather to general kinds of E. coli (like K12), and it wasn’t specified whether or not it infected CFT073 as well, which was less likely.
Phages attacking a bacterium. Credit
Putting the CRISPR-Cas3 Complex Into a Plasmid
To get the CRISPR system into the phage genome for it to successfully carry it and infect the bacteria, the CRISPR-Cas3 system needed to be packaged and made more secure. I had to get a plasmid backbone so I could insert my CRISPR system into it, and that plasmid would act like a carrier that could be inserted into the phage genome.
A plasmid backbone provides the necessary functions for replication, survival in the host cell, and expression of the inserted genes. I chose the backbone pHAGE-GFP. This plasmid has a GFP reporter (Green Fluorescent Protein), which is a protein that glows under green UV or blue light. This acts as a visual indicator as to whether the plasmid was inserted successfully into a cell and if the Cas genes are being expressed.
The plasmid also includes an ampicillin-resistant gene (AmpR). It helps break down ampicillin, a common antibiotic. This is important for wet lab experiments to help identify which bacteria were infected by your plasmid.
- You would expose all the tested bacteria to ampicillin in petri dishes.
- After inserting the plasmid into bacteria, only cells with the plasmid survive on ampicillin plates because they have the resistant gene.
- Bacteria without the plasmid get killed off, so only the ones carrying and expressing your CRISPR-Cas3 system survive in culture.
In the plasmid backbone, I inserted the CRISPR-Cas3 system between two restriction enzymes:
- XhoI: CTCGAG
- NheI: GCTAGC
Restriction enzymes act as a designated zone for different inserts inside a plasmid backbone, protecting the rest of the important parts of the backbone from being disrupted. Both of the plasmids create “sticky ends” on either side of the insert that help bind with the DNA of the insert. The enzymes I used in the plasmid backbone are well characterized, meaning they are reliable as an insert location for this specific plasmid.
You can see the constructed plasmid here.
Inserting the Plasmid Into the Phage Vector
The entire plasmid (with the CRISPR insert) is 15393 bp long. The unedited Killian phage genome was 169905 bp long. I had to first clear any unnecessary parts of the Killian phage genome to make space for the length of my plasmid insert. To do this, I deleted:
- Unannotated space
- Anything labelled “Hypothetical Protein”
- Some endonuclease or endonuclease-like genes
After going through the genome and doing a sweep to clear things, I had made enough space to insert my plasmid into the Killian phage. But now I had to find a space in the genome that was safe to disrupt, and preferably away from any necessary genes. This task was also one of my harder tasks, since it took a lot of patience and time to make sure I did it right.
I disrupted an hnh endonuclease gene, which is a gene commonly used for insertions and can be disrupted without causing major problems to the phage. Once the insert was in the Killian phage, I had basically finished my phage vector, and I could now move on to testing.
You can see the edited Killian phage genome here.
Running the Simulations
If you want to see in detail the codes and parameters for each of the simulations run, you can read from here.
How Effective is the Killian Phage at Delivering CRISPR-Cas3 to E. coli CFT073?
This was the first simulation I ran to test my phage vector. I wanted to see how many of the bacteria from the population were actually able to receive the CRISPR-Cas3 plasmid from the Killian phage.
This is the graph measuring the amount of uninfected, infected without CRISPR, and infected with CRISPR bacteria across 60 minutes after initial infection begins.
What You’re Seeing
This graph shows the results of running the simulation 5 times with slight variations in parameters to account for biological variability:
Solid Lines: The average (mean) values across all simulations
Shaded Areas: The standard deviation, showing the range of possible outcomes
Key Takeaways
- Infection Rate: The phage infects bacteria rapidly, with approximately 98.5% of bacteria infected within 60 minutes at the specified MOI of 10:1.
- CRISPR-Cas3 Delivery: Of the infected bacteria, about 85% successfully receive the CRISPR-Cas3 system, resulting in approximately 83.7% of the total bacterial population containing the CRISPR-Cas3 system after 60 minutes.
- Efficiency: The simulation suggests that the modified Killian phage is an effective delivery vehicle for the CRISPR-Cas3 system to E. coli CFT073.
How Effective is the CRISPR-Cas3 System at Knocking Out the Essential Genes of E. coli CFT073?
In this simulation, I wanted to test the actual efficacy of the Cas3 complex at degrading the essential genes, which was going to be very important for this treatment to work.
This is the graph that shows the cleavage timeline and results of the CRISPR-Cas3 system for each of the essential target genes in E. coli CFT073.
What You’re Seeing
X-axis: Time in minutes post-infection (122–138 minutes shown)
Y-axis: Target genes (ftsZ, gyrA, murA, recA, rpoB)
Red horizontal lines: Duration of cleavage for each gene
Blue dots: Cleavage start times (labelled with exact minutes)
Green dots: Cleavage completion times (labelled with exact minutes)
Key Takeaways
- The graph reveals the sequential nature of gene targeting by the CRISPR-Cas3 system
- The total gene disruption process spans approximately 15.5 minutes from start to finish
How Effective is the Entire Treatment Process (From Initial Infection to Gene Knockout) at Degrading the Populations of E. coli CFT073?
This final simulation was to test the entire treatment process using the CRISPR-Cas3 phage vector, taking the results and parameters from both simulations 1 and 2 to do one final, cohesive test.
This graph shows the clearance rate of the bacteria over a 25 hour span from the beginning of the infection process.
What You’re Seeing
X-axis: Time in hours (0–25 hours shown)
Y-axis (left, blue): Bacterial load in CFU/mL (logarithmic scale)
Y-axis (right, red): Phage concentration in PFU/mL (logarithmic scale)
Blue line: Total bacterial load
Red line: Phage concentration
Key Takeaways
- Clearance begins immediately after phage administration
- 90%, 99%, and 100% clearance are all achieved within 0.5 hours (30 minutes)
- All simulation runs show similar clearance patterns with minimal variability, which is why there is no noticeable standard deviation to be shown
- Complete clearance is maintained throughout the 24-hour simulation period
This graph shows the different phases that the bacterial cells go through from the beginning of infection till they are eliminated.
Detailed E. coli CFT073 Population Dynamics
This high-resolution version includes additional annotations marking key events like the phage infection peak and CRISPR-Cas3 activation, with more detailed explanatory text.
Time to 90% clearance: 22 minutes
Time to 99% clearance: 29 minutes
Time to 100% clearance: 30 minutes
Maximum clearance: 100%
Regarding Stochastic Variability Across the Multiple Simulation Runs
The system is so efficient that:
- All five simulation runs achieved 100% clearance at almost exactly the same time (0.5 hours)
- The biological variability that would normally exist is overwhelmed by the system’s efficiency
- The transition from 0% to 100% clearance is extremely rapid and consistent
Key Takeaways From This Project
Was the Project “Successful”
The vector I built to knockout E. coli CFT073 was successful. The phage vector itself was able to transport the CRISPR-Cas3 system into CFT073. The CRISPR-Cas3 system, once inside CFT073, was able to knock out each of the target genes successfully, with the entire process from start to finish taking roughly 15 minutes. When the entire treatment simulation was run, E. coli CFT073 was killed in roughly 30 minutes, with consistent results across 5 simulation runs with biological variability.
Unexpected Patterns and Behaviours
The fact that the final simulation run had roughly “perfect results” across the different runs was unexpected. But I made sure there weren’t any bugs or mistakes in the code, and it just proved that the system worked really well.
Another thing that was unexpected was that the E. coli died 30 minutes into the simulation, but the CRISPR-Cas3 system hadn’t kicked in until after the population died out. After a little more research, I found that the natural phage infection acted much quicker than the CRISPR-Cas3 system, so the E. coli was dying off because of that. At first, I thought that meant the phage vector I built wasn’t actually working, and that there was no point to it. But what the new CRISPR payload does is it acts like an extra line of defence. Because there’s a chance that not all the bacteria will die out, or that the population will build up again, the CRISPR-Cas3 system shuts down any of those possibilities by knocking out the essential genes and making the bacteria non-functional and damaged, providing extra reassurance that the bacteria are eliminated.
Besides that, the CRISPR-Cas3 system is a way to make all the phages consistent and uniform throughout treatment, because all the phages will be carrying the same payload. It will reduce any variability in the treatment, so having the payload is still necessary.
Comparing to Real-World Data
Locus Biosciences — CRISPR-Cas3 Phage Therapy Against E. coli and Klebsiella
Locus Biosciences is creating phage vectors to carry CRISPR-Cas3 payloads that are able to target both E. coli and Klebsiella bacteria. This is where I based my project off of. Locus has been able to take their CRISPR phage vectors to human trials, currently at phase 1b.
Zhang et al. 2022 — CRISPR-Cas9 Phages Targeting E. coli K-12 and Pathogenic Strains
In this research, they created a modified T4 phage vector to deliver CRISPR-Cas9 systems to E. coli K-12. The Cas9 systems were targeting antibiotic resistance genes as well as essential genes for E. coli. With the modified vector, they found that there was efficient killing of the bacteria and less off-target risk.
Yosef et al. 2015 — CRISPR Phage Re-Sensitization of E. coli
This research designed a phagemid carrying CRISPR-Cas9 to cut antibiotic-resistant genes in resistant E. coli and Staphylococcus. A phagemid is a plasmid that contains elements from both a plasmid and a phage. It’s a small, circular DNA molecule that replicates inside a bacterium, like a plasmid, and can be packaged into a phage. The key finding was that it proved they could re-sensitize bacteria to antibiotics by deleting the antibiotic-resistant genes they developed.
What’s Next?
There are a few next steps for this project that I’m excited to explore. I’m going to refine the code that I used to identify the required gRNAs for this project. I want to take the code and turn it into something everyone can use, and see if others will find it useful as well for their research. I’ve also been developing a modified phage vector that will function the same way as the one I created for this project, only it will be a T4 phage made to infect the E. coli strain B. I’m making this modified vector so that I can test my phage vector in the lab using non-pathogenic substances.
Phage therapy X gene editing has a lot of potential, and more research like mine and like the many other people investing their time in this area are what’s going to get this treatment FDA-approved and accessible for those who need it.
I’m Isavella
I’m 15 years old and am currently researching phage therapy and the obstacles in its implementation. Phage therapy has been known for over 100 years, but is only being researched now, and that unknown potential is what interests me most. In my spare time, I love to write, do Muay Thai, sing, and learn something new. Feel free to connect with me on LinkedIn, visit my website, or subscribe to my newsletter!
How Does CRISPR Work (Explained Very Briefly)
For those who don’t know, CRISPR is one of the most highly researched and used methods for gene editing. The most commonly used CRISPR type is CRISPR-Cas9.
All CRISPR systems have 2–3 main components. The 2 that are always found are the Cas protein and the guide RNA (gRNA). If someone is using CRISPR to replace a current gene, then the third component would include a template DNA sequence for replacement.
The gRNA is what allows the CRISPR system to find the target gene, and the Cas protein is what makes the actual modification of the gene. The Cas protein will make a cut at the target site, and the loose ends will either be filled in by the template DNA or they will repair themselves, causing the deletion of the cut gene (since nothing replaced it).
What makes different CRISPR systems different is based on what Cas protein they use. The most commonly used Cas protein is Cas9. But another popular one is Cas12a.
Cas3 is actually not one but 6 different Cas proteins, which are very different in structure compared to the singular Cas9 or Cas12a.
gRNA Identifier Software Code
Here is the entire code I used to find the gRNA:
#!/usr/bin/env python3.11
import re
import csv # Keep for potential future use, though current output is print
# --- Configuration Parameters ---
VALID_PAMS = ["GAA", "TTC"]
GRNA_LENGTH = 32 # Standard gRNA length, can be adjusted
MAX_OFF_TARGET_MISMATCHES = 3 # For evaluating specificity
# Gene sequences provided by the user
GENE_SEQUENCES = {
"ftsZ": "TTAATCAGCTTGCTTACGCAGGAATGCTGGGATATCCAGATAATCCGGCTCTTTCGCAGTTTGCGGCGCATTGTCATTCACGACTTTAGCAACCGGCTTCTGCTCCTGGGTCAGCGGAGCCATCCCGTGCTGCTGGTAGCGATCCATCACTGGCTGCTGAACCTGCTTATTAGTCACCAGAGTGATTTCAGGACGTTTGTCCATGCCGATACCTGTCGCAACAACAGTTACGCGCAGCTCGTCATTCATATCCGGGTCAAGAGAAGTACCGATAACCACGGTCGCGTTGTCGGAAGCAAATGCACGGATGGTGTTACCTACCGTTTCGAACTCATCCAGACGCAGGTCGAAGCCCGCTGTGATGTTTACCAGTACGCCGCGTGCGCCAGACAGGTCGATATCTTCCAGCAGCGGAGAAGAGATAGCCATTTCAGCAGCTTCTTCCGCACGATCTTCACCGCTCGCCACACCAGAACCCATCATTGCGTAGCCCATTTCAGACATTACGGTGCGTACGTCTGCAAAGTCCACGTTCATCAAACCCGGACGAGTAATCAGTTCGGCGATACCTTGCACAGCGCCTTTTAGTACATCGTTCGCTGCGCCAAACGCATCCAGCAGGGAGATACCACGGCCCAGAACTTTCAGCAGTTTGTCGTTCGGGATAGTGATCAGAGAGTCCACATGCTTGGACAGTTCAGTGATCCCCTGCTCCGCGAATGCCATACGCTTCTTGCCTTCAAAGTTGAAAGGCTTAGTGACGACAGCAACGGTCAGGATACCCAAATCTTTTGCCACTTCAGCGACGACTGGTGCTGCACCTGTACCGGTACCACCACCCATACCCGCAGCAATAAAGACCATGTCTGCACCTTCCAGCGCCGCACGCAGTGCATCGCGATCCTCATCAGCCGCATTGCGGCCAACTTCTGGATTAGCGCCAGCGCCCAGTCCTTTGGTGATACCGCTACCGATTTGAATCGTCTGTCCAACCGCTGTTTTACGCAGCGCTTGTGCATCGGTATTTACCGCGAAGAATTCAACACCTTCAATGCGCTCACGCACCATGTGTTCAACAGCATTACCGCCGCCGCCGCCGACGCCGATGACTTTAATCACCGCGTCATTGGTAAGTTCCATTGGTTCAAACAT",
"gyrA": "ATGAGCGACCTTGCGAGAGAAATTACACCGGTCAACATTGAGGAAGAGCTGAAGAGCTCCTATCTGGATTATGCGATGTCGGTCATTGTTGGCCGTGCGCTGCCAGATGTCCGAGATGGCCTGAAGCCGGTACACCGTCGCGTACTTTACGCCATGAACGTACTAGGCAATGACTGGAACAAAGCCTATAAAAAATCTGCCCGTGTCGTTGGTGACGTAATCGGTAAATACCATCCCCATGGTGACTCGGCGGTTTATGACACGATCGTCCGTATGGCGCAGCCATTCTCGCTGCGTTACATGCTGGTAGACGGTCAGGGTAACTTCGGTTCCATCGACGGCGACTCTGCGGCGGCAATGCGTTATACGGAAATCCGTCTGGCGAAAATTGCCCATGAACTGATGGCCGATCTCGAAAAAGAGACGGTCGATTTCGTTGATAACTATGACGGTACGGAAAAAATTCCGGACGTCATGCCAACCAAAATTCCTAACCTGCTGGTGAACGGTTCTTCCGGTATCGCCGTAGGTATGGCAACCAACATCCCGCCGCACAACCTGACGGAAGTCATCAACGGTTGTCTGGCGTATATCGATGATGAAGACATCAGCATTGAAGGGCTGATGGAACACATCCCGGGGCCAGACTTCCCGACGGCGGCAATCATTAACGGTCGTCGCGGTATTGAAGAAGCTTACCGTACCGGTCGCGGCAAGGTATATATCCGCGCCCGTGCCGAAGTAGAAGTTGACGCCAAAACCGGACGTGAAACCATTATCGTCCACGAAATTCCTTATCAGGTGAACAAAGCGCGCCTGATCGAGAAGATTGCGGAACTGGTAAAAGAAAAACGCGTGGAAGGCATCAGCGCGCTGCGTGACGAGTCTGACAAAGACGGTATGCGCATCGTGATTGAAGTGAAACGCGATGCGGTCGGTGAAGTTGTGCTCAACAACCTCTACTCCCAGACCCAGTTGCAGGTTTCTTTCGGTATCAACATGGTGGCATTGCACCATGGTCAGCCGAAGATCATGAACCTGAAAGACATCATCGCGGCGTTTGTTCGTCACCGCCGTGAAGTAGTGACCCGTCGTACTATTTTCGAACTGCGTAAAGCTCGCGATCGTGCTCATATCCTTGAAGCATTAGCCGTGGCGCTGGCGAACATCGACCCGATCATCGAACTGATCCGTCATGCGCCGACGCCTGCAGAAGCGAAAACCGCGCTGGTTGCAAATCCGTGGCAGCTGGGCAACGTTGCTGCGATGCTGGAACGTGCGGGCGACGATGCTGCGCGTCCGGAATGGCTGGAGCCAGAGTTCGGCGTGCGTGATGGTCTGTACTACCTGACCGAACAGCAAGCTCAGGCGATTCTGGATCTGCGTCTGCAGAAACTGACCGGCCTTGAGCACGAAAAACTGCTCGACGAATACAAAGAGCTGCTGGATCAGATTGCGGAACTGTTGCGTATTCTTGGTAGCGCCGATCGTCTGATGGAAGTGATCCGTGAAGAGCTGGAGTTGGTTCGTGAACAGTTCGGTGACAAACGTCGTACCGAAATCACCGCTAACAGCGCCGACATCAACCTGGAAGATCTGATCACCCAGGAAGATGTGGTGGTCACGCTTTCTCACCAGGGTTACGTGAAGTATCAGCCGCTTTCTGAATACGAAGCGCAGCGTCGTGGCGGGAAAGGTAAATCTGCCGCACGTATTAAAGAAGAAGACTTTATCGACCGACTGCTGGTGGCGAACACTCACGACCATATTCTGTGCTTCTCCAGCCGTGGCCGCGTCTATTCGATGAAGGTCTACCAGTTGCCGGAAGCCACTCGTGGCGCGCGCGGTCGTCCGATCGTCAACCTGCTGCCGCTGGAACAGGACGAACGTATCACCGCGATCCTGCCGGTGACCGAGTTTGAAGAAGGCGTGAAAGTCTTCATGGCGACCGCTAACGGTACTGTGAAGAAAACCGTCCTCACCGAGTTCAACCGTCTGCGTACCGCCGGTAAAGTGGCGATCAAACTGGTTGAAGGCGATGAGCTGATCGGCGTTGACCTGACCAGCGGCGAAGACGAAGTAATGCTGTTCTCCGCTGAAGGTAAAGTGGTGCGCTTTAAAGAGTCTTCTGTCCGTGCGATGGGCTGCAACACCACCGGTGTTCGTGGTATTCGTTTAGGCGAAGGCGATAAAGTCGTCTCTCTGATCGTGCCTCGTGGCGATGGCGCAATCCTCACCGCAACGCAGAACGGTTATGGTAAACGTACCGCAGTGGCGGAATACCCAACCAAGTCGCGTGCGACGAAAGGGGTTATCTCCATCAAGGTTACCGAACGTAACGGTTTAGTGGTTGGCGCGGTGCAGGTAGATGACTGCGACCAGATCATGATGATCACCGATGCCGGTACGCTGGTACGTACTCGCGTTTCGGAAATCAGCATCGTAGGCCGTAACACCCAGGGCGTAATCCTTATCCGTACTTCGGAAGATGAAAACGTGGTGGGTCTGCAACGTGTTGCTGAACCGGTTGACGAGGAAGATCTGGATACCATCGACGGCAGTGCCGCGGAAGGGGACGATGAAATCGCTCCGGAAGTGGACGTTGACGACGAGCCAGAAGAAGAATAA",
"murA": "ATGGATAAATTTCGTGTTCAGGGGCCAACGAAGCTCCAGGGCGAAGTCACAATTTCCGGCGCTAAAAATGCTGCTCTGCCTATCCTTTTTGCCGCACTACTGGCGGAAGAACCGGTAGAGATCCAGAACGTCCCGAAACTAAAAGACGTCGATACATCAATGAAGCTGCTAAGCCAGCTGGGTGCGAAAGTAGAACGTAATGGTTCTGTGCATATTGATGCCCGCGACGTTAATGTATTCTGCGCACCTTACGATCTGGTTAAAACCATGCGTGCTTCTATCTGGGCGCTGGGGCCGCTGGTAGCGCGCTTTGGTCAGGGGCAAGTTTCACTGCCTGGCGGTTGTACGATCGGCGCACGTCCGGTTGATCTACACATTTCTGGCCTCGAACAATTAGGCGCGACCATCAAACTGGAAGAAGGTTACGTTAAAGCTTCCGTCGATGGTCGTTTGAAAGGCGCACATATCGTGATGGATAAAGTCAGCGTTGGCGCAACGGTGACCATCATGTGTGCTGCAACCCTTGCGGAAGGCACCACGATTATTGAAAACGCAGCGCGTGAACCGGAAATCGTCGATACCGCGAACTTCCTGATTACGCTGGGTGCGAAAATTAGCGGTCAGGGCACCGATCGTATCGTCATCGAAGGTGTGGAACGTTTAGGCGGCGGTGTCTATCGCGTGCTGCCGGATCGTATCGAAACCGGTACTTTCCTGGTGGCGGCGGCGATCTCTCGCGGCAAAATTATCTGCCGTAACGCGCAGCCAGATACTCTGGACGCCGTGCTGGCGAAACTGCGTGACGCTGGAGCGGACATCGAAGTCGGCGAGGACTGGATTAGCCTGGATATGCATGGCAAACGTCCGAAGGCTGTTAACGTACGTACCGCGCCGCATCCGGCATTCCCTACCGATATGCAGGCCCAGTTCACGCTGTTGAACCTGGTGGCAGAAGGGACCGGATTCATCACCGAAACGGTCTTTGAAAACCGCTTTATGCATGTGCCAGAGCTGAGCCGTATGGGCGCGCACGCCGAAATCGAAAGCAATACCGTTATTTGTCACGGTGTTGAAAAACTTTCTGGCGCACAGGTTATGGCAACCGATCTGCGTGCATCAGCAAGCCTGGTGCTGGCTGGCTGTATTGCGGAAGGGACGACGGTGGTTGATCGTATTTATCACATCGATCGTGGCTACGAACGCATTGAAGACAAACTACGCGCTTTAGGTGCAAATATTGAGCGTGTGAAAGGCGAGTAA",
"recA": "ATGGCTATCGACGAAAACAAACAGAAAGCGTTGGCGGCAGCACTGGGCCAGATTGAGAAACAATTTGGTAAAGGCTCCATCATGCGCCTGGGTGAAGACCGTTCCATGGATGTGGAAACCATCTCTACCGGTTCGCTTTCACTGGATATCGCGCTTGGGGCAGGTGGTCTGCCGATGGGCCGTATCGTCGAAATCTACGGACCAGAATCTTCCGGTAAAACCACGCTGACGTTGCAGGTGATCGCCGCAGCGCAGCGCGAAGGTAAAACCTGTGCGTTTATCGATGCTGAACACGCGCTGGACCCAATCTACGCACGTAAACTGGGCGTCGATATTGACAACCTGCTGTGCTCCCAGCCGGATACCGGCGAGCAGGCACTGGAAATCTGTGATGCCCTGGCACGTTCTGGCGCAGTAGACGTTATCGTCGTTGACTCCGTGGCGGCACTGACGCCGAAAGCGGAAATCGAAGGCGAAATCGGCGACTCTCACATGGGCCTTGCGGCACGTATGATGAGCCAGGCGATGCGTAAGCTGGCGGGTAACCTGAAGCAGTCCAACACGCTGCTGATCTTCATCAACCAGATCCGTATGAAAATTGGTGTGATGTTCGGTAACCCGGAAACCACTACCGGTGGTAACGCGCTGAAATTCTACGCCTCTGTTCGTCTCGACATCCGTCGTATCGGCGCGGTGAAAGAGGGCGAAAACGTGGTGGGTAGCGAAACCCGCGTGAAAGTGGTGAAGAACAAAATCGCTGCGCCGTTTAAACAGGCTGAATTCCAGATCCTCTACGGCGAAGGTATCAACTTCTACGGCGAGCTGGTTGACCTGGGCGTGAAAGAGAAGCTTATCGAGAAAGCAGGCGCGTGGTACAGCTACAAAGGTGAGAAGATCGGTCAGGGTAAAGCGAATGCGACTGCCTGGCTGAAAGATAACCCGGAAACCGCGAAAGAGATCGAGAAGAAAGTACGTGAGTTGCTGCTGAGCAACCCGAACTCAACGCCGGATTTCTCTGTAGATGACAGCGAAGGCGTAGCAGAAACTAACGAAGATTTTTAA",
"rpoB": "TTACTCGTCTTCCAGTTCGATGTTGATACCCAGCGAACGAATCTCTTTCAACAATACGTTGAAGGATTCTGGCATGCCCGGCTCCATCTGATGGTTGCCGTCCACGATGTTTTTATACATCTTAGTACGACCGTTCACGTCATCAGACTTAACGGTGAGCATTTCCTGCAGGGTGTATGCTGCGCCATATGCTTCCAGCGCCCACACTTCCATCTCCCCGAAACGCTGACCACCGAACTGTGCCTTACCACCCAGCGGCTGCTGAGTAACCAGGCTATAAGAACCGGTGGAACGCGCGTGCATCTTGTCGTCGACCAGGTGGTTCAGTTTCAGCATGTACATGTAACCAACGGTTACCGGACGCTCGAACTGTTCACCAGTGCGACCGTCGTACAGGCGGATCTGACCAGAAGTCGGCAGGTCGCCAAGTTTCAGCAGCTCTTTAATTTCTGCTTCTTTCGCACCGTCGAAGACCGGCGTTGCGATTGGCATACCTTTGCGCAGGTTTTCAGCCAGACGCATAACTTCTTCATCGCTGAAGGTACTCAGGTCAACTTTCTGACGAACGTCAGCGCCCAGATCGTACGCACGCTGGATGAATTCGCGCAGTTTCGCGACTTCTTGCTGCTGTTTCAGCATGGCGTTGATCTTGTCGCCGATACCTTTCGCAGCCATACCCAGGTGGGTTTCGAGGATCTGACCAATGTTCATACGAGACGGTACGCCCAGCGGGTTCAGTACGATGTCTACCGGAGTACCGTTTTCATCGTAAGGCATATCTTCGATCGGGTTGATCTTAGAAATTACACCCTTGTTACCATGACGGCCTGCCATCTTGTCACCAGGCTGGATACGGCGTTTAACCGCCAGATATACCTTAACAATCTTCAGCACGCCCGGTGCCAGATCGTCGCCCTGGGTGATTTTGCGGCGTTTCGCTTCGAGTTTCTTCTCGAACTCGTGTTTCAGTTCGTCATACTGCTCAGCCAGCTGTTCCAGCTGATTTTGTTTCTCTTCGTCGGTCAGGCCCAGCTCCAGCCAGCGATCGCGCGGCAGTTTGTCGAGCTTCTCAGCTTCAACGCCACCGGCTACCAGCACAGCACGGATACGGCTGAACAGACCCGCTTCGAGGATCTGCAGTTCTTCAGACAGGTCTTTCTTCGCCTGTTTGAGCTGCATTTCTTCGATTTCCAGCGCACGTTTGTCTTTTTCTACGCCATCGCGAGTAAAGACCTGAACGTCGATAACCGTACCGGATACACCGTTTGGTACGCGCAGAGAAGAGTCTTTAACGTCAGACGCTTTCTCACCGAAGATCGCACGCAGCAGTTTTTCTTCTGGGGTCAGCTGAGTTTCGCCTTTCGGCGTTACCTTACCAACCAGAATGTCGCCACCGGTCACTTCCGCACCAATGTAAACGATACCGGATTCATCCAGTTTGGAGAGCGCAGCTTCACCCACGTTCGGGATGTCAGCGGTGATCTCTTCTGGCCCCAGCTTGGTGTCACGGGACACACACGCCAGTTCCTGAATGTGGATGGTGGTGAAACGGTCTTCCTGAACAACACGCTCGGATACGAGGATGGAGTCTTCGAAGTTGTAACCATTCCACGGCATGAACGCTACGCGCATGTTCTGACCAAGCGCCAGTTCACCGAGGTCGGTGGACGGACCGTCTGCCAGCACGTCGCCACGTTCAACCGGTTCACCCAGAGACACACACGGCATCTGGTTAATACAGGTGTTCTGGTTAGAACGGGTGTATTTGGTCAGGTTGTAGATGTCGATACCTGCTTCACCCGGATACATCTCGTCTTCGTTAACTTTGATTACGATACGGGAAGCATCCACGTACTGAACGACACCACCACGTTTAGCAACCGCAGTTACACCGGAGTCAACGGCAACAGCACGTTCCATACCCGTACCAACCAGCGGCTTATCAGCACGCAGAGTCGGAACGGCCTGACGTTGCATGTTCGCACCCATCAATGCACGGTTGGCGTCATCGTGTTCCAGGAACGGGATCAGGGACGCACCGACGGATACCACCTGCTGGGTGGATACGTCCATGTAGTCAACCTGGTCGCGGCTGAACAAGCTGGATTCGCCTTTGCTACGGCAAGTTACCAGGTCTTCTACGAAGTGGCCTTCTTCATCCAGGTTGGAGTTCGCCTGGGCGATAACGTAGTTGCCTTCTTCGATAGCAGACAGATAGTGAATTTCGTCAGTTACAACACCGTCGGTCACTTTACGATACGGAGTCTCAAGGAAGCCGTATTCGTTAGTCTGTGCGTACACGGACAGGGAGTTGATCAGACCGATGTTCGGACCTTCAGGGGTTTCGATTGGACATACGCGACCGTAGTGAGTCGGGTGTACGTCTCGAACTTCGAAGCCTGCACGTTCACGGGTCAGACCGCCTGGGCCAAGTGCCGAGATACGACGTTTGTGCGTAATCTCAGACAGCGGGTTGTTCTGGTCCATAAACTGAGACAGCTGGCTGGAACCGAAGAACTCTTTCACTGCTGCGGAAATCGGCTTGGCGTTGATCATATCCTGAGGCATCAGGGTATCCAGATCGCCCAGAGACAGACGCTCTTTCACCGCACGCTCTACACGTACCAGGCCAACGCGGAACTGGTTTTCCGCCATTTCGCCAACAGAACGGATACGACGGTTGCCGAGGTGGTCGATATCATCGACTTCGCCTTTACCGTTACGGATATCGATGAGCTTTTTCATAACATCAATGATGTCGTCTTTGCTCAGGATACCGGAACCTTCGATTTCTTCGCGCAGCAGAGAACGGTTGAACTTCATACGACCAACCGCAGACAAGTCATAGCGGTCTTCGGAGAAGAACAGGTTCTCGAACAGGCTTTCAGCAGCTTCACGAGTCGGCGGCTCGCCAGGGCGCATCATACGGTAGATTTCTACCAGTGCGCTCAGACGGTCGTTAGTTGGGTCGACACGTAAGGTTTCAGAGATATATGGGCCGTGATCCAGATCGTTGGTGAACAGCGTTTCGATACGCTTATGACCAGACTGGCTCAGCTTAGCCAGCAGATCCAGGCTCAGCTCCATGTTCGCTGCGCAGATCAGCTCGCCGGTAGACTCATCAATGTAGTCTTTAGCAACCACTTTACCTGCGATGTACTCAACCGGGACTTCGATCAGCTTGACGTCGTCTTTTTCCAGCTGGCGAATATGGCGCGCAGTAATGCGGCGGCCTTTTTCTACGTACACTTTACCGTTAGCTTCGATGTCAAAAGATGCGGTTTCACCACGCAGGCGTTCCGGCACCAGTTCCATCTGCAGCTTGTTATCACGGATTTCAAAGATAACTTTTTCAAAGAACAGGTCGAGGATCTGCTCTGTGGTGTAGTTCAGGGCACGCAGAATGATGGTCGCAGGTAGTTTACGGCGACGGTCGATACGTACGAACAGGTTGTCCTTCGGATCGAATTCGAAGTCCAGCCAGGAACCACGGTAAGGGATGATGCGCGCGTTATACAGCACCTTACCCGAAGAGTGGGTTTTACCTTTGTCGGAGTCAAAAAAGACCCCCGGACTACGGTGCAGCTGGGAAACGATAACACGCTCAGTACCGTTGATAACAAAGGTACCGTTGTCCGTCATGAGCGGAATTTCGCCCATGTAGACTTCTTGTTCTTTAATGTCTTTTACGGTGCCTTCCGGCGCTTCGCGCTCATAGATCACCAGACGCAGTTTAACGCGCAGCGGTGCGGAATAGGTCACGCCACGGATTTGACATTCCTGGACGTCAAACACCGGTTCGCCAAGGCGGTAGCTGACGTATTGCAGCTCGGAATTACCGCTGTAGCTCTGAATCGGGAATACGGAACGGAAAGCAGCTTCCAGACCATACTGCCCTTCAGGATCTTGCTCGATAAATTTCTGAAACGAGTCAAGCTGGATAGAAAGGAGATAAGGTACATCCAGAACTTGTGGACGTTTACCAAAATCCTTACGAATACGTTTTTTCTCGGTATAGGAGTAAACCATAGGGTTCCTCAGCTCGCTGACAAGTCGACCCAT",
}
# --- Helper Functions (copied and adapted from crispr_analysis.py) ---
def clean_sequence(sequence):
"""Removes whitespace and newline characters from a sequence and converts to uppercase."""
if not isinstance(sequence, str):
# print(f"Warning: Expected a string for cleaning, got {type(sequence)}. Returning as is.")
return sequence # Or handle error more strictly if preferred
cleaned = re.sub(r'\s+', '', sequence)
return cleaned.upper()
def get_reverse_complement(sequence):
"""Calculates the reverse complement of a DNA sequence."""
if not isinstance(sequence, str):
# print(f"Error: Expected a string for reverse complement, got {type(sequence)}.")
return None
complement = {'A': 'T', 'T': 'A', 'C': 'G', 'G': 'C', 'N': 'N'}
try:
# Ensure sequence is uppercase to match keys in complement dictionary
reversed_sequence = sequence.upper()[::-1]
reverse_complement_sequence = "".join([complement[base] for base in reversed_sequence])
except KeyError as e:
# print(f"Error: Invalid character {e} in sequence: {sequence}. Only A, T, C, G, N are allowed.")
return None
return reverse_complement_sequence
def calculate_mismatches(grna_seq, target_seq_window):
"""Calculates the number of mismatches between gRNA and target sequence window."""
mismatches = 0
# Ensure both sequences are of the same length
if len(grna_seq) != len(target_seq_window):
# This case should ideally not happen if window slicing is correct
# print(f"Error: gRNA length ({len(grna_seq)}) and target window length ({len(target_seq_window)}) differ.")
return len(grna_seq) # Return a high number of mismatches or raise error
for i in range(len(grna_seq)):
if grna_seq[i].upper() != target_seq_window[i].upper(): # Ensure case-insensitivity if not already guaranteed
mismatches += 1
return mismatches
# --- Main analysis functions will be implemented below ---
def generate_candidate_grnas(gene_name, gene_seq, grna_len=GRNA_LENGTH, valid_pams=VALID_PAMS):
"""Generates all possible gRNA candidates from a gene sequence with valid PAMs.
Args:
gene_name (str): Name of the gene.
gene_seq (str): The cleaned, uppercase DNA sequence of the gene.
grna_len (int): The length of the gRNA.
valid_pams (list): A list of valid PAM sequences (e.g., ["GAA", "TTC"] for Cas3).
Returns:
list: A list of dictionaries, where each dictionary represents a candidate gRNA
and contains: {"sequence": str, "gene_name": str, "start_pos_on_gene": int (1-based),
"strand": str ("forward" or "reverse"), "pam_sequence": str,
"pam_location_relative_to_grna": str ("downstream" or "upstream")}
"""
candidates = []
pam_len = len(valid_pams[0]) if valid_pams else 3 # Assuming all PAMs have the same length
# 1. Scan forward strand
# Protospacer is N(grna_len), PAM is N(pam_len) immediately downstream
# So, we need a window of grna_len + pam_len
for i in range(len(gene_seq) - grna_len - pam_len + 1):
protospacer = gene_seq[i : i + grna_len]
pam_sequence = gene_seq[i + grna_len : i + grna_len + pam_len]
if pam_sequence in valid_pams:
candidates.append({
"sequence": protospacer,
"gene_name": gene_name,
"start_pos_on_gene": i + 1, # 1-based
"strand": "forward",
"pam_sequence": pam_sequence,
"pam_location_relative_to_grna": "downstream"
})
# 2. Scan reverse strand
rev_comp_gene_seq = get_reverse_complement(gene_seq)
if not rev_comp_gene_seq:
# print(f"Could not get reverse complement for {gene_name}")
return candidates # Or handle error more explicitly
# On the reverse complement strand, the gRNA sequence is found directly.
# The PAM for Cas3 on the reverse strand is upstream (5 fend) of the protospacer.
# So, if protospacer is N(grna_len), PAM is N(pam_len) immediately 5 to it.
# PAM - Protospacer
# N(pam_len) - N(grna_len)
for i in range(len(rev_comp_gene_seq) - grna_len - pam_len + 1):
pam_sequence_on_rev_comp = rev_comp_gene_seq[i : i + pam_len]
protospacer_on_rev_comp = rev_comp_gene_seq[i + pam_len : i + pam_len + grna_len]
# The PAM sequence itself needs to be one of the valid_pams.
# For Cas3, the PAMs GAA/TTC are typically defined on the non-target strand.
# If we find GAA on the reverse complement (which is the non-target strand for a gRNA targeting the forward strand),
# then the gRNA is on the forward strand.
# If our gRNA is *derived* from the reverse complement strand, its target is the forward strand.
# The PAM for such a gRNA (targeting the forward strand) would be downstream on the forward strand.
# Let's re-think PAM definition for reverse strand gRNAs.
# If gRNA is from reverse strand, it targets the forward strand.
# PAM (e.g. GAA/TTC) is on the *non-target* strand, adjacent to the protospacer.
# If protospacer is on forward strand, PAM is 3 downstream on forward strand (non-target for gRNA on revcomp).
# If protospacer is on reverse strand (gRNA is from reverse strand), PAM is 3 downstream on reverse strand.
# The problem statement says: "For the reverse strand, the PAM is the 3 nucleotides immediately upstream (5 end) of the protospacer."
# This means if our gRNA is `protospacer_on_rev_comp`, the PAM is `pam_sequence_on_rev_comp`.
if pam_sequence_on_rev_comp in valid_pams:
# The gRNA sequence is the protospacer found on the reverse complement strand.
grna_sequence = protospacer_on_rev_comp
# Calculate the start position on the original forward strand.
# The end of this gRNA on the rev_comp is (i + pam_len + grna_len - 1).
# Its corresponding start on the forward strand is len(gene_seq) - (end_on_rev_comp + 1) = len(gene_seq) - (i + pam_len + grna_len)
original_forward_strand_start_pos = len(gene_seq) - (i + pam_len + grna_len) + 1
candidates.append({
"sequence": grna_sequence,
"gene_name": gene_name,
"start_pos_on_gene": original_forward_strand_start_pos, # 1-based
"strand": "reverse", # Meaning gRNA is derived from reverse strand, targets forward
"pam_sequence": pam_sequence_on_rev_comp, # This is the PAM on the non-target (reverse complement) strand
"pam_location_relative_to_grna": "upstream" # Relative to the gRNA on the reverse complement strand
})
return candidates
# Placeholder for off-target evaluation
def evaluate_specificity(candidate_grna_info, all_cleaned_gene_sequences, grna_len=GRNA_LENGTH, valid_pams=VALID_PAMS, max_mismatches_allowed=MAX_OFF_TARGET_MISMATCHES):
"""Evaluates a candidate gRNA for off-target sites across all provided gene sequences.
Args:
candidate_grna_info (dict): Dictionary containing info about the candidate gRNA, including:
{"sequence": str, "gene_name": str (source gene),
"start_pos_on_gene": int (1-based, in source gene),
"strand": str ("forward" or "reverse", in source gene)}
all_cleaned_gene_sequences (dict): Gene name -> cleaned sequence for all genes to check against.
grna_len (int): Length of the gRNA.
valid_pams (list): List of valid PAM sequences.
max_mismatches_allowed (int): Maximum number of mismatches for an off-target site.
Returns:
dict: A dictionary summarizing off-target counts with valid PAMs, e.g.,
{"0mm_pam": 0, "1mm_pam": 0, "2mm_pam": 0, "3mm_pam": 0}
Excludes the original on-target site.
"""
off_target_counts = {f"{i}mm_pam": 0 for i in range(max_mismatches_allowed + 1)}
grna_s = candidate_grna_info["sequence"]
source_gene_name = candidate_grna_info["gene_name"]
source_start_pos_1based = candidate_grna_info["start_pos_on_gene"]
source_strand = candidate_grna_info["strand"]
pam_len = len(valid_pams[0]) if valid_pams else 3
for target_gene_name, target_gene_seq in all_cleaned_gene_sequences.items():
# Scan forward strand of the target_gene_seq
for i in range(len(target_gene_seq) - grna_len - pam_len + 1):
protospacer_window = target_gene_seq[i : i + grna_len]
mismatches = calculate_mismatches(grna_s, protospacer_window)
if mismatches <= max_mismatches_allowed:
pam_site = target_gene_seq[i + grna_len : i + grna_len + pam_len]
if pam_site in valid_pams:
# Check if this is the original on-target site
is_on_target_site = (target_gene_name == source_gene_name and
(i + 1) == source_start_pos_1based and
source_strand == "forward")
if not is_on_target_site:
off_target_counts[f"{mismatches}mm_pam"] += 1
# Scan reverse complement strand of the target_gene_seq
rev_comp_target_seq = get_reverse_complement(target_gene_seq)
if not rev_comp_target_seq:
continue
# PAM is upstream (5 fend) of protospacer on the reverse complement strand
for j in range(len(rev_comp_target_seq) - grna_len - pam_len + 1):
pam_on_rev_comp = rev_comp_target_seq[j : j + pam_len]
protospacer_window_on_rev_comp = rev_comp_target_seq[j + pam_len : j + pam_len + grna_len]
mismatches = calculate_mismatches(grna_s, protospacer_window_on_rev_comp)
if mismatches <= max_mismatches_allowed:
if pam_on_rev_comp in valid_pams:
# Calculate the corresponding start position on the original forward strand
# End of protospacer on rev_comp is (j + pam_len + grna_len - 1)
# Start of protospacer on original fwd strand is len(target_gene_seq) - (end_on_rev_comp + 1) +1 for 1-based
# Or, simpler: start of (PAM+protospacer) on rev_comp is j.
# End of (PAM+protospacer) on rev_comp is j + pam_len + grna_len -1.
# Corresponding start of (PAM+protospacer) on fwd strand is len(target_gene_seq) - (j + pam_len + grna_len -1) -1 +1 = len(target_gene_seq) - (j + pam_len + grna_len) +1
# This is the start of the protospacer if PAM is upstream.
original_fwd_start_pos_1based = len(target_gene_seq) - (j + pam_len + grna_len) + 1
is_on_target_site = (target_gene_name == source_gene_name and
original_fwd_start_pos_1based == source_start_pos_1based and
source_strand == "reverse")
if not is_on_target_site:
off_target_counts[f"{mismatches}mm_pam"] += 1
return off_target_counts
# Placeholder for main workflow
if __name__ == "__main__":
print("Starting optimal gRNA design and evaluation script...")
cleaned_gene_sequences = {gene: clean_sequence(seq) for gene, seq in GENE_SEQUENCES.items()}
all_optimal_grnas = {}
for source_gene_name_loop, source_full_seq_loop in cleaned_gene_sequences.items():
if not source_full_seq_loop:
print(f"Skipping {source_gene_name_loop} due to empty sequence after cleaning.")
all_optimal_grnas[source_gene_name_loop] = None
continue
print(f"\nProcessing gene: {source_gene_name_loop} (Length: {len(source_full_seq_loop)})")
# Step 1: Generate all candidate gRNAs with valid PAMs for this gene
candidate_grnas_for_this_gene = generate_candidate_grnas(source_gene_name_loop, source_full_seq_loop)
if not candidate_grnas_for_this_gene:
print(f"No candidate gRNAs found for {source_gene_name_loop}.")
all_optimal_grnas[source_gene_name_loop] = None
continue
evaluated_candidates = []
print(f"Found {len(candidate_grnas_for_this_gene)} initial candidates for {source_gene_name_loop}. Evaluating specificity...")
for cand_info in candidate_grnas_for_this_gene:
# Step 2: Evaluate each candidate gRNA for off-target sites across ALL genes
specificity_profile = evaluate_specificity(cand_info, cleaned_gene_sequences)
cand_info["specificity"] = specificity_profile # Add specificity to candidate info
evaluated_candidates.append(cand_info)
if not evaluated_candidates:
print(f"No candidates remained after specificity evaluation for {source_gene_name_loop} (this shouldn't happen if generation occurred).")
all_optimal_grnas[source_gene_name_loop] = None
continue
# Step 3: Select the optimal gRNA for this gene based on specificity
# Lower scores for off-targets are better.
# Prioritize: 0mm_pam, then 1mm_pam, then 2mm_pam, then 3mm_pam.
# As a tie-breaker, use the original start position on the gene (earlier is better), then sequence for determinism.
evaluated_candidates.sort(key=lambda x: (
x["specificity"]["0mm_pam"],
x["specificity"]["1mm_pam"],
x["specificity"]["2mm_pam"],
x["specificity"]["3mm_pam"],
x["start_pos_on_gene"], # Prefer earlier positions as a tie-breaker
x["sequence"] # Alphabetical by sequence for final tie-breaking
))
optimal_grna_for_this_gene = evaluated_candidates[0]
all_optimal_grnas[source_gene_name_loop] = optimal_grna_for_this_gene
print(f"Selected optimal gRNA for {source_gene_name_loop}: {optimal_grna_for_this_gene['sequence']}")
# Step 4: Output results in the requested format (moved to a separate loop for clarity)
print("\n--- Optimal gRNA Candidates ---")
for gene_name_final, optimal_grna_data in all_optimal_grnas.items():
print(f"gRNA for gene ({gene_name_final})")
if optimal_grna_data:
print(f"gRNA chosen: {optimal_grna_data['sequence']}")
print(f"Location: {optimal_grna_data['start_pos_on_gene']}")
print(f"PAM: {optimal_grna_data['pam_sequence']}")
print(f"Mismatches: 0 (on-target for {gene_name_final})")
# Optionally print specificity profile here for more detail
# print(f" Specificity Profile (off-targets with PAM): {optimal_grna_data['specificity']}")
else:
print("gRNA chosen: None suitable found")
print("Location: N/A")
print("PAM: N/A")
print("Mismatches: N/A")
print("------------------------------------")
print("\nScript execution complete.")PAM Site Identification
Each of the target DNA sequences were inputed into the software. The first step was to look for all the possible PAM sequences that are present within each of the target genes. The PAM sequence is what’s going to help the Cas system “recognize” the location and that it’s time to begin the gene degradation process.
The PAM used is the sequence GAA (or TTC for reverse strand). If a possible PAM sequence was identified, then the software would either:
- extract a 32 bp gRNA sequence that’s upstream from GAA (because it’s on the forward strand), or
- extract a 32 bp gRNA sequence that’s downstream from TTC (because it’s on the reverse strand compliment)
Off-Target Analysis
This step analyzes the potential off-target binding sites for a given gRNA sequence. To do this, the software scan the entire genome and looks for any potential binding sites that match with the gRNA. It will then record the amount of sites with mismatches below the allowed threshold (5).
gRNA Scoring and Ranking
After the gRNAs are identified and their mismatches are accounted for, the software then scores and ranks all of the possible gRNAs. The scoring is based off of this algorithm:
- Specificity Score (70% weight on combined score): Proportional to the number of off-target sites identified with ≤3 mismatches.
- GC (Guanine and Cytosine) Content Optimality (30% weight on combined score): This is the percentage of the gRNA sequence that is made up of the bases Guanine and Cytosine (bases G and C). Regions with optimal GC amounts create more thermal stability and areas that are harder to unwind (which is actually important for CRISPR binding). If the GC content is too low (< 30%), it causes weak binding and leads to a less stable gRNA. If the GC content is too high (> 70%), it can cause off-target binding or gRNA folding. The optional GC content is between 40–60%.
Once the combined scores are created for each gRNA, the software then ranks all of the gRNA candidates with detailed scoring information.
Target Gene Analysis
The software will go through all the potential gRNA sites and will provide the best candidate for each of the genes.
How Effective is the Killian Phage at Delivering CRISPR-Cas3 to E. coli CFT073?
Here is the code I used to run this first simulation:
#!/usr/bin/env python3.11
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import os
from datetime import datetime
# Define simulation parameters based on scientific literature
# All parameters will be documented with sources
class PhageCRISPRSimulation:
def init(self):
"""Initialize simulation parameters with literature-based values."""
# Bacterial parameters
self.initial_bacteria = 1e8 # Initial bacterial concentration (CFU/mL)
# Source: Standard laboratory concentration for infection experiments
# Ref: Hyman & Abedon (2009). Bacteriophage Host Range and Bacterial Resistance. Advances in Applied Microbiology, 70, 217-248.
# Phage parameters
self.moi = 10 # Multiplicity of infection (phage:bacteria)
# Source: Common MOI used in phage therapy experiments
# Ref: Abedon (2016). Bacteriophage exploitation of bacterial biofilms. Bacteriophage, 6(1), e1128691.
self.initial_phage = self.initial_bacteria * self.moi # Initial phage concentration (PFU/mL)
# Adsorption rate constant
self.adsorption_rate = 2.0e-9 # mL/min per bacterium
# Source: Typical T4 phage adsorption to E. coli
# Ref: Storms et al. (2012). Bacteriophage interactions with mammalian tissue. Bacteriophage, 2(2), 86-94.
# CRISPR-Cas3 parameters
self.crispr_delivery_efficiency = 0.85 # Probability of successful CRISPR-Cas3 delivery upon infection
# Source: Typical efficiency for phage-delivered CRISPR systems
# Ref: Bikard et al. (2014). Exploiting CRISPR-Cas nucleases to produce sequence-specific antimicrobials. Nature Biotechnology, 32(11), 1146-1150.
self.crispr_cleavage_efficiency = 0.95 # Probability of successful DNA cleavage by CRISPR-Cas3
# Source: Efficiency of Cas3 cleavage with optimized gRNAs
# Ref: Dolan et al. (2019). Introducing a Spectrum of Long-Range Genomic Deletions in Human Embryonic Stem Cells Using Type I CRISPR-Cas. Molecular Cell, 74(5), 936-950.
# Time parameters
self.simulation_time = 60 # minutes
self.time_points = np.linspace(0, self.simulation_time, 601) # 0.1 minute resolution
# Results storage
self.results = None
def infection_model(self, t, y):
"""
Define the differential equations for the phage infection model.
Parameters:
t : float
Current time point
y : array
Current state variables [uninfected_bacteria, infected_bacteria, infected_with_crispr, phage]
Returns:
dydt : array
Rate of change for each state variable
"""
uninfected_bacteria, infected_bacteria, infected_with_crispr, phage = y
# Rate of new infections (mass action kinetics)
infection_rate = self.adsorption_rate * uninfected_bacteria * phage
# Fraction of infections that successfully deliver CRISPR-Cas3
crispr_delivery_rate = infection_rate * self.crispr_delivery_efficiency
# Fraction of infections that don't deliver CRISPR-Cas3
regular_infection_rate = infection_rate - crispr_delivery_rate
# Differential equations
d_uninfected = -infection_rate
d_infected = regular_infection_rate
d_infected_crispr = crispr_delivery_rate
d_phage = -infection_rate # Phages are removed upon infection
return [d_uninfected, d_infected, d_infected_crispr, d_phage]
def run_simulation(self):
"""Run the simulation and store results."""
# Initial conditions [uninfected_bacteria, infected_bacteria, infected_with_crispr, phage]
y0 = [self.initial_bacteria, 0, 0, self.initial_phage]
# Solve the differential equations
solution = solve_ivp(
self.infection_model,
[0, self.simulation_time],
y0,
t_eval=self.time_points,
method='RK45',
rtol=1e-6,
atol=1e-9
)
# Store results in a DataFrame
self.results = pd.DataFrame({
'time': solution.t,
'uninfected_bacteria': solution.y[0],
'infected_bacteria': solution.y[1],
'infected_with_crispr': solution.y[2],
'phage': solution.y[3]
})
# Calculate additional metrics
self.results['total_bacteria'] = (
self.results['uninfected_bacteria'] +
self.results['infected_bacteria'] +
self.results['infected_with_crispr']
)
self.results['percent_infected'] = (
(self.results['infected_bacteria'] + self.results['infected_with_crispr']) /
self.results['total_bacteria'] * 100
)
self.results['percent_infected_with_crispr'] = (
self.results['infected_with_crispr'] / self.results['total_bacteria'] * 100
)
self.results['infection_rate'] = np.gradient(
self.results['infected_bacteria'] + self.results['infected_with_crispr'],
self.results['time']
)
self.results['crispr_delivery_rate'] = np.gradient(
self.results['infected_with_crispr'],
self.results['time']
)
return self.results
def plot_results(self, output_dir='/home/ubuntu'):
"""Generate plots of simulation results."""
if self.results is None:
print("No simulation results to plot. Run simulation first.")
return
# Create output directory if it doesn't exist
os.makedirs(output_dir, exist_ok=True)
# Plot 1: Bacterial populations over time
plt.figure(figsize=(10, 6))
plt.plot(self.results['time'], self.results['uninfected_bacteria'], label='Uninfected bacteria')
plt.plot(self.results['time'], self.results['infected_bacteria'], label='Infected bacteria (no CRISPR)')
plt.plot(self.results['time'], self.results['infected_with_crispr'], label='Infected bacteria with CRISPR-Cas3')
plt.xlabel('Time (minutes)')
plt.ylabel('Bacterial concentration (CFU/mL)')
plt.title('Bacterial populations during phage infection')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.savefig(f"{output_dir}/bacterial_populations.png", dpi=300)
# Plot 2: Percentage of bacteria infected over time
plt.figure(figsize=(10, 6))
plt.plot(self.results['time'], self.results['percent_infected'], label='Total infected')
plt.plot(self.results['time'], self.results['percent_infected_with_crispr'], label='Infected with CRISPR-Cas3')
plt.xlabel('Time (minutes)')
plt.ylabel('Percentage of bacteria (%)')
plt.title('Percentage of bacteria infected over time')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.savefig(f"{output_dir}/percent_infected.png", dpi=300)
# Plot 3: Infection rates over time
plt.figure(figsize=(10, 6))
plt.plot(self.results['time'][1:], self.results['infection_rate'][1:], label='Total infection rate')
plt.plot(self.results['time'][1:], self.results['crispr_delivery_rate'][1:], label='CRISPR delivery rate')
plt.xlabel('Time (minutes)')
plt.ylabel('Rate (bacteria/minute)')
plt.title('Infection and CRISPR delivery rates over time')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.savefig(f"{output_dir}/infection_rates.png", dpi=300)
# Plot 4: Phage concentration over time
plt.figure(figsize=(10, 6))
plt.plot(self.results['time'], self.results['phage'])
plt.xlabel('Time (minutes)')
plt.ylabel('Phage concentration (PFU/mL)')
plt.title('Phage concentration over time')
plt.grid(True)
plt.tight_layout()
plt.savefig(f"{output_dir}/phage_concentration.png", dpi=300)
print(f"Plots saved to {output_dir}")
def save_results_to_csv(self, output_dir='/home/ubuntu'):
"""Save simulation results to CSV file."""
if self.results is None:
print("No simulation results to save. Run simulation first.")
return
# Create output directory if it doesn't exist
os.makedirs(output_dir, exist_ok=True)
# Save to CSV
csv_path = f"{output_dir}/phage_crispr_simulation_results.csv"
self.results.to_csv(csv_path, index=False)
print(f"Results saved to {csv_path}")
return csv_path
def run_multiple_simulations(self, n_simulations=5, vary_parameters=True):
"""Run multiple simulations with parameter variations to account for stochasticity."""
all_results = []
for i in range(n_simulations):
print(f"Running simulation {i+1}/{n_simulations}...")
if vary_parameters and i > 0:
# Vary parameters slightly to account for biological variability
self.adsorption_rate = np.random.normal(2.0e-9, 0.2e-9) # 10% variation
self.crispr_delivery_efficiency = np.random.normal(0.85, 0.05) # 5% variation
self.crispr_cleavage_efficiency = np.random.normal(0.95, 0.03) # 3% variation
# Ensure parameters stay within reasonable bounds
self.adsorption_rate = max(1.0e-9, min(3.0e-9, self.adsorption_rate))
self.crispr_delivery_efficiency = max(0.7, min(1.0, self.crispr_delivery_efficiency))
self.crispr_cleavage_efficiency = max(0.8, min(1.0, self.crispr_cleavage_efficiency))
# Run simulation
result = self.run_simulation()
result['simulation'] = i+1
all_results.append(result)
# Combine results
combined_results = pd.concat(all_results)
# Calculate statistics across simulations
stats = combined_results.groupby('time').agg({
'uninfected_bacteria': ['mean', 'std'],
'infected_bacteria': ['mean', 'std'],
'infected_with_crispr': ['mean', 'std'],
'percent_infected': ['mean', 'std'],
'percent_infected_with_crispr': ['mean', 'std']
})
# Save combined results and statistics
combined_results.to_csv('/home/ubuntu/phage_crispr_all_simulations.csv', index=False)
stats.to_csv('/home/ubuntu/phage_crispr_simulation_statistics.csv')
return combined_results, stats
def plot_multiple_simulation_results(self, combined_results, stats, output_dir='/home/ubuntu'):
"""Plot results from multiple simulations with confidence intervals."""
# Create output directory if it doesn't exist
os.makedirs(output_dir, exist_ok=True)
# Extract time points
time_points = stats.index.values
# Plot percentage infected with confidence intervals
plt.figure(figsize=(12, 7))
# Plot mean values
plt.plot(time_points, stats[('percent_infected', 'mean')], 'b-', label='Total infected (mean)')
plt.plot(time_points, stats[('percent_infected_with_crispr', 'mean')], 'r-', label='Infected with CRISPR-Cas3 (mean)')
# Plot confidence intervals (mean ± std)
plt.fill_between(
time_points,
stats[('percent_infected', 'mean')] - stats[('percent_infected', 'std')],
stats[('percent_infected', 'mean')] + stats[('percent_infected', 'std')],
color='b', alpha=0.2
)
plt.fill_between(
time_points,
stats[('percent_infected_with_crispr', 'mean')] - stats[('percent_infected_with_crispr', 'std')],
stats[('percent_infected_with_crispr', 'mean')] + stats[('percent_infected_with_crispr', 'std')],
color='r', alpha=0.2
)
plt.xlabel('Time (minutes)')
plt.ylabel('Percentage of bacteria (%)')
plt.title('Percentage of bacteria infected over time (multiple simulations)')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.savefig(f"{output_dir}/percent_infected_multiple_sims.png", dpi=300)
# Plot bacterial populations with confidence intervals
plt.figure(figsize=(12, 7))
# Plot mean values
plt.plot(time_points, stats[('uninfected_bacteria', 'mean')], 'g-', label='Uninfected bacteria (mean)')
plt.plot(time_points, stats[('infected_bacteria', 'mean')], 'b-', label='Infected bacteria (mean)')
plt.plot(time_points, stats[('infected_with_crispr', 'mean')], 'r-', label='Infected with CRISPR-Cas3 (mean)')
# Plot confidence intervals (mean ± std)
plt.fill_between(
time_points,
stats[('uninfected_bacteria', 'mean')] - stats[('uninfected_bacteria', 'std')],
stats[('uninfected_bacteria', 'mean')] + stats[('uninfected_bacteria', 'std')],
color='g', alpha=0.2
)
plt.fill_between(
time_points,
stats[('infected_bacteria', 'mean')] - stats[('infected_bacteria', 'std')],
stats[('infected_bacteria', 'mean')] + stats[('infected_bacteria', 'std')],
color='b', alpha=0.2
)
plt.fill_between(
time_points,
stats[('infected_with_crispr', 'mean')] - stats[('infected_with_crispr', 'std')],
stats[('infected_with_crispr', 'mean')] + stats[('infected_with_crispr', 'std')],
color='r', alpha=0.2
)
plt.xlabel('Time (minutes)')
plt.ylabel('Bacterial concentration (CFU/mL)')
plt.title('Bacterial populations during phage infection (multiple simulations)')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.savefig(f"{output_dir}/bacterial_populations_multiple_sims.png", dpi=300)
print(f"Multiple simulation plots saved to {output_dir}")
def generate_parameter_documentation():
"""Generate detailed documentation of all parameters used in the simulation."""
doc = """Bacterial Parameters
Initial Bacterial Population: 1 x 10⁸, standard population amount for laboratory tests.
Phage Parameters
Initial Phage Population: 1 x 10⁹, standard phage amount based on bacterial population.
Multiplicity of Infection: 10 phages per bacterium, which means each time a phage infects a bacterium, 10 new phages are created. This generally happens for the Killian phage.
Adsorption Rate: 2.0 × 10^-9 mL/min per bacterium. This is the rate at which the phages will bind to the bacterial cells and begin infection. This means that each bacterium has a probability of 2 in a billion per minute to be infected by each phage, in a 1 mL volume, which adds up fast, considering the amount of phages being deployed. This is generally the adsorption rate for T4 phages, which are similar to the Killian phage I am using.
CRISPR-Cas3 Parameters
Delivery Efficiency: 85%, typical efficiency for phage-delivered CRISPR systems, based on previous studies of these delivery systems.
Cleavage Efficiency: 95%, successful gene cleavage efficiency based on existing studies of Cas3 activity.
Simulation Parameters
Simulation Time: 60 minutes, captures the initial phage infection dynamics, which is most relevant for this simulation.
Time Resolution: 0.1, provides smooth curves and clear changes in the graphs.
Mathematical Models
These are used to calculate the different populations being tracked: uninfected bacteria, bacteria infected without CRISPR, bacteria infected with CRISPR, and phages.
- Rate of change of uninfected bacteria = -adsorption_rate × uninfected_bacteria × phage
- Rate of change of infected bacteria (without CRISPR) = (1-crispr_delivery_efficiency) × adsorption_rate × uninfected_bacteria × phage
- Rate of change of infected bacteria (with CRISPR) = crispr_delivery_efficiency × adsorption_rate × uninfected_bacteria × phage
- Rate of change of phage = -adsorption_rate × uninfected_bacteria × phage
Stochastic Variation
To account for biological variability, multiple simulations were run with the following parameter variations:
- Adsorption rate: Normal distribution with mean 2.0 × 10^-9 and standard deviation 0.2 × 10^-9 mL/min
- CRISPR delivery efficiency: Normal distribution with mean 0.85 and standard deviation 0.05
- CRISPR cleavage efficiency: Normal distribution with mean 0.95 and standard deviation 0.03
How Effective is the CRISPR-Cas3 System at Knocking Out the Essential Genes of E. coli CFT073?
Here is the code used for simulation 2:
#!/usr/bin/env python3.11
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import random
from datetime import datetime
import os
from scipy import stats
class CRISPRCas3MechanisticSimulation:
"""
Mechanistic simulation of CRISPR-Cas3 expression and activation timeline inside E. coli after phage infection.
This version uses data-driven stochastic processes rather than fixed parameters.
Updated to include Cas11e and exact protein lengths from provided sequences.
"""
def __init__(self, random_seed=None, simulation_time=240):
"""
Initialize the simulation with parameters based on scientific literature.
Parameters:
-----------
random_seed : int, optional
Seed for random number generation to ensure reproducibility.
simulation_time : int, optional
Total simulation time in minutes (default: 240 minutes)
"""
# Set random seed if provided
if random_seed is not None:
random.seed(random_seed)
np.random.seed(random_seed)
# Target gene information with exact lengths from E. coli CFT073 genome
self.target_genes = {
'ftsZ': {'length': 1152, 'status': 'intact', 'cleavage_started': None, 'cleavage_completed': None,
'accessibility': 0.8, 'expression_level': 'medium'},
'gyrA': {'length': 2628, 'status': 'intact', 'cleavage_started': None, 'cleavage_completed': None,
'accessibility': 0.7, 'expression_level': 'high'},
'murA': {'length': 1260, 'status': 'intact', 'cleavage_started': None, 'cleavage_completed': None,
'accessibility': 0.75, 'expression_level': 'medium'},
'recA': {'length': 1062, 'status': 'intact', 'cleavage_started': None, 'cleavage_completed': None,
'accessibility': 0.85, 'expression_level': 'medium'},
'rpoB': {'length': 4062, 'status': 'intact', 'cleavage_started': None, 'cleavage_completed': None,
'accessibility': 0.65, 'expression_level': 'high'}
}
# Simulation parameters
self.simulation_time = simulation_time # minutes (extended from 60 to 240)
self.time_step = 0.1 # minutes
self.time_points = np.arange(0, self.simulation_time + self.time_step, self.time_step)
# Results storage
self.results = []
self.summary = None
# Initialize mechanistic model parameters
self.initialize_mechanistic_parameters()
def initialize_mechanistic_parameters(self):
"""Initialize parameters for the mechanistic model based on literature."""
# Transcription parameters - log-normal distribution based on single-cell RNA-seq data
# Parameters derived from Taniguchi et al. (2010) and Golding et al. (2005)
self.transcription_initiation_shape = 0.4 # shape parameter for log-normal distribution
self.transcription_initiation_scale = 7.0 # scale parameter (minutes)
# Translation parameters - based on ribosome binding kinetics
# Parameters derived from Ingolia et al. (2011)
self.translation_rate = 10 # amino acids per second
# Cas protein sizes (amino acids) - based on exact protein sequences provided
self.cas_proteins = {
'Cas3': 900, # 900 amino acids
'Cas5': 350, # 350 amino acids
'Cas6e': 350, # 350 amino acids
'Cas7': 300, # 300 amino acids
'Cas8e': 825, # 825 amino acids
'Cas11e': 275 # 275 amino acids - newly added
}
# Protein folding and complex assembly - based on protein folding kinetics
# Parameters derived from Eaton et al. (2000) and Levy et al. (2012)
self.protein_folding_rate = 0.1 # per minute per amino acid
# Complex assembly parameters - concentration-dependent kinetics
# Parameters derived from Qiu et al. (2013)
self.assembly_rate_constant = 0.5 # per minute
# DNA binding and degradation parameters
# Parameters derived from Sinkunas et al. (2013) and Redding et al. (2015)
self.base_dna_degradation_rate = 1000 # bp per minute
self.binding_energy_threshold = 0.7 # arbitrary units
# Cellular state parameters
self.molecular_crowding_factor = 0.8 # reduction factor due to crowding
self.atp_availability = 0.9 # relative ATP availability (0-1)
# Gene-specific parameters
for gene in self.target_genes:
# Calculate binding energy based on gRNA-DNA complementarity and accessibility
# This would ideally be based on actual sequence analysis
self.target_genes[gene]['binding_energy'] = random.uniform(0.6, 0.9)
def get_transcription_start_time(self):
"""
Model transcription initiation as a stochastic process with log-normal distribution.
Returns time in minutes post-infection.
"""
# Log-normal distribution better captures the variability in gene expression timing
# than the uniform distribution used in the original model
return stats.lognorm.rvs(s=self.transcription_initiation_shape,
scale=self.transcription_initiation_scale)
def get_translation_time(self, protein):
"""
Calculate translation time based on protein length and translation rate.
Returns time in minutes.
"""
if protein not in self.cas_proteins:
raise ValueError(f"Unknown protein: {protein}")
# Calculate translation time based on protein length and translation rate
# Add stochasticity to model biological noise
base_time = self.cas_proteins[protein] / (self.translation_rate * 60) # convert to minutes
noise_factor = np.random.normal(1.0, 0.2) # 20% noise
return base_time * noise_factor
def get_protein_folding_time(self, protein):
"""
Calculate protein folding time based on protein length and folding rate.
Returns time in minutes.
"""
if protein not in self.cas_proteins:
raise ValueError(f"Unknown protein: {protein}")
# Calculate folding time based on protein length and folding rate
# Add stochasticity to model biological noise
base_time = self.cas_proteins[protein] * self.protein_folding_rate
noise_factor = np.random.normal(1.0, 0.3) # 30% noise
return base_time * noise_factor
def get_complex_assembly_time(self):
"""
Model complex assembly as a concentration-dependent process.
Returns time in minutes.
"""
# Model assembly time as dependent on protein concentrations
# In a real cell, this would depend on the concentration of each Cas protein
# Here we use a simplified model with stochasticity
# Base assembly time is now adjusted to account for the additional Cas11e protein
# More proteins in the complex can increase assembly time due to more interactions
num_proteins = len(self.cas_proteins)
complexity_factor = 1.0 + (num_proteins - 5) * 0.1 # 10% increase for each protein beyond 5
base_time = 1.0 / self.assembly_rate_constant * complexity_factor
# Add effects of molecular crowding and ATP availability
adjusted_time = base_time / (self.molecular_crowding_factor * self.atp_availability)
# Add stochasticity
noise_factor = np.random.normal(1.0, 0.25) # 25% noise
return adjusted_time * noise_factor
def get_dna_degradation_rate(self, gene):
"""
Calculate DNA degradation rate for a specific gene based on its properties.
Returns rate in bp per minute.
"""
# Base rate adjusted for gene-specific factors
gene_data = self.target_genes[gene]
# Adjust rate based on gene accessibility (chromatin state, etc.)
accessibility_factor = gene_data['accessibility']
# Adjust rate based on gene expression level (affects DNA compaction)
expression_factors = {'low': 0.8, 'medium': 1.0, 'high': 1.2}
expression_factor = expression_factors[gene_data['expression_level']]
# Adjust rate based on ATP availability (Cas3 is ATP-dependent)
atp_factor = self.atp_availability
# Calculate adjusted rate
adjusted_rate = self.base_dna_degradation_rate * accessibility_factor * expression_factor * atp_factor
# Add stochasticity
noise_factor = np.random.normal(1.0, 0.15) # 15% noise
return adjusted_rate * noise_factor
def calculate_binding_probability(self, gene):
"""
Calculate the probability of Cas complex binding to a specific gene.
Returns probability (0-1).
"""
# Calculate binding probability based on binding energy and threshold
gene_data = self.target_genes[gene]
binding_energy = gene_data['binding_energy']
# Binding is more likely if binding energy exceeds threshold
if binding_energy > self.binding_energy_threshold:
# Sigmoid function to model probability
probability = 1.0 / (1.0 + np.exp(-(binding_energy - self.binding_energy_threshold) * 10))
else:
probability = 0.1 # low but non-zero probability for off-target binding
# Adjust for gene accessibility
probability *= gene_data['accessibility']
# Cas11e can enhance target recognition and binding
# This is modeled as a slight increase in binding probability
cas11e_enhancement_factor = 1.15 # 15% enhancement due to Cas11e
probability *= cas11e_enhancement_factor
# Cap probability at 1.0
probability = min(probability, 1.0)
return probability
def select_next_gene(self, available_genes):
"""
Select the next gene to cleave based on binding probabilities.
Returns gene name.
"""
if not available_genes:
return None
# Calculate binding probability for each available gene
probabilities = {gene: self.calculate_binding_probability(gene) for gene in available_genes}
# Normalize probabilities
total_probability = sum(probabilities.values())
if total_probability > 0:
normalized_probabilities = {gene: prob/total_probability for gene, prob in probabilities.items()}
else:
# Equal probability if all are zero
normalized_probabilities = {gene: 1.0/len(available_genes) for gene in available_genes}
# Convert to list for random selection
genes = list(normalized_probabilities.keys())
probs = list(normalized_probabilities.values())
# Select gene based on probabilities
selected_gene = np.random.choice(genes, p=probs)
return selected_gene
def run_simulation(self):
"""Run the CRISPR-Cas3 expression and activation simulation using mechanistic model."""
# Reset target genes status
for gene in self.target_genes:
self.target_genes[gene]['status'] = 'intact'
self.target_genes[gene]['cleavage_started'] = None
self.target_genes[gene]['cleavage_completed'] = None
# Clear previous results
self.results = []
# Model transcription initiation as a stochastic process
transcription_start = self.get_transcription_start_time()
print(f"Transcription starts at: {transcription_start:.2f} minutes post-infection")
# Model translation of each Cas protein
translation_times = {protein: self.get_translation_time(protein) for protein in self.cas_proteins}
total_translation_time = max(translation_times.values())
print(f"Translation completed at: {transcription_start + total_translation_time:.2f} minutes")
print("Individual protein translation times:")
for protein, time in translation_times.items():
print(f" {protein}: {time:.2f} minutes")
# Model protein folding
folding_times = {protein: self.get_protein_folding_time(protein) for protein in self.cas_proteins}
total_folding_time = max(folding_times.values())
print(f"Protein folding completed at: {transcription_start + total_translation_time + total_folding_time:.2f} minutes")
print("Individual protein folding times:")
for protein, time in folding_times.items():
print(f" {protein}: {time:.2f} minutes")
# Model complex assembly
assembly_time = self.get_complex_assembly_time()
cas_activation_time = transcription_start + total_translation_time + total_folding_time + assembly_time
print(f"Cas complex assembly time: {assembly_time:.2f} minutes")
print(f"Cas complex assembled and activated at: {cas_activation_time:.2f} minutes")
# Initialize current time and active cleavage
current_time = 0
active_cleavage = None
# Simulate each time step
for time in self.time_points:
current_time = time
# Record state at this time point for all genes
for gene_name, gene_data in self.target_genes.items():
self.results.append({
'time': time,
'gene': gene_name,
'cleavage_started': gene_data['cleavage_started'],
'cleavage_completed': gene_data['cleavage_completed'],
'status': gene_data['status']
})
# Check if Cas3 is activated
if time < cas_activation_time:
continue
# If no active cleavage, select next gene to cleave
if active_cleavage is None:
# Get list of intact genes
intact_genes = [gene for gene, data in self.target_genes.items()
if data['status'] == 'intact']
if intact_genes:
# Select next gene based on binding probabilities
next_gene = self.select_next_gene(intact_genes)
# Start cleavage of this gene
self.target_genes[next_gene]['status'] = 'cleaving'
self.target_genes[next_gene]['cleavage_started'] = time
# Calculate degradation rate for this gene
degradation_rate = self.get_dna_degradation_rate(next_gene)
# Calculate completion time based on gene length and degradation rate
completion_time = time + (self.target_genes[next_gene]['length'] / degradation_rate)
self.target_genes[next_gene]['cleavage_completed'] = completion_time
# Set as active cleavage
active_cleavage = next_gene
print(f"Started cleaving {next_gene} at {time:.2f} minutes, expected completion at {completion_time:.2f} minutes")
# Check if active cleavage is completed
if active_cleavage and time >= self.target_genes[active_cleavage]['cleavage_completed']:
print(f"Completed cleaving {active_cleavage} at {time:.2f} minutes")
self.target_genes[active_cleavage]['status'] = 'disrupted'
(Content truncated due to size limit. Use line ranges to read in chunks)Initial Parameter Setup
Initial Bacterial Population: 1.0 × 10⁸ CFU/mL
Initial Phage Concentration: 1.0 × 10⁹ PFU/mL
Phage Adsorption Rate: 2.0 × 10^-9 mL/min per bacterium
CRISPR-Cas3 Delivery Efficiency: 85%
CRISPR-Cas3 Cleavage Efficiency: 95%
Simulation Start: 0 minutes (post-infection)
Simulation End: 180 minutes
Time Step: 0.1 minutes
Data Recording Interval: 1 minute
Essential Genes Tracked: ftsZ, murA, rpoB, recA, gyrA
Cas Protein Parameters
Protein Sizes (amino acids):
- Cas3: 900
- Cas5: 350
- Cas6e: 350
- Cas7: 300
- Cas8e: 825
- Cas11e: 275
Protein Folding Efficiency: 0.92 ± 0.05, which is the rate at which the chains of amino acids for each of the Cas genes are able to fold into their correct 3D shape and become functional and active.
Complex Assembly Efficiency: 0.95 ± 0.03, which is the rate at which the Cas genes can assemble to form the entire CRISPR-Cas3 complex.
Target Gene Parameters
Target Gene Lengths:
- ftsZ: 1152 bp
- gyrA: 2628 bp
- murA: 1260 bp
- recA: 1062 bp
- rpoB: 4062 bp
Gene Targeting Order: This is based on how accessible they are in the genome and their binding levels.
- ftsZ
- murA
- rpoB
- recA
- gyrA
Bacterial Cell Death Parameters
Death Rate With One Essential Gene Disrupted: 0.03 ± 0.01 per minute
Death Rate With Multiple Genes Disrupted: 0.23 ± 0.05 per minute
Stochastic Parameters
The simulations were run 5 times with different variations for biological variability and for more realism.
How Effective is the Entire Treatment Process (From Initial Infection to Gene Knockout) at Degrading the Populations of E. coli CFT073?
Here is the code used for simulation 3:
#!/usr/bin/env python3.11
"""
Full Infection Clearance Simulation with Engineered Killian Phages Carrying CRISPR-Cas3
This script simulates the complete process of infection clearance in a host treated with
engineered Killian phages carrying CRISPR-Cas3. The simulation integrates phage infection dynamics,
CRISPR-Cas3 activation, and bacterial death kinetics to model the full 24-hour treatment process.
All parameters are based on literature values and previous simulations for consistency and accuracy.
Updated to include Cas11e and new Cas protein lengths from simulation 3.
"""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import os
from datetime import datetime
import random
from collections import Counter
import matplotlib.dates as mdates
from matplotlib.ticker import PercentFormatter
import matplotlib.gridspec as gridspec
# Set random seed for reproducibility
np.random.seed(42)
class InfectionClearanceSimulation:
def __init__(self,
initial_bacteria=1e8, # CFU/mL
initial_phage=1e9, # PFU/mL
simulation_time=1440, # minutes (24 hours)
time_step=5, # minutes
output_dir="/home/ubuntu",
configurable_params=None):
"""
Initialize the infection clearance simulation with parameters from literature and previous simulations.
Parameters:
-----------
initial_bacteria : float
Initial bacterial concentration (CFU/mL)
initial_phage : float
Initial phage concentration (PFU/mL)
simulation_time : int
Total simulation time in minutes
time_step : int
Time step for simulation in minutes
output_dir : str
Directory to save output files
configurable_params : dict or None
Optional dictionary of configurable parameters to override defaults
"""
self.initial_bacteria = initial_bacteria
self.initial_phage = initial_phage
self.simulation_time = simulation_time
self.time_step = time_step
self.output_dir = output_dir
# Initialize parameters with default values
self.params = self._initialize_parameters()
# Override defaults with configurable parameters if provided
if configurable_params:
for key, value in configurable_params.items():
if key in self.params:
self.params[key] = value
print(f"Parameter {key} set to {value}")
else:
print(f"Warning: Parameter {key} not recognized and will be ignored")
# Initialize data structures
self.time_points = np.arange(0, simulation_time + time_step, time_step)
self.results = None
self.detailed_results = None
self.clearance_metrics = None
def _initialize_parameters(self):
"""
Initialize simulation parameters based on literature and previous simulations.
Returns:
--------
dict
Dictionary of simulation parameters
"""
params = {
# Bacterial parameters
"bacterial_doubling_time": 22.4, # minutes, from Forsyth et al. (2018)
# Phage parameters
"adsorption_rate": 2.0e-9, # mL/min per bacterium
"adsorption_duration": 20, # minutes
"burst_size": 139, # phage particles per infected cell
"latent_period": 30, # minutes
"phage_half_life": 24 * 60, # minutes (24 hours)
# CRISPR-Cas3 parameters
"crispr_delivery_efficiency": 0.85, # probability of successful delivery
"crispr_cleavage_efficiency": 0.95, # probability of successful cleavage
"transcription_start_min": 5, # minutes post-infection
"transcription_start_max": 10, # minutes post-infection
"cas_assembly_time": 5, # minutes
"dna_degradation_rate": 1000, # bp per minute
# Updated gene parameters based on simulation 3 with Cas11e
# Times are updated to reflect earlier activation (~25 minutes sooner)
"target_genes": {
# Updated cleavage times based on simulation 3 with Cas11e
"ftsZ": {"length": 1152, "cleavage_start": 122.4, "cleavage_end": 123.8},
"recA": {"length": 1062, "cleavage_start": 133.2, "cleavage_end": 135.2},
"murA": {"length": 1260, "cleavage_start": 124.0, "cleavage_end": 126.4},
"gyrA": {"length": 2628, "cleavage_start": 135.4, "cleavage_end": 137.9},
"rpoB": {"length": 4062, "cleavage_start": 126.6, "cleavage_end": 133.1}
},
# Bacterial death parameters
"death_rate_one_gene": {
"mean": 0.03, # per minute (90% death within ~30 minutes)
"std": 0.01 # standard deviation to account for cell-to-cell variability
},
"death_rate_multiple_genes": {
"mean": 0.23, # per minute (100% death within ~10 minutes)
"std": 0.05 # standard deviation
},
"death_distribution": "weibull", # Options: "exponential", "weibull"
"weibull_shape": 2.2, # Shape parameter for Weibull distribution
# Stochasticity parameters
"cell_to_cell_variability": 0.15, # Coefficient of variation
# Cas protein parameters (for documentation)
"cas_proteins": {
"Cas3": 900, # 900 amino acids
"Cas5": 350, # 350 amino acids
"Cas6e": 350, # 350 amino acids
"Cas7": 300, # 300 amino acids
"Cas8e": 825, # 825 amino acids
"Cas11e": 275 # 275 amino acids - newly added
},
# Enhanced binding due to Cas11e
"binding_enhancement": 1.15 # 15% enhancement in binding probability
}
return params
def _bacterial_growth_rate(self):
"""
Calculate bacterial growth rate based on doubling time.
Returns:
--------
float
Growth rate in min^-1
"""
return np.log(2) / self.params["bacterial_doubling_time"]
def _phage_decay_rate(self):
"""
Calculate phage decay rate based on half-life.
Returns:
--------
float
Decay rate in min^-1
"""
return np.log(2) / self.params["phage_half_life"]
def _calculate_death_probability(self, disrupted_genes_count, time_since_disruption):
"""
Calculate the probability of cell death based on the number of disrupted genes
and time since disruption.
Parameters:
-----------
disrupted_genes_count : int
Number of disrupted essential genes
time_since_disruption : float
Time in minutes since gene disruption
Returns:
--------
float
Probability of cell death in the current time step
"""
if disrupted_genes_count == 0:
return 0.0
# Select appropriate death rate based on number of disrupted genes
if disrupted_genes_count == 1:
base_rate = np.random.normal(
self.params["death_rate_one_gene"]["mean"],
self.params["death_rate_one_gene"]["std"]
)
else: # 2 or more genes
base_rate = np.random.normal(
self.params["death_rate_multiple_genes"]["mean"],
self.params["death_rate_multiple_genes"]["std"]
)
# Ensure rate is positive
base_rate = max(0.001, base_rate)
# Calculate probability based on distribution type
if self.params["death_distribution"] == "exponential":
# Simple exponential decay model
prob = 1 - np.exp(-base_rate * self.time_step)
elif self.params["death_distribution"] == "weibull":
# Weibull distribution - more realistic for biological systems
# Cumulative distribution function: 1 - exp(-(t/λ)^k)
scale = 1.0 / base_rate # λ parameter (scale)
shape = self.params["weibull_shape"] # k parameter (shape)
# Probability of death between t and t+Δt
prob = (1 - np.exp(-(((time_since_disruption + self.time_step) / scale) ** shape))) - \
(1 - np.exp(-((time_since_disruption / scale) ** shape)))
# Ensure probability is within bounds
prob = min(1.0, max(0.0, prob))
else:
# Default to exponential if distribution not recognized
prob = 1 - np.exp(-base_rate * self.time_step)
return prob
def _determine_gene_disruption_status(self, infection_time, current_time, cell_variability_factor=1.0):
"""
Determine which genes are disrupted at a given time point after infection.
Parameters:
-----------
infection_time : float
Time when the bacterium was infected (minutes)
current_time : float
Current simulation time (minutes)
cell_variability_factor : float
Factor to account for cell-to-cell variability
Returns:
--------
dict
Dictionary with gene names as keys and disruption status as values
"""
# Time since infection
time_since_infection = current_time - infection_time
# If not enough time has passed for CRISPR-Cas3 activation, no genes are disrupted
if time_since_infection < 0:
return {gene: "intact" for gene in self.params["target_genes"]}
# Determine transcription start time (stochastic)
transcription_start = infection_time + random.uniform(
self.params["transcription_start_min"],
self.params["transcription_start_max"]
) * cell_variability_factor
# Calculate Cas complex assembly completion time
cas_assembly_complete = transcription_start + self.params["cas_assembly_time"] * cell_variability_factor
gene_status = {}
for gene, gene_data in self.params["target_genes"].items():
# Apply cell-to-cell variability to cleavage timing
cleavage_start = infection_time + gene_data["cleavage_start"] * cell_variability_factor
cleavage_end = infection_time + gene_data["cleavage_end"] * cell_variability_factor
if current_time < cleavage_start:
# Gene not yet targeted
gene_status[gene] = "intact"
elif cleavage_start <= current_time < cleavage_end:
# Gene being cleaved
gene_status[gene] = "cleaving"
else:
# Gene fully disrupted
gene_status[gene] = "disrupted"
return gene_status
def _count_disrupted_genes(self, gene_status):
"""
Count the number of disrupted genes.
Parameters:
-----------
gene_status : dict
Dictionary with gene names as keys and disruption status as values
Returns:
--------
int
Number of disrupted genes
"""
return sum(1 for status in gene_status.values() if status == "disrupted")
def run_simulation(self):
"""
Run a single simulation of infection clearance with engineered phages.
Returns:
--------
pandas.DataFrame
Simulation results
"""
# Initialize tracking variables
time_series_data = []
# Initialize bacterial population
bacteria = []
for i in range(int(min(100000, self.initial_bacteria))): # Cap for computational efficiency
scaling_factor = self.initial_bacteria / min(100000, self.initial_bacteria)
bacterium = {
"id": i,
"state": "uninfected", # uninfected, infected, infected_with_crispr, non_dividing, dead
"infection_time": None, # Time when bacterium was infected
"disrupted_genes": {}, # Will track which genes are disrupted
"time_since_first_disruption": 0, # Time since first gene disruption
"variability_factor": np.random.normal(1.0, self.params["cell_to_cell_variability"]), # Cell-to-cell variability
"progeny_release_time": None, # Time when phage progeny will be released
"scaling": scaling_factor # How many real bacteria this simulated one represents
}
bacteria.append(bacterium)
# Initialize phage population
phage_count = self.initial_phage
# Simulation loop
for time_point in self.time_points:
# Count bacteria in each state
uninfected_count = 0
infected_count = 0
infected_with_crispr_count = 0
non_dividing_count = 0
dead_count = 0
# Track phage progeny release in this time step
new_phages = 0
# Update each bacterium
for bacterium in bacteria:
# Skip dead bacteria
if bacterium["state"] == "dead":
dead_count += bacterium["scaling"]
continue
# Check for phage progeny release
if bacterium["state"] in ["infected", "infected_with_crispr"] and \
bacterium["progeny_release_time"] is not None and \
time_point >= bacterium["progeny_release_time"]:
# Release phage progeny
new_phages += self.params["burst_size"] * bacterium["scaling"]
# Bacterium dies after lysis
bacterium["state"] = "dead"
dead_count += bacterium["scaling"]
continue
# Process bacteria with CRISPR-Cas3
if bacterium["state"] == "infected_with_crispr":
# Determine gene disruption status
gene_status = self._determine_gene_disruption_status(
bacterium["infection_time"],
time_point,
bacterium["variability_factor"]
)
# Update bacterium's disrupted genes
bacterium["disrupted_genes"] = gene_status
# Count disrupted genes
disrupted_count = self._count_disrupted_genes(gene_status)
# Update bacterium state based on disrupted genes
if disrupted_count == 0:
# No genes disrupted yet, still infected with CRISPR
infected_with_crispr_count += bacterium["scaling"]
else:
# Update time since first disruption
if bacterium["time_since_first_disruption"] == 0 and disrupted_count > 0:
bacterium["time_since_first_disruption"] = self.time_step
else:
bacterium["time_since_first_disruption"] += self.time_step
(Content truncated due to size limit. Use line ranges to read in chunks)For this simulation, the parameters used in simulations 1 and 2 were reused, as well as the results from simulations 1 and 2. That’s because we want to put the entire treatment process together from start to finish and see how effective and efficient it is at degrading the bacterial population.
References
Sinkunas, T., Gasiunas, G., Fremaux, C., Barrangou, R., Horvath, P., & Siksnys, V. (2013). Cas3 is a single-stranded DNA nuclease and ATP-dependent helicase in the CRISPR/Cas immune system. The EMBO Journal, 32(3), 385–394.
Redding, S., Sternberg, S. H., Marshall, M., Gibb, B., Bhat, P., Guegler, C. K., Wiedenheft, B., Doudna, J. A., & Greene, E. C. (2015). Surveillance and Processing of Foreign DNA by the Escherichia coli CRISPR-Cas System. Cell, 163(4), 854–865.
Hochstrasser, M. L., Taylor, D. W., Bhat, P., Guegler, C. K., Sternberg, S. H., Nogales, E., & Doudna, J. A. (2014). CasA mediates Cas3-catalyzed target degradation during CRISPR RNA-guided interference. Proceedings of the National Academy of Sciences, 111(17), 6618–6623.
Sternberg, S. H., Redding, S., Jinek, M., Greene, E. C., & Doudna, J. A. (2014). DNA interrogation by the CRISPR RNA-guided endonuclease Cas9. Nature, 507(7490), 62–67.
Qiu, P., Shandilya, H., D’Alessio, J. M., O’Connor, K., Durocher, J., & Gerard, G. F. (2013). The CRISPR/Cas9 system: a powerful tool for genome engineering. Chinese Journal of Biotechnology, 29(7), 957–968.
Levy, Y., Onuchic, J. N., & Wolynes, P. G. (2012). Quantitative analysis of protein folding and stability in the lens capsule. Journal of Biological Chemistry, 287(46), 38085–38093.
Eaton, W. A., Muñoz, V., Hagen, S. J., Jas, G. S., Lapidus, L. J., Henry, E. R., & Hofrichter, J. (2000). Fast kinetics and mechanisms in protein folding. Annual Review of Biophysics and Biomolecular Structure, 29, 327–359.
Ellis, R. J. (2001). Macromolecular crowding: obvious but underappreciated. Trends in Biochemical Sciences, 26(10), 597–604.
Taniguchi, Y., Choi, P. J., Li, G. W., Chen, H., Babu, M., Hearn, J., Emili, A., & Xie, X. S. (2010). Quantifying E. coli proteome and transcriptome with single-molecule sensitivity in single cells. Science, 329(5991), 533–538.
Golding, I., Paulsson, J., Zawilski, S. M., & Cox, E. C. (2005). Real-time kinetics of gene activity in individual bacteria. Cell, 123(6), 1025–1036.
Ingolia, N. T., Lareau, L. F., & Weissman, J. S. (2011). Ribosome profiling of mouse embryonic stem cells reveals the complexity and dynamics of mammalian proteomes. Cell, 147(4), 789–802.
Kahramanoglou, C., Seshasayee, A. S. N., Prieto, A. I., Ibberson, D., Schmidt, S., Zimmermann, J., Benes, V., Fraser, G. M., & Luscombe, N. M. (2011). Direct and indirect effects of H-NS and Fis on global gene expression control in Escherichia coli. Nucleic Acids Research, 39(6), 2073–2091.
Schink, S. J., Biselli, E., Ammar, C., & Gerland, U. (2019). Death Rate of E. coli during Starvation Is Set by Maintenance Cost and Biomass Recycling. Cell Systems, 9(1), 64–73.
Gruber, J. D., Chen, W., Nilsen, S., Grandjean, J. M. D., Li, X., Tan, Y., Burroughs, A. M., Gao, Y., Borkovich, K. A., & Xu, W. (2022). Degradation of the Escherichia coli Essential Proteins DapB and Dxr Results in Oxidative Stress, which Contributes to Lethality through Incomplete Base Excision Repair. mBio, 13(1), e03756–21.
Forsyth, V. S., Armbruster, C. E., Smith, S. N., Pirani, A., Springman, A. C., Walters, M. S., Nielubowicz, G. R., Himpsl, S. D., Snitkin, E. S., & Mobley, H. L. T. (2018). Factors influencing bacterial growth rates in the human gut microbiome. Nature Microbiology, 3(8), 979–988.
Milo, R., & Phillips, R. (2015). Cell Biology by the Numbers. Garland Science.
UniProt Consortium (2021). UniProt: the universal protein knowledgebase in 2021. Nucleic Acids Research, 49(D1), D480-D489.
E. coli CFT073 genome sequence (GenBank accession number AE014075.1) — Used for direct sequence analysis of target gene lengths.
Gillespie, D. T. (1977). Exact stochastic simulation of coupled chemical reactions. The Journal of Physical Chemistry, 81(25), 2340–2361.
Karr, J. R., Sanghvi, J. C., Macklin, D. N., Gutschow, M. V., Jacobs, J. M., Bolival, B., Assad-Garcia, N., Glass, J. I., & Covert, M. W. (2012). A whole-cell computational model predicts phenotype from genotype. Cell, 150(2), 389–401.
Hucka, M., Finney, A., Sauro, H. M., Bolouri, H., Doyle, J. C., Kitano, H., Arkin, A. P., Bornstein, B. J., Bray, D., Cornish-Bowden, A., et al. (2003). The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models. Bioinformatics, 19(4), 524–531.
Baba, T., Ara, T., Hasegawa, M., Takai, Y., Okumura, Y., Baba, M., Datsenko, K. A., Tomita, M., Wanner, B. L., & Mori, H. (2006). Construction of Escherichia coli K-12 in-frame, single-gene knockout mutants: the Keio collection. Molecular Systems Biology, 2(1), 2006.0008.
Jiang, W., Bikard, D., Cox, D., Zhang, F., & Marraffini, L. A. (2013). RNA-guided editing of bacterial genomes using CRISPR-Cas systems. Nature Biotechnology, 31(3), 233–239.
Datsenko, K. A., & Wanner, B. L. (2000). One-step inactivation of chromosomal genes in Escherichia coli K-12 using PCR products. Proceedings of the National Academy of Sciences, 97(12), 6640–6645.
Brouns, S. J., Jore, M. M., Lundgren, M., Westra, E. R., Slijkhuis, R. J., Snijders, A. P., Dickman, M. J., Makarova, K. S., Koonin, E. V., & van der Oost, J. (2008). Small CRISPR RNAs guide antiviral defense in prokaryotes. Science, 321(5891), 960–964.
메타데이터
- post_id
- 280a0c9de34a
- slug
- 20-min-read-creating-a-crispr-cas3-phage-vector-to-knock-out-the-essential-genes-of-e-280a0c9de34a
- url
- https://medium.com/@IsavellaT/20-min-read-creating-a-crispr-cas3-phage-vector-to-knock-out-the-essential-genes-of-e-280a0c9de34a
- canonical_url
- https://medium.com/@IsavellaT/20-min-read-creating-a-crispr-cas3-phage-vector-to-knock-out-the-essential-genes-of-e-280a0c9de34a
- author_url
- https://medium.com/@IsavellaT
- status
- ok
- fetched_at
- 2026-07-19 00:05:24
