Automatic screening method, system, device and storage medium for gene editing sites

Abstract

The present invention discloses an automatic screening method, system, device and storage medium for gene editing sites, the method comprising: obtaining the gene name and species name, and obtaining the sequence, transcript information and position information of the gene on the chromosome from a gene database; determining the candidate knocking region according to the position information and transcript information; performing GC content and complexity analysis on the candidate knocking region to determine the knocking region; sorting the knocking region according to the number of common transcripts and position information; according to the sorting, intercepting the knocking region sequence from the sequence of the gene through the position information of the knocking region; obtaining all gRNAs and corresponding specificity scores and cutting efficiency scores of the knocking region sequence from the gRNA online design software, and determining and displaying the editing site of the gRNA according to the specificity score and cutting efficiency score. The embodiment of the present invention can automatically and efficiently screen gene editing sites with good knockout effect, strong specificity and high cutting efficiency, and can be widely used in the field of gene editing technology.

Description

Technical Field

[0001] The present invention relates to the field of gene editing technology, and in particular to an automatic screening method, system, device and storage medium for gene editing sites. Background Art

[0002] CRISPR / Cas9 uses a specific guide RNA to recognize the target sequence and guide the Cas9 endonuclease to cut the PAM upstream of the target sequence, causing a double-strand break in the target DNA. After the DNA double-strand break, the cleavage site is repaired by non-homologous end joining (NHEJ) or homologous recombination (HDR) in the cell to achieve gene knockout, knock-in or point mutation at the DNA level. CRISPR / Cas9 has become one of the hottest technologies in the field of life sciences due to its advantages such as simple operation, low cost and editing efficiency. It has been widely used in gene editing-related functional studies of many model organisms, such as mammals (rats, mice, pigs, rabbits, monkeys, etc.), zebrafish, stem cells, tumor cell lines, bacteria and fungi, etc.

[0003] The design of CRISPR / Cas9 knockout scheme seems to be easy, but if you want to achieve efficient cutting and good subsequent protein-level knockout effect, you need to consider many influencing factors. For gene editing beginners, it takes several days and the effect may not be ideal; for researchers with gene editing experience, they may not choose gRNA with high cutting efficiency or protein expression remains unchanged after DNA base deletion in the selected knockout region due to incomplete information sources. Summary of the invention

[0004] In view of this, the purpose of the embodiments of the present invention is to provide an automatic screening method, system, device and storage medium for gene editing sites, which can automatically and efficiently screen gene editing sites with good knockout effect, strong specificity and high cutting efficiency.

[0005] In a first aspect, an embodiment of the present invention provides a method for automatically screening gene editing sites, comprising the following steps:

[0006] Obtaining a gene name and a species name, and obtaining the sequence, transcript information, and chromosome location information of the gene from a gene database according to the gene name and the species name;

[0007] Determine a candidate knocking area according to the position information and the transcript information;

[0008] Performing GC content analysis and complexity analysis on the candidate knocking region to determine the knocking region;

[0009] sorting the knockout regions according to the number of shared transcripts and the position information;

[0010] According to the sorting, the knock region sequence is intercepted from the sequence of the gene through the position information of the knock region;

[0011] All gRNAs of the knock region sequence and the corresponding specificity scores and cutting efficiency scores are obtained from the gRNA online design software, and the editing site of the gRNA is determined and displayed based on the specificity score and the cutting efficiency score.

[0012] Optionally, determining the candidate knocking area according to the position information and according to the transcript information comprises the steps of:

[0013] Determine the range of the tapping area according to the position information;

[0014] Preprocessing the transcripts in the knockout region to obtain reference transcripts;

[0015] The alternative knockout region is determined based on the position and length of the protein coding region of the reference transcript and / or the position and length of the exon.

[0016] Optionally, the alternative knock-in region is determined based on the position and length of the protein coding region and / or the position and length of the exon of the reference transcript, including: selecting a region at a preset position and preset length of the reference transcript as the alternative knock-in region.

[0017] Optionally, the alternative knock-in region is determined based on the position and length of the protein coding region of the reference transcript and / or the position and length of the exon, including: determining the alternative knock-in region based on the position of the gRNA in the exon.

[0018] Optionally, the alternative knock region is determined based on the position and length of the protein coding region of the reference transcript and / or the position and length of the exon, including: after merging other exons within a preset range of the exon, the alternative knock region is determined based on the position and length of the merged exons in the coding region of the reference transcript.

[0019] Optionally, performing GC content analysis and complexity analysis on the candidate knocking region to determine the knocking region includes:

[0020] The average GC content within the preset region length selected in the candidate knocking area is within the preset interval;

[0021] An area within the preset interval without a preset complex sequence is selected as the tapping area.

[0022] Optionally, the sorting of the knocking regions according to the number of shared transcripts and the position information comprises:

[0023] The knocked regions are ranked according to the number of shared transcripts;

[0024] If the number of shared transcripts is the same, the knockout regions are sorted according to the position information.

[0025] Optionally, determining the editing site of the gRNA according to the specificity score and the cleavage efficiency score comprises:

[0026] Selecting a gRNA whose specificity score and the cleavage efficiency score are within a preset range;

[0027] According to the ranking of the specificity score and the cutting efficiency score, the editing sites of several gRNAs are determined within the preset range.

[0028] In a second aspect, an embodiment of the present invention provides an automatic screening device for gene editing sites, comprising:

[0029] at least one processor;

[0030] at least one memory for storing at least one program;

[0031] When the at least one program is executed by the at least one processor, the at least one processor implements the automatic screening method for gene editing sites described in the embodiment of the first aspect.

[0032] In a third aspect, an embodiment of the present invention provides a storage medium storing a program executable by a processor, wherein the program executable by the processor is used to perform the automatic screening method for gene editing sites as described in the embodiment of the first aspect when executed by the processor.

[0033] In a fourth aspect, an embodiment of the present invention provides an automatic screening system for gene editing sites, comprising a computer device, a gene database connected to the computer device, and gRNA online design software; wherein,

[0034] The gene database is used to obtain the sequence and information of the gene, the transcript information and the location information on the chromosome;

[0035] The gRNA online design software is used to obtain all gRNAs of the knock region sequence and the corresponding specificity scores and cutting efficiency scores;

[0036] The computer device comprises:

[0037] at least one processor;

[0038] at least one memory for storing at least one program;

[0039] When the at least one program is executed by the at least one processor, the at least one processor implements an automatic screening method for gene editing sites as described in the first aspect embodiment.

[0040] Implementation of the embodiments of the present invention includes the following beneficial effects: the embodiments of the present invention determine candidate knocking regions based on the position information and transcript information of the gene on the chromosome, perform GC content analysis and complexity analysis on the candidate knocking regions to determine the knocking regions, sort the knocking regions according to the number of common transcripts and the position information, thereby selecting knocking regions with good effects; according to the sorting, the knocking region sequence is intercepted from the sequence of the gene through the position information of the knocking region, all gRNAs of the knocking region sequence and the corresponding specificity scores and cutting efficiency scores are obtained from the gRNA online design software, and the editing sites of the gRNA are determined and displayed according to the specificity scores and the cutting efficiency scores, thereby determining gene editing sites with strong specificity and high cutting efficiency; in addition, the entire gene screening process is automated, efficient and fast. BRIEF DESCRIPTION OF THE DRAWINGS

[0041] Figure 1 This is a schematic diagram of the steps of an automatic screening method for gene editing sites provided by an embodiment of the present invention;

[0042] Figure 2 It is a structural block diagram of an automatic screening device for gene editing sites provided by an embodiment of the present invention;

[0043] Figure 3 It is a structural block diagram of an automatic screening system for gene editing sites provided in an embodiment of the present invention. DETAILED DESCRIPTION

[0044] The present invention is further described in detail below in conjunction with the accompanying drawings and specific embodiments. The step numbers in the following embodiments are only provided for the convenience of explanation and description, and the order between the steps is not limited in any way. The execution order of each step in the embodiment can be adaptively adjusted according to the understanding of those skilled in the art.

[0045] like Figure 1 As shown, an embodiment of the present invention provides an automatic screening method for gene editing sites, which includes the following steps.

[0046] S1. Obtaining a gene name and a species name, and obtaining the sequence, transcript information and chromosome location information of the gene from a gene database according to the gene name and the species name;

[0047] S2, determining a candidate knocking area according to the position information and the transcript information;

[0048] S3, performing GC content analysis and complexity analysis on the candidate knocking region to determine the knocking region;

[0049] S4, sorting the knocking regions according to the number of shared transcripts and the position information;

[0050] S5. According to the sorting, extracting the knock region sequence from the gene sequence through the position information of the knock region;

[0051] S6. Obtain all gRNAs of the knockout region sequence and the corresponding specificity scores and cutting efficiency scores from the gRNA online design software, and determine and display the editing site of the gRNA based on the specificity scores and the cutting efficiency scores.

[0052] It should be noted that the gene database may be the Ensembl database or other databases, and the embodiment of the present invention does not make any specific limitation.

[0053] It should be noted that the website of the gRNA online design software is http: / / crispor.tefor.net / , and other software may also be used, and the embodiment of the present invention does not make specific limitations.

[0054] Implementation of the embodiments of the present invention includes the following beneficial effects: the embodiments of the present invention determine candidate knocking regions based on the position information of the gene on the chromosome and the transcript information, and perform GC content analysis and complexity analysis on the candidate knocking regions to determine the knocking regions, and sort the knocking regions according to the number of common transcripts and the position information, thereby selecting knocking regions with good effects; according to the sorting, the knocking region sequence is intercepted from the sequence of the gene through the position information of the knocking region, all gRNAs of the knocking region sequence and the corresponding specificity scores and cutting efficiency scores are obtained from the gRNA online design software, and the editing sites of the gRNA are determined and displayed according to the specificity scores and the cutting efficiency scores, thereby determining gene editing sites with strong specificity and high cutting efficiency; in addition, the entire gene editing site screening process is automated, efficient and fast.

[0055] Optionally, determining the candidate knocking area according to the position information and according to the transcript information comprises the steps of:

[0056] Determine the range of the tapping area according to the position information;

[0057] Preprocessing the transcripts in the knockout region to obtain reference transcripts;

[0058] The alternative knockout region is determined based on the position and length of the protein coding region of the reference transcript and / or the position and length of the exon.

[0059] Specifically, the position information refers to the position information of the target gene on the chromosome and the position information of the exon on the transcript, the knockout zone range refers to the range suitable for the knockout plan, pretreatment refers to the exclusion of non-protein coding transcripts and incomplete transcripts, the position of the protein coding region refers to the preset position of the protein coding region in the reference transcript, and the length is calculated based on the starting position and the ending position.

[0060] Optionally, the alternative knock-in region is determined based on the position and length of the protein coding region and / or the position and length of the exon of the reference transcript, including: selecting a region at a preset position and preset length of the reference transcript as the alternative knock-in region.

[0061] Optionally, the alternative knock-in region is determined based on the position and length of the protein coding region of the reference transcript and / or the position and length of the exon, including: determining the alternative knock-in region based on the position of the gRNA in the exon.

[0062] It should be noted that the position of the gRNA in the exon includes the gRNA being located on the exon and the gRNA being located on the introns at both ends of the exon. The method for calculating the length of the alternative knock-in region is different when the gRNA is located at different positions of the exon.

[0063] Optionally, the alternative knock region is determined based on the position and length of the protein coding region of the reference transcript and / or the position and length of the exon, including: after merging other exons within a preset range of the exon, the alternative knock region is determined based on the position of the merged exons in the coding region of the reference transcript.

[0064] It should be noted that after merging other exons within the preset range of the exon, the duplicated regions need to be deleted.

[0065] Optionally, performing GC content analysis and complexity analysis on the candidate knocking region to determine the knocking region includes:

[0066] The average GC content within the preset region length selected in the candidate knocking area is within the preset interval;

[0067] An area within the preset interval without a preset complex sequence is selected as the tapping area.

[0068] It should be noted that the preset area length and the preset length range are set according to specific circumstances, and the embodiments of the present invention do not impose specific restrictions; in addition, complex areas include a variety of situations and can be set according to specific circumstances, and the embodiments of the present invention do not impose specific restrictions.

[0069] Optionally, the sorting of the knocking regions according to the number of shared transcripts and the position information comprises:

[0070] The knocked regions are ranked according to the number of shared transcripts;

[0071] If the number of shared transcripts is the same, the knockout regions are sorted according to the position information.

[0072] It should be noted that when sorting, the number of shared transcripts is given priority, followed by the preceding and following position information.

[0073] Optionally, determining the editing site of the gRNA according to the specificity score and the cleavage efficiency score comprises:

[0074] Selecting a gRNA whose specificity score and the cleavage efficiency score are within a preset range;

[0075] According to the ranking of the specificity score and the cutting efficiency score, the editing sites of several gRNAs are determined within the preset range.

[0076] It should be noted that the number of gRNAs to be determined is determined according to the specific circumstances, such as two or one gRNA, and the embodiments of the present invention do not impose specific limitations.

[0077] like Figure 2 As shown, an embodiment of the present invention provides an automatic screening device for gene editing sites, comprising:

[0078] at least one processor;

[0079] at least one memory for storing at least one program;

[0080] When the at least one program is executed by the at least one processor, the at least one processor implements the steps of the automatic screening method for gene editing sites described in the above method embodiment.

[0081] It can be seen that the contents of the above method embodiments are all applicable to the present device embodiments, the functions specifically implemented by the present device embodiments are the same as those of the above method embodiments, and the beneficial effects achieved are also the same as those achieved by the above method embodiments.

[0082] In addition, the present application embodiment also discloses a computer program product or a computer program, which is stored in a computer-readable storage medium. The processor of the computer device can read the computer program from the computer-readable storage medium, and the processor executes the computer program so that the computer device performs the above-mentioned automatic screening method for gene editing sites. Similarly, the contents in the above-mentioned method embodiment are all applicable to the present storage medium embodiment, and the functions specifically implemented by the present storage medium embodiment are the same as those in the above-mentioned method embodiment, and the beneficial effects achieved are also the same as those achieved by the above-mentioned method embodiment.

[0083] like Figure 3 As shown, an embodiment of the present invention provides an automatic screening system for gene editing sites, including a computer device, a gene database connected to the computer device, and gRNA online design software; wherein,

[0084] The gene database is used to obtain the sequence and information of the gene, the transcript information and the location information on the chromosome;

[0085] The gRNA online design software is used to obtain all gRNAs of the knock region sequence and the corresponding specificity scores and cutting efficiency scores;

[0086] The computer device comprises:

[0087] at least one processor;

[0088] at least one memory for storing at least one program;

[0089] When the at least one program is executed by the at least one processor, the at least one processor implements the above-mentioned automatic screening method for gene editing sites.

[0090] Specifically, the computer device may be different types of electronic devices, including but not limited to desktop computers, laptops and other terminals.

[0091] It can be seen that the contents of the above method embodiments are all applicable to the present system embodiments, the functions specifically implemented by the present system embodiments are the same as those of the above method embodiments, and the beneficial effects achieved are also the same as those achieved by the above method embodiments.

[0092] The automatic screening method of the above-mentioned gene editing sites is described below with a specific embodiment. The automatic screening method of the gene editing sites first determines the knockout region, then determines the gRNA according to the knockout region, and finally displays the determined knockout region and gRNA in a drawing. Among them, according to the different knockout schemes, the methods for determining the knockout region and gRNA are different. This embodiment provides three knockout schemes and their corresponding gRNA determination schemes. The three knockout schemes include a frameshift knockout scheme, a small fragment knockout scheme, and a large fragment knockout scheme. The specific implementation method is as follows.

[0093] It should be noted that before automatically editing genes, it is necessary to obtain the target gene name and species name entered by the user, and the interface provided by the Ensembl database can be used to obtain relevant information of the gene from Ensembl.

[0094] Implementation Method 1: Frameshift Knockout Solution

[0095] (1) Determine the knockout region

[0096] S111, determine whether the target gene is suitable for designing a knockout scheme. The specific method is as follows: obtain the position of the target gene on the chromosome from Ensembl, and search on the chromosome to see if there are other genes in the segment where the target gene is located according to this position; if there are other genes, continue to obtain the types of these other genes overlapping with the target gene, if the other genes are protein-coding genes, then the target gene is not suitable for designing a knockout scheme.

[0097] S112, preprocessing the transcript of the target gene. The specific method is: obtaining the transcript information of the target gene from Ensembl, and excluding the transcript of non-coding protein and incomplete transcript.

[0098] S113, finding the protein coding region of each transcript, i.e., CDS. The specific method is: finding the CDS of the transcript through the start position and end position of each transcript, the start position and end position of the coding region, and the position of the exon of the transcript, wherein the overlapping part of the exon and the coding region is the CDS.

[0099] S114, calculating the common situation of CDS regions in all transcripts, that is, how many transcripts have this CDS region. The specific method is: put the CDS regions of all transcripts into a table, each CDS region has a start position and an end position, traverse the table, and judge according to the positional relationship between CDSs of different transcripts. If the CDS of a transcript overlaps with the CDS of another transcript, then this CDS region is shared by the two transcripts.

[0100] S115, select candidate knockout regions. Select the longest transcript as the reference transcript, and use the first 50% of the CDS of the reference transcript as the candidate knockout region. Then calculate the length of the region by the start and end positions of this region, and select a region with a length of 100-300 bp as the candidate knockout region.

[0101] S116. GC content analysis, where GC content refers to the sum of the proportions of G and C bases in a certain sequence. The specific requirement is that the average GC content within 800 bp upstream and downstream of the knockout region must be between 40% and 60%.

[0102] Complexity analysis: No complex regions can appear upstream and downstream of the knockout region. The following regions need to be excluded: 1. Forward or reverse repeats greater than 50bp. Sequences of 50bp are cut out from the DNA sequence of the entire knockout region and compared with the complete DNA sequence. If it appears twice in the DNA sequence, it is a repeat and the region is excluded; 2. 20 consecutive identical bases. Compare 20 A bases, 20 T bases, 20 G bases and 20 C bases with the DNA sequence respectively. If such a sequence appears, the region is excluded; 3. Hairpin structures with two loops less than 25bp. If a certain sequence is the reverse complementary sequence of another sequence, , the DNA single-stranded molecule folds back on itself to form a hairpin structure, and 10bp sequences are taken from the sequence in turn (record the subscript of the last base) for reverse complementary conversion (A becomes T, T becomes A, G becomes C, C becomes G, and then the sequence is reversed), and then it is used as a pattern to search for the same sequence in the entire sequence, and the index subscript of the sequence identical to the pattern is obtained. If the subscript differs from the subscript of the pattern within 25, it is a hairpin structure with a loop less than 25bp. If such a structure appears twice, the iteration is terminated and the region is excluded.

[0103] S117. Sort the remaining candidate knockout regions according to the number of shared transcripts and the order of their positions on the transcripts, with the ones with more shared transcripts given priority, followed by those with earlier positions, and design gRNAs in that order.

[0104] (2) Determine gRNA

[0105] S211. Obtain the sequence of the target gene from the Ensembl database, extract the sequence of the knockout region from the gene sequence through the start and end points of the knockout region, and then obtain all gRNAs in this region and their specificity scores and off-target conditions from the gRNA online design software (http: / / crispor.tefor.net / ).

[0106] S212. Perform a preliminary screening of gRNA. Select gRNAs with a specificity score (MIT Specificity Score) of more than 80 and a Doench'16 score of more than 50; then align each gRNA with the sequence of the target gene to obtain the position of the gRNA on the gene; perform a GC content analysis again, the GC content of the gRNA adjacent region cannot be higher than 80% or lower than 20%, such as using 50bp as the window size, that is, calculating the GC content every 50bp (1-50, 2-51...), calculate the local GC content of the knockout region and the region 800bp upstream and downstream, find all regions with GC content greater than 80% and less than 20%, and then merge the overlapping regions, such as merging 1-50 and 21-70 into 1-70, if the gRNA falls in these regions, it will be excluded.

[0107] S213, predict the cutting efficiency score of gRNA through online software, the cutting score is required to be above 0.5, and then judge whether the latter gRNA overlaps with the former gRNA in order of the specificity score, that is, judge whether it overlaps by the position of the gRNA obtained in the previous step on the gene, and delete the latter gRNA if there is overlap; then sort the gRNA according to the specificity score and cutting efficiency score, and take the best two gRNAs to complete the design. If two gRNAs cannot be found in the knockout region after screening, change the next region to continue the design.

[0108] Implementation Method 2: Small Fragment Knockout Solution

[0109] (1) Determine the knockout region

[0110] S121, determine whether the target gene is suitable for designing a knockout scheme. The specific method is as follows: obtain the position of the target gene on the chromosome from Ensembl, and search on the chromosome to see if there are other genes in the segment where the target gene is located according to this position; if there are other genes, continue to obtain the types of these other genes overlapping with the target gene, if the other genes are protein-coding genes, then the target gene is not suitable for designing a knockout scheme.

[0111] S122, pre-processing the transcript of the target gene. The specific method is: obtaining the transcript information of the target gene from Ensembl, and excluding the transcript of non-coding protein and incomplete transcript.

[0112] S123, finding the protein coding region of each transcript, i.e., CDS. The specific method is: finding the CDS of the transcript through the start position and end position of each transcript, the start position and end position of the coding region, and the position of the exon of the transcript, wherein the overlapping part of the exon and the coding region is the CDS.

[0113] S124, calculate the common situation of CDS region in all transcripts, that is, how many transcripts have this CDS region. The specific method is: put the CDS regions of all transcripts into a table, each CDS region has a starting position and an ending position, traverse the table, and judge according to the positional relationship between CDSs of different transcripts. If the CDS of a transcript overlaps with the CDS of another transcript, then this CDS region is shared by the two transcripts.

[0114] S125. Select all CDSs located in the first 50% of the reference transcript as candidates for knockout regions. The selection of knockout regions for small-fragment knockout schemes is relatively complex and can be divided into two major categories of regions.

[0115] S125A, set gRNA in exon. Find the first 50% of the transcript by the position of the exon on the transcript, subtract the start position from the end position of the transcript to get the length of the transcript, and add half of the length from the start position to get the 50% position of the transcript. If the end position of the exon is less than the 50% position of the transcript, it is in compliance, and the exon length is 300-500bp, the end position of the exon minus the start position is the length of the exon.

[0116] S125B, set gRNA on the introns at both ends of the exon. Determine whether there are other exons within 200bp before and after each exon. If there are, merge and knock out together, and calculate the length of the merged region and the total length of the exons in this region. Since each exon performs this operation, repeated merging will occur, so the repeated region can be excluded based on the starting and ending positions of the region.

[0117] S126. GC content analysis, where GC content refers to the sum of the proportions of G and C bases in a certain sequence. The specific requirement is that the average GC content within 800bp upstream and downstream of the knockout region must be between 40% and 60%.

[0118] Complexity analysis: No complex regions can appear upstream and downstream of the knockout region. The following regions need to be excluded: 1. Forward or reverse repeats greater than 50bp. Sequences of 50bp are cut out from the DNA sequence of the entire knockout region and compared with the complete DNA sequence. If it appears twice in the DNA sequence, it is a repeat and the region is excluded; 2. 20 consecutive identical bases. Compare 20 A bases, 20 T bases, 20 G bases and 20 C bases with the DNA sequence respectively. If such a sequence appears, the region is excluded; 3. Hairpin structures with two loops less than 25bp. If a certain sequence is the reverse complementary sequence of another sequence, , the DNA single-stranded molecule folds back on itself to form a hairpin structure, and 10bp sequences are taken from the sequence in turn (record the subscript of the last base) for reverse complementary conversion (A becomes T, T becomes A, G becomes C, C becomes G, and then the sequence is reversed), and then it is used as a pattern to search for the same sequence in the entire sequence, and the index subscript of the sequence identical to the pattern is obtained. If the subscript differs from the subscript of the pattern within 25, it is a hairpin structure with a loop less than 25bp. If such a structure appears twice, the iteration is terminated and the region is excluded.

[0119] S127. Sort the remaining candidate knockout regions according to the number of shared transcripts and the order of their positions on the transcripts. The regions with more shared transcripts are given priority, followed by those with earlier positions, and gRNAs are designed in that order.

[0120] (2) Determine gRNA

[0121] According to the classification of knockout regions, the methods for designing gRNA for different types of knockout regions are different. The design steps for setting gRNA in exons are as follows.

[0122] S221A, use the sequence of the exon to obtain all gRNAs in this region and their specificity scores and off-target conditions from the gRNA online design software.

[0123] S222A, perform preliminary screening of gRNAs. Select gRNAs with a specificity score (MIT Specificity Score) of 80 or above and a Doench'16 score of 50 or above, then align each gRNA with the sequence of the target gene to obtain the position of the gRNA on the gene. Perform another GC content analysis, and the GC content of the gRNA adjacent region cannot be higher than 80% or lower than 20%.

[0124] S223A, predict the gRNA cutting efficiency score, exclude the gRNA with a cutting score lower than 0.5, and then sort by specificity score. Then find a pair of gRNAs that meet the following conditions in order: the knockout size is ≥ 200bp and the specificity score difference is no more than 10 points, and one of them is located in the first 1 / 2 of the coding region. If there is no gRNA that meets the requirements, replace the next knockout region.

[0125] The design steps for placing gRNA on the introns at both ends of the exon are as follows.

[0126] S221B, use the upstream and downstream sequences of the exon to obtain their gRNAs, their specificity scores, and off-target conditions from the gRNA online design software.

[0127] S222B, perform preliminary screening of gRNA: select gRNAs with a specificity score (MIT Specificity Score) above 80 and a Doench'16 score above 50, then align each gRNA with the sequence of the target gene to obtain the position of the gRNA on the gene; perform a GC content analysis again, the GC content of the gRNA adjacent region cannot be higher than 70% or lower than 30%, such as using 30bp as the window size, that is, calculating the GC content every 30bp (1-30, 2-31...), calculate the local GC content of the knockout region and the region 800bp upstream and downstream, find all regions with GC content greater than 70% and less than 30%, and then merge the overlapping regions, such as merging 1-30 and 21-50 into 1-50, if the gRNA falls in these regions, exclude it.

[0128] S223B, predict the cutting efficiency score of gRNA through online software, exclude gRNA with cutting score lower than 0.5 points, and then sort gRNA according to specificity score, select the gRNA with the highest upstream specificity score and the gRNA with the highest downstream score, and the specificity score difference is within 10 points. First, select a gRNA with the highest specificity score from the upstream, and then select a gRNA with a difference of less than 10 points from the downstream gRNA according to the specificity score sorting; if there is no matching gRNA in the downstream, replace an upstream gRNA according to the specificity score order, and then find a downstream gRNA that meets the requirements.

[0129] Implementation method 3: large fragment knockout solution

[0130] (1) Determine the knockout region

[0131] S131, determine whether the target gene is suitable for designing a knockout scheme. The specific method is as follows: obtain the position of the target gene on the chromosome from Ensembl, and search on the chromosome to see if there are other genes in the segment where the target gene is located according to this position; if there are other genes, continue to obtain the types of these other genes overlapping with the target gene, if the other genes are protein-coding genes, then the target gene is not suitable for designing a knockout scheme.

[0132] S132, preprocessing the transcript of the target gene. The specific method is: obtaining the transcript information of the target gene from Ensembl, and excluding the transcript of non-coding protein and incomplete transcript.

[0133] S133, finding the protein coding region of each transcript, i.e., CDS. The specific method is: finding the CDS of the transcript through the start position and end position of each transcript, the start position and end position of the coding region, and the position of the exon of the transcript, wherein the overlapping part of the exon and the coding region is the CDS.

[0134] S134, calculate the common situation of CDS region in all transcripts, that is, how many transcripts have this CDS region. The specific method is: put the CDS regions of all transcripts into a table, each CDS region has a starting position and an ending position, traverse the table, and judge according to the positional relationship between CDSs of different transcripts. If the CDS of a transcript overlaps with the CDS of another transcript, then this CDS region is shared by the two transcripts.

[0135] S135. Select all CDSs of the reference transcripts, determine whether there are other exons within 200 bp before and after each exon, merge them if there are, and calculate the length of the entire region and the total length of the exons in this region, and then delete the duplicated regions based on the starting and ending positions of the region.

[0136] S136, calculate the 1 / 2 position of the coding region of the reference transcript, the calculation formula: (start+(end-start) / 2), where start represents the starting position and end represents the ending position. If the front end of the exon is located in the front 1 / 2, the front end of these exons can be used as the starting point of the knockout region, and the back end of all exons is used as the end point, and then paired in pairs, requiring the following conditions to be met: the size of the entire region must be between 5kb-15kb. All knockout regions obtained are sorted in the order of position.

[0137] S137. GC content analysis, where GC content refers to the sum of the proportions of G and C bases in a certain sequence. The specific requirement is that the average GC content within 800bp upstream and downstream of the knockout region must be between 40% and 60%.

[0138] Complexity analysis: No complex regions can appear upstream and downstream of the knockout region. The following regions need to be excluded: 1. Forward or reverse repeats greater than 50bp. Sequences of 50bp are cut out from the DNA sequence of the entire knockout region and compared with the complete DNA sequence. If it appears twice in the DNA sequence, it is a repeat and the region is excluded; 2. 20 consecutive identical bases. Compare 20 A bases, 20 T bases, 20 G bases and 20 C bases with the DNA sequence respectively. If such a sequence appears, the region is excluded; 3. Hairpin structures with two loops less than 25bp. If a certain sequence is the reverse complementary sequence of another sequence, , the DNA single-stranded molecule folds back on itself to form a hairpin structure, and 10bp sequences are taken from the sequence in turn (record the subscript of the last base) for reverse complementary conversion (A becomes T, T becomes A, G becomes C, C becomes G, and then the sequence is reversed), and then it is used as a pattern to search for the same sequence in the entire sequence, and the index subscript of the sequence identical to the pattern is obtained. If the subscript differs from the subscript of the pattern within 25, it is a hairpin structure with a loop less than 25bp. If such a structure appears twice, the iteration is terminated and the region is excluded.

[0139] (2) Determine gRNA

[0140] S231. Use the upstream and downstream sequences of the exon to obtain their gRNAs and their specificity scores from the gRNA online design software.

[0141] S232. Perform preliminary screening of gRNAs: select gRNAs with a specificity score (MIT Specificity Score) of more than 80 and a Doench'16 score of more than 50; then align each gRNA with the sequence of the target gene to obtain the position of the gRNA on the gene; perform another GC content analysis, the GC content of the gRNA adjacent region cannot be higher than 70% or lower than 30%, such as using 30bp as the window size, that is, calculating the GC content every 30bp (1-30, 2-31...), calculate the local GC content of the knockout region and the region 800bp upstream and downstream, find all regions with GC content greater than 70% and less than 30%, and then merge the overlapping regions (such as 1-30 and 21-50 are merged into 1-50), if the gRNA falls in these regions, it is excluded.

[0142] S233. Predict the cleavage efficiency scores of gRNAs through online software, exclude gRNAs with cleavage scores below 0.5, and then sort the gRNAs according to the specificity scores, select the gRNA with the highest upstream specificity score and the gRNA with the highest downstream score, and the specificity scores should differ by less than 10 points. First, select a gRNA with the highest specificity score from the upstream, and then select a gRNA with a difference of less than 10 points from the downstream gRNAs according to the specificity score sorting; if there is no matching gRNA in the downstream, replace an upstream gRNA according to the specificity score order, and then find a downstream gRNA that meets the requirements.

[0143] (3) Draw and display

[0144] After determining the knocking area and gRNA according to the knocking schemes of the above three implementation methods, draw a picture for display using the following method. You can choose any one of the display pictures for display, or choose multiple display pictures for display.

[0145] S301, a schematic diagram showing all transcripts and boxing the knockout regions.

[0146] The display diagram mainly uses the data information obtained from Ensembl and processes it to draw a display diagram of all transcripts, and then frame them according to the coordinates of the designed gRNA.

[0147] First, determine the coordinates of the display graph. Specifically, the following method can be used: determine the X-axis range based on the length of the gene. The range of each transcript on the X-axis is the length of each transcript. Each transcript is assigned a Y value, where the length of the gene represents the number of bases. In addition, the color can be determined based on the information provided on Ensembl: if the source of the transcript is "ensembl_havana", the color is determined to be orange; if the transcript is a non-coding protein type, it is determined to be blue; the rest of the transcripts are red.

[0148] Then, draw their exons on each transcript. Considering the location of exons and the location of coding regions, exons can be represented by different marks, such as hollow boxes for exons not in coding regions and solid boxes for exons in coding regions; in addition, the name and direction of the transcript can be marked below each transcript, and the color is the same as that of the transcript.

[0149] Finally, based on the determined position of the gRNA, outline the region to be knocked out.

[0150] S302, schematic diagram showing reference transcripts.

[0151] The schematic diagram of the reference transcript is screened among the schematic diagrams of all transcripts, and only the schematic diagram of the reference transcript is drawn, and the knockout region is framed in the diagram.

[0152] S303, display an enlarged image of the knockout area.

[0153] Based on the location of the exons on the reference transcript and the determined gRNA position, only the exons within the gRNA knockout region are drawn.

[0154] S304, displaying the dot matrix analysis graph and GC content graph.

[0155] The upstream and downstream sequences of the selected knockout region are compared by themselves, and the window size is set to a preset number, such as 20, that is, every 20bp is used as a unit for comparison analysis, that is, a 20bp sequence of a sequence and its reverse complementary sequence are compared with another sequence, and the comparison results are marked on the coordinate axis. In addition, the forward sequence and the reverse complementary sequence can be distinguished by color, such as red for forward and green for reverse complementary.

[0156] The GC graph sets the window size to a preset number, such as 30, that is, the GC content is calculated every 30 bp, and then a graph is drawn with the base position as the X-axis and the GC content as the Y-axis.

[0157] The above is a specific description of the preferred implementation of the present invention, but the invention is not limited to the embodiments. Those skilled in the art may make various equivalent modifications or substitutions without violating the spirit of the present invention. These equivalent modifications or substitutions are all included in the scope defined by the claims of this application.

Claims

1. A method for automatically screening gene editing sites, characterized in that: include: Obtaining a gene name and a species name, and obtaining the sequence, transcript information, and chromosome location information of the gene from a gene database according to the gene name and the species name; Determine a candidate knocking area according to the position information and the transcript information; Performing GC content analysis and complexity analysis on the candidate knocking region to determine the knocking region; sorting the knockout regions according to the number of shared transcripts and the position information; According to the sorting, the knock region sequence is intercepted from the sequence of the gene through the position information of the knock region; Obtain all gRNAs of the knock region sequence and the corresponding specificity scores and cutting efficiency scores from the gRNA online design software, and determine and display the editing site of the gRNA according to the specificity score and the cutting efficiency score; Determining the candidate knocking area according to the position information and the transcript information includes the steps of: Determine the range of the tapping area according to the position information; Preprocessing the transcripts in the knockout region to obtain reference transcripts; Determining the candidate knockout region according to the position and length of the protein coding region and / or the position and length of the exon of the reference transcript; The step of performing GC content analysis and complexity analysis on the candidate knocking region to determine the knocking region comprises: The average GC content within the preset region length selected in the candidate knocking area is within the preset interval; Selecting an area within the preset interval without a preset complex sequence as a tapping area; The region of the preset complex sequence includes: a region where the forward repeat sequence or the reverse repeat sequence is greater than 50 bp, or a region where a sequence of 20 identical bases appears continuously, or a region where a reverse complementary sequence of a hairpin structure with a loop portion less than 25 bp appears twice; The step of determining the range of the knocking area according to the position information comprises the steps of: Based on the location information of the gene on the chromosome, the knocking region ranges corresponding to several target genes are retrieved and determined on the chromosome, and there is no protein-coding gene overlapping with the gene in the region where the target gene is located; The step of preprocessing the transcripts in the knocking area to obtain reference transcripts comprises the following steps: Based on the transcript information of several target genes within the knockout region, non-protein-coding and incomplete transcripts are screened to obtain several first target transcripts; Determining the presence of a plurality of second target transcripts in the protein coding region based on the target position information, wherein the target position information includes the start and end positions of the plurality of first target transcripts, the start and end positions of the coding region, and the positions of exons in the plurality of first target transcripts; Determining the second target transcript with the longest length among the plurality of second target transcripts as a reference transcript; Determining the candidate knockout region according to the position and length of the protein coding region of the reference transcript and / or the position and length of the exon comprises the steps of: If the knockout type is a frameshift knockout, a protein coding region in a preset region is obtained as the first knockout region, wherein the preset region is a region from the start position to half the length of the start and end positions in the reference transcript; The target protein coding region in the first knock region whose length is within the first length range is used as an alternative knock region; If the knockout type is a small fragment knockout, the exon position that satisfies the termination position being within the preset region and the length being within the second length range is used as the second knockout region; Merge adjacent exons within a preset length and delete duplicate regions, and use the exon position that satisfies the termination position within the preset region and the length within the second length range after merging as the third knock region; The second striking area and the third striking area are used as candidate striking areas; If the knockout type is a large-fragment knockout, the adjacent exons within a preset length are merged and the repeated regions are deleted, and the target exon whose starting position after merging is located within the preset region is determined; The starting position of the target exon is taken as the starting point, and the ending position of all exons in the reference transcript is taken as the end point for pairing in pairs to form alternative knockout regions within a preset size range and arranged in order.

2. The automatic screening method for gene editing sites according to claim 1, characterized in that: The step of sorting the knocking regions according to the number of shared transcripts and the position information comprises: The knocked regions are ranked according to the number of shared transcripts; If the number of shared transcripts is the same, the knockout regions are sorted according to the position information.

3. The automatic screening method for gene editing sites according to claim 1, characterized in that: Determining the editing site of the gRNA according to the specificity score and the cutting efficiency score comprises: Selecting a gRNA whose specificity score and the cleavage efficiency score are within a preset range; According to the ranking of the specificity score and the cutting efficiency score, the editing sites of several gRNAs are determined within the preset range.

4. An automatic screening device for gene editing sites, characterized in that: include: at least one processor; at least one memory for storing at least one program; When the at least one program is executed by the at least one processor, the at least one processor implements the automatic screening method for gene editing sites as described in any one of claims 1 to 3.

5. A storage medium storing a program executable by a processor, characterized in that: The processor-executable program is used to perform the automatic screening method for gene editing sites as described in any one of claims 1-2 when executed by the processor.

6. An automatic screening system for gene editing sites, characterized in that: It includes a computer device and a gene database and gRNA online design software connected to the computer device; wherein, The gene database is used to obtain the sequence, transcript information and chromosome location information of the gene; The gRNA online design software is used to obtain all gRNAs of the knock region sequence and the corresponding specificity scores and cutting efficiency scores; The computer device comprises: at least one processor; at least one memory for storing at least one program; When the at least one program is executed by the at least one processor, the at least one processor implements an automatic screening method for gene editing sites as described in any one of claims 1 to 3.

Images