WO2013097061A1 - Bs and rrbs sequencing-based bioinformatics analysis method and device - Google Patents

Abstract

Provided in the present invention are a method and a device for analyzing methylation in a genome by bisulfite sequencing or reduced representation bisulfite sequencing. Also provided in the present invention is a bioinformatics analysis method based on bisulfite sequencing or reduced representation bisulfite sequencing. The method comprises: detecting DNA methylation, then performing the bioinformatics analysis, wherein the analysis is an analysis for one or more items selected from the following： data output, sequencing fragments alignment, the coverage situation of methylation sites, the methylation level of the methylation sites, the distribution trend of the whole genome methylation data, the whole genome DNA methylation map, the differentially methylated regions, the statistic for the length of insert fragments, and the coverage of CpG sites.

Description

 Biological information analysis method and device based on BS and RRBS sequencing

Technical field

 The invention belongs to the field of bioinformatics, and performs bioinformatics analysis on the sequencing results of two important DNA thiolation detection methods in epigenetics. Background technique

 Bisulfite-sequencing (BS) is a powerful tool for DNA thiolation research. The unthiolated C thiol group in the genome can be converted to U by bisulfite treatment, and further converted to T by PCR to distinguish it from the thiolated C base. Combining bisulfite treatment with high-throughput sequencing technology, it is possible to map DNA methylation profiles with single-subtraction resolution. The analysis of high-accuracy thiolation modification patterns of specific species is bound to be a milestone in the study of epigenetics, and lays a foundation for cell differentiation, tissue development, and animal and plant breeding research.

 The reduced representative ')·reduced representation bisulfite sequencing (RRBS) is a new method based on BS. The cost of BS sequencing is high. Based on this, RRBS selects parts of the whole genome by enzymatic cleavage technology and then performs BS sequencing. The treatment of DNA by heavy tartaric acid converts ungerminated C in the genome to U and is further converted to T by PCR to distinguish it from the thiolated C base. This method is capable of plotting a single-base resolution DNA thiolated map with a cost advantage.

 In this field, two methods of DNA thiolation detection, such as BS method and RRBS method, are needed for genomic thiolation analysis, and methods for bioinformatics analysis using the sequencing results are needed. Summary of the invention

 In one aspect, the invention relates to a method for analyzing genomic thiolation by heavy bisulfite hydrogenate sequencing or reduced representative bisulfite sequencing, the method comprising the steps of:

 1) subjecting the sample to bisulfite sequencing or reduced representative bisulfite sequencing, and preferably filtering the sequencing sequence to remove unqualified sequences;

2) Determine the thiolation site represented by the above sequencing sequence, the steps are as follows: i) a first sequencing sequence obtained by converting all c (cytosine) in the respective sequencing sequences into T (thymidine) bases; converting all G (guanine) bases in the complementary sequence of the respective sequencing sequences into a second sequencing sequence obtained from the A (adenine) base;

 Ii) converting all C (cytosine;) 1⁄2 in the reference genome sequence of the corresponding species into the first reference sequence generated by T (thymidine) subtraction; all G (guanine) bases in the reference genome sequence of the corresponding species Converting to a second reference sequence generated by A (adenine) subtraction;

 Iii) aligning each of the first sequencing sequence and the second sequencing sequence with the first reference sequence and the second reference sequence, respectively, and obtaining each of the first sequencing sequence and the second sequencing sequence in the first reference sequence and the second Selecting, in the position on the sequence, each of the first sequencing sequence and the second sequencing sequence having unique alignments on the first reference sequence and the second reference sequence;

 IV) determining the thiolation information of the first sequencing sequence and the second sequencing sequence by using the sequence information of the first sequencing sequence and the second sequencing sequence aligned by the second sequencing sequence selected above, as follows: if in the first sequencing sequence and The corresponding positions of the aligned reference genomic sequences are C (cytosine) bases, indicating that they are thiolated at this position of the first sequencing sequence; if at the corresponding position of the second sequencing sequence and the aligned reference genome sequence All of which are G (guanine) bases, indicating that they are thiolated at this position of the second sequencing sequence;

3) For each C site on the reference genomic sequence, the number of sequencing sequences that are thiolated at the corresponding position by the first or second sequencing sequence covering the C site and the number of sequencing sequences that are not thiolated The number of minimum expected sequencing sequences for thiolation of the C site is determined according to the following inequality, ie, the minimum k value satisfying the following formula: 匪C *: Q)p ^k (l - p ~ ^k < 0.01 * mC; wherein

 mC is the number of sequencing sequences in which the first sequencing sequence or the second sequencing sequence covering the C site is thiolated at the corresponding position,

 nmC is the number of sequencing sequences that cover the C-site of the first sequencing sequence or the second sequencing sequence that is not thiolized at the corresponding position,

 n is the total number of sequencing sequences for this C site,

p is the error rate, for example, the above steps 1) - 2) are carried out using unmodified control DNA such as lamda DNA, and the p value is determined according to the conversion of C to T 1⁄2, p is usually 0.001-0.01, if covering The first sequencing sequence or the second sequencing sequence of the C site is thiolated at the corresponding position The number of sequencing sequences is greater than the number of minimal desired sequencing sequences at the site, which is the site of thiolation, preferably forming a result file comprising the following information:

 Chr-a b +/- c d e f g h

 First column Chr-a: chromosome number,

 The second column b: the position of the C site,

 The third column +/-: positive and negative chain information,

 The fourth column c: the type of C site,

 The fifth column d: specific bases,

 Column 6 e: average alignment position of the sequencing sequence,

 Column 7 f: the number of sequencing sequences that are methylated at the corresponding position by the first sequencing sequence or the second sequencing sequence covering the C site,

 The eighth column g: the number of sequencing sequences covering the first or the second sequencing sequence of the C site that are not thiolated at the corresponding position,

 Ninth column h: Number of minimum expected sequencing sequences for thiolation of the C site.

 In another aspect, the invention relates to a method of bioinformatics analysis by detecting a method of DNA thiolation selected from:

 1. Coverage of the thiolation site, the analysis can be performed in the following manner:

 The C site can be divided into three types, namely CG, CHG, and CHH, where H represents a non-G base. From the result file, it can be obtained whether the site is a methylation site, and the site supporting the thiolated reads number greater than the expected value is the thiolation site. The ratio of mC sites to total sites under different coverage depths (the number of times this site is covered by reads) is counted. At the same time, the coverage of the loci under different coverage depths on the genotype is also counted, and the coverage is the proportion of the loci on the whole genome.

 2. The level of thiolation of the thiolation site, which can be performed by: To understand the overall characteristics of the genomic thiolation map and the average thiolation level of the genome, count each thiolated C base The average thiolation level of the base is calculated as:

 The number of reads that support merging / (the number of reads that support merging + the number of reads that do not support methylation) χ100.

At the same time, the statistics reflect the effective coverage of the reads of the thiolated information. The calculation of the effective coverage is based on three aspects: chromosome, gene interval and gene function elements. After the reads with ambiguous ambiguity are excluded, based on the support-supported reads and the reads data that does not support priming, The effective coverage calculation is performed, and the statistical result is obtained.

3. Trends in genome-wide thiolation data distribution, which can be performed in the following manner: In order to obtain the trend of genome-wide thiolation distribution, the number and distribution ratio of CG, CHG and CHH in the thiolated C, ie each type The percentage of thiolated C to the total thiolization C, for example: proportion of mCHG = number of mCHG / total number of mC. At the same time, different distribution types of C 3⁄4 ( CG, CHG and CHH ) have different levels of thiolation between different species and even different cell types of the same species. Therefore, the distribution of the C base thiolation levels of each type (CG, CHG, and CHH) thiolation was counted. Similarly, the thiolation levels of C-magnetic groups (CG, CHG, and CHH) for each distribution type in each region are calculated by chromosome, gene interval, and gene functional elements. The calculation formula is: Average thiolation level = support for methylation The _reads / ⁽ supports simplification of reads+ does not support methylated reads) χ100. Finally, the distribution of bases near the thiolated C at non-CG sites was analyzed, and the probability of occurrence of different thiolation patterns (CHG and CHH) was counted.

 4. Whole genome DNA thiolation map, which can be performed by: Delineating the distribution of thiolated C bases from the chromosomal level. The density of the thiolated C bases on the chromosome was counted by a 10 kb window. At the same time, the density distributions of different types of thiolated C tt (CG, CHG, and CHH) in each lOOkb window were counted. Different genomic regions have different biological functions, and the thiolation distribution is preferably expressed in the form of a heat map to help to further understand the thiolated graphical features of these regions. The transcription of genes is regulated by different regulatory elements, dividing all coding gene sequences into seven different regions of transcriptional elements, here! ^ Statistics on the level of thiolation of different regions of the transcriptional element.

5. Differential thiolated region (DMR) analysis, which can be performed by: Determining the thiolation difference region by the thiolation level of the two samples, as follows: First, identify a CG site as a seed At the locus, two samples were obtained to support the thiolation and the unsupported methylation values at the locus, and then the chi-squared level of the two samples was tested by the chi-square test to see whether the difference was significant or not. The horizontal threshold p is generally 0.05, preferably 0.01, more preferably 0.001, most preferably 0.0001; if the difference is significant, a CG site is extended backward, and two samples are obtained to support thiolation and unsupporting at these two sites, respectively. The 数值 化 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个This site is then used as a seed site to cycle the above process to obtain a differential thiolated region. After finding the differential methylation regions, the corresponding genes, namely the thiolated differential genes, are found through these regions, and the GO clustering function analysis is performed on these genes to analyze whether the differentially related genes have obvious functional clustering. That is, whether to perform certain types of functions in a targeted manner. 6. Insert fragment length statistics, which can be done in the following ways:

 Compare the length of the reads on the pair to determine whether the comparison result matches the length of the built-in fragment, and count the percentage of the total number of reads in the length of the different inserts.

 7. Coverage of CpG sites, the analysis can be performed in the following ways:

 In the RRBS, the fragments were obtained by restriction enzyme digestion and then subjected to BS sequencing. The theoretical fragment of the restriction enzyme digestion (40-220 bp) was counted in the chromosomal spots and their proportions in each chromosome and genome, and then counted. The coverage of the actual data obtained by the two samples at the theoretical site.

 The invention also provides a device for bioinformatics analysis based on heavy bisulfite hydrogenation sequencing or reduced representative bisulfite sequencing, the device comprising bisulfite sequencing and/or reduced representation a bisulfite sequencing module for performing bisulfite sequencing or reduced representative bisulfite sequencing of the sample, and preferably filtering the sequencing sequence to remove unqualified sequences; data output analysis module, Statistics are made on the output of the sequencing data, including the conversion rate, the length of the sequencing fragment, the number of sequencing fragments, and the total amount of the sample; the sequencing fragment comparison module is used to compare the fragments obtained by sequencing to the reference. On the genome, and compare the results, the statistical content includes the original number of original fragments, the number of comparison fragments, the amount of comparison data, the average comparison ratio, the average coverage depth; the coverage analysis module of the thiolation site Used to calculate the ratio of the number of C-sites covered by the thiolated sequencing sequence to the total number of sites at different depths of coverage; Point methylation level analysis module for counting the average thiolation level of each thiolated C base; genome-wide thiolation data distribution trend analysis module, used to obtain the whole genome thiolation distribution trend, statistics The number and distribution ratio of CG, CHG and CHH in the thiolated C base, where H stands for non-G 3⁄4; a genome-wide DNA thiolation map analysis module for describing the distribution of thiolated C bases from the chromosomal level Situation; a differential thiolation region analysis module for determining the thiolation difference region by the thiolation level of the two samples; the insert length statistical analysis module, comparing the length of the sequencing sequence on the pair to determine the ratio Whether the result conforms to the length of the built-in fragment; the coverage analysis module of the CpG locus is used to count the number of CG sites and their proportions covered in each chromosome and genome, and then statistically obtain the actual data at the theoretical site. Coverage on. DRAWINGS

 Figure 1. Percentage of loci at different depths of coverage based on BS sequencing.

 Figure 2. Degrees of loci at different depths of coverage based on BS sequencing.

 Figure 3. Distribution ratios of different types of thiolated C-based groups based on BS sequencing.

Figure 4. The thiolation level distribution of thiolated C based on BS sequencing. Figure 5. Quaternized C-attached features of non-CG sites based on BS sequencing.

 Figure 6. Illustratively shows the thiolated C density distribution on chromosomes based on BS sequencing.

 Figure 7. C 1⁄2 曱-based distribution of different genomic regions based on BS sequencing, and the corresponding CpG density.

 Figure 8. The horizontal distribution of thiolation in different functional component regions of the entire base group based on BS sequencing. Figure 9. Insert length statistics based on RRBS sequencing.

 Figure 10. Percentage of loci at different depths of coverage based on RRBS sequencing.

 Figure 11. Coverage of covered CG loci in the genome based on RRBS sequencing.

 Figure 12. Coverage of CpG loci at different depths based on RRBS sequencing. detailed description

 In the present invention, bisulfite sequencing or reduced representative bisulfite sequencing can be performed by high throughput sequencing techniques, such as Illumina GA sequencing technology, or other existing high throughput sequencing techniques. Removal of unacceptable sequences can be performed by those skilled in the art based on the actual conditions of sequencing, as needed. For example, the unqualified sequence includes: The number of bases whose sequencing quality is lower than the preset threshold exceeds 50% of the entire sequence and is considered to be a failed sequence. The low-quality threshold is determined by the specific sequencing technology and sequencing environment; the number of subtractions (such as N in Illumina GA sequencing results) whose sequence is uncertain in sequencing is more than 10% of the number of bases in the entire sequence is considered to be inconsistent. Each sequence; in addition to the sample linker sequence, is aligned with other exogenous sequences introduced by the experiment, such as various linker sequences. If there is a foreign sequence in the sequence, it is considered to be a non-conforming sequence.

 In the present invention, the short sequence mapping program may be a block program such as SOAP (available from http://soap.genomics.org.cn/).

 In the present invention, the screening results of the thiolation can be selected by, for example, unique comparison, comparison of comparison lengths, comparison of mismatches, comparison of the number of comparisons, etc., to screen out the alignment results of each sequence. For comparison information, the selection criteria should be based on the selected software and sequence background.

 In the present invention, reads refers to a sequence of fragments obtained from a sequencer, or a read.

In the present invention, the number of minimum expected sequencing sequences of the site thiolation is determined according to the following inequality, ie, the minimum k value satisfying the following formula: leg C * Q) p ^k l - p) ⁿ ~ ^k < 0.01 * mC; Its towel,

k is the expected value to support the number of singulated reads, which reflects the minimum support at a certain error rate曱 Number of base sequencing sequences,

 mC is the number of sequencing sequences in which the first sequencing sequence or the second sequencing sequence covering the site is thiolated at the corresponding position,

 nmC is the number of sequencing sequences that cover the first sequencing sequence or the second sequencing sequence of the site that are not deuterated at the corresponding position,

 n is the total number of sequencing sequences at this site,

 p is the error rate, for example, the above steps 1) - 2) are carried out using unmodified control DNA such as lamda DNA, and the p value is determined according to the conversion rate of C to T tt, p is usually 0.001-0.01, The number of minimum expected sequencing sequences that determine the thiolation of the site is determined based on the binomial distribution. The present invention will be described in detail below with reference to the accompanying drawings and embodiments. invention. Example 1 Bio-information analysis method based on BS sequencing Sample: Sample of two people

Record

First, data preprocessing

 1. Filter the sequencing sequence;

 In this example, two human blood sample DNAs numbered NB and 103 were subjected to bisulfite sequencing or reduced representative bisulfite sequencing at Shenzhen Huada Gene Research Institute, wherein the sequencing technology was Illumina GA sequencing technology. . Sequencing yielded 657.2 trillion and 660.5 trillion sequences, respectively, and then the sequencing sequence was filtered to remove the unqualified sequences, yielding 575.46 megabytes and 576.12 megabytes of qualified sequences, respectively. The unqualified sequence includes: The number of sequencing quality below the preset threshold (66 in this example) exceeds 50% of the entire sequence minus the number of bases. Bases with undefined sequencing results in the sequence (in this example, N in Illumina GA sequencing results) are 10% of the entire number of sequences; exogenous sequences introduced by other experiments except sample linker sequences ; There is an exogenous sequence in the sequence.

 2. Determine the thiolation information represented by the above-described filtered sequencing sequence

The reference genomic sequence of the species used in this example is the human genome hgl8, from the database UCSC (http://genome.ucsc.edu/). The filtered sequencing sequence and the reference sequence are preprocessed as follows: i) a sequence obtained by converting all C (cytosine) bases in the respective sequencing sequences into T (thymidine) bases is defined as fql; all G (guanine) bases in the complementary sequence of the respective sequencing sequences The sequence obtained by conversion to A (adenine) tt is defined as fq2, and the schematic is as follows: fql

Iq2

In the figure, C ^m represents methylated C (cytosine)

 Ii) converting the sequence generated by converting all c (cytosine) in the reference genome sequence to τ (thymidine) tt as refl; converting all G (guanine) 1⁄2 in the reference genome sequence into A (adenine) The sequence generated is defined as re 2;

 Iii) aligning each sequence fql and fq2 with reference sequences refl and ref2, respectively, and obtaining positional information of each sequence fql and fq2 on the transformed reference genomic sequence, which is retained on the transformed reference genome sequence. The unique sequences fql and fq2 are as follows:

Refl fql

 M M

 i;ef2 iv) combining the sequence information of a sequence fql and fq2 on the transformed reference genome sequence with the sequence information at the corresponding position to determine the thiolation information of the sequences fql and fq2: if in the sequence Both fql and C (cytosine) bases at the corresponding positions in the reference sequence indicate that the position of the sequence fql is thiolated; if the sequence fq2 and the corresponding position in the reference sequence are G (bird)嘌呤) base, which means that the position of the sequence fq2 is methylated, the schematic is as follows:

Fq2. T CC ACC C TO C S3⁄4 A^ 3⁄4CC Gft^AT CC GAAiSC TG G T AC AT CACj.

T TCCACCC TC CT A & C GSSG CC TGA 1⁄2C CCT AC AT C3⁄4G Sequence In the figure, * and # indicate undensylated, and C ^ra denotes thiolated C (cytosine). 3. After the above treatment, the sample NB obtained the 498.82 megabyte effective sequence, and the sample 103 obtained 49-1.28 megabytes with 3⁄4 ^ column, which can cover the entire genome to determine the thiolation site of the reference genome sequence.

For each locus on the reference genome, the number of sequencing sequences supporting the thiolation at this site, mC, and the number of sequencing sequences that do not support methylation, are counted. The following inequality is obtained from the binomial distribution to determine whether the site is methyl. Site: nmC * ( )p ^k (t - p) ⁿ - ^k < 0.01 * mC; where n is the total number of sequencing sequences at this site, p is the error rate (o.oi), and k is the supporting thiol group. The expected number of reads, mC is the number of sequencing sequences that support thiolation, and nmC is the number of sequencing sequences that do not support thiolation;

 The expected value k of the number of thiolated reads supported by each site on the reference genome sequence is obtained by the above formula. For a certain site, if the number of sequencing sequences supporting thiolation is greater than the expected value of the site, the site is a thiol group. The resulting site preferably forms a result file *.cout including the following information:

 Chr-a b +/- c d e f g h First column Chr-a: chromosome number,

 The second column b: the position of the C site,

 The third column +/-: positive and negative chain information,

 The fourth column c: the type of C site,

 The fifth column d: specific,

 Column 6 e: average alignment position of the sequencing sequence,

 Column 7 f: the number of sequencing sequences that cover the first or the second sequencing sequence of the C site and are thiolated at the corresponding position,

 The eighth column g: the number of sequencing sequences covering the first or the second sequencing sequence of the C site that are not thiolated at the corresponding position,

 Ninth column h: Number of minimum expected sequencing sequences for thiolation of the C site.

 The results are as follows:

 The NB methylation CG site of sample NB was 36.47 megabytes, accounting for 96.88% of the total CG site; the thiolated CHG site was 0.19 megabytes, accounting for 0.51% of the total CHG site; the thiolated CHH site was 0.98 megabytes. It accounts for 2.62% of the total CHH locus.

Sample 103 thiolated CG site was 35.01 megabytes, accounting for 97.56% of the total CG site; thiolated CHG site was 0.15 megabytes, accounting for 0.43% of the total CHG site; thiolated CHH site was 0.72 megabytes, Total 2.02% of the CHH locus.

 Second, biological information analysis

 Data output

 Statistics were made on the sequencing data of the DNA of the two jk liquid samples, including the conversion rate, the length of the sequencing fragment read, the number of sequencing fragments, and the total amount of data in the sample. See Table 1 for the specific results.

 2. Sequencing fragment alignment

 The sequenced fragments are aligned to the genome of the species. For specific alignment methods, see Step 2 of Data Preprocessing: Determination of the position of the sequence on the genome. After the data comparison is completed, the comparison results of the two human samples are counted, and the statistics include Raw reads, Raw bases, Mapped reads, and comparisons t. See Figure 2 for the results of the Mapped bases, Average mapped rate, and Average coverage depth (X), the average number of times the entire genome is covered by reads.

 3. Coverage of thiolated sites

 C 1⁄2 sites can be divided into three types, namely CG, CHG, CHH, where H represents a non-G base. From the *.cout file we can get whether the site is a thiolated site, and the site that supports the thiolized reads greater than the expected value is the thiolation site. The ratio of mC locus to the total locus under different coverage depths (the number of times the locus is covered by reads) is counted. The statistical result is shown in Figure 1. The horizontal axis represents greater than or equal to the coverage depth, and the vertical axis represents the corresponding depth. proportion. At the same time, the genomic coverage of the sites under different coverage depths is also counted, and the coverage is the proportion of the number of sites on the full genomics. The statistical results are shown in Fig. 2, the horizontal axis ^ the depth of coverage, and the vertical axis represents the coverage of the site under the corresponding depth.

 4. The thiolation level of the thiolated site

 To understand the overall characteristics of the genomic thiolation map and the average thiolation level of the whole genome, the average thiolation level of each methylated C base is calculated by:

 The number of reads that support merging / (the number of reads that support merging + the number of reads that do not support merging) xl'00. The specific results are shown in Table 3.

 At the same time, the statistics reflect the reads of the priming information by 3⁄4^ degrees, and the coverage calculation is based on three aspects: chromosomes, gene intervals, and gene functional elements. After the reads with ambiguous state are excluded, the effective coverage calculation is performed based on the read-supported reads and the reads data that does not support the support, and the statistical results are obtained, as shown in Table 4.

 5. Whole genome methylation data distribution trend

In order to obtain the trend of genome-wide thiolation distribution, statistics of thiolated C bases in CG, CHG and The number and distribution ratio of CHH, that is, the percentage of each type of methylated C to the total thiolation C, for example: the proportion of mCHG = the number of mCHG / the total number of mC. The statistical results are shown in Table 5 and Figure 3.

 At the same time, different distribution types of C CG, CHG and CHH ) have different levels of thiolation between different species and even different cell types of the same species. Therefore, the distribution of the thiolation of each type (CG, CHG, and CHH) is statistically calculated. The statistical results are shown in Figure 4: where the horizontal axis represents the thiolation level from low to high, for every 10 % is a file, and the vertical axis represents the proportion of C 3⁄4 ^ all thiolation C of a particular thiolation level. Similarly, the thiolation levels of C bases (CG, CHG, and CHH) for each distribution type in each region are calculated by chromosome, gene interval, and gene function elements, and the formula is: Average thiolation level = support thiolation The read/ (supports simplistic reads+ does not support simplified reads) xl00. The statistical results are shown in Table 6.

 Finally, the distribution of bases near the thiolated C of non-CG sites was analyzed, and the probability of occurrence of different thiolation patterns (CHG and CHH) was statistically analyzed. The results are shown in Fig. 5.

6. Whole genome DNA thiolation map

 The distribution of thiolated C-velocines is described from the chromosomal level. The distribution of density on the chromosome was counted by a 10 kb window; at the same time, the density distribution of different types of thiolated C minus groups (CG, CHG, and CHH) in each lOOkb window was counted, and then plotted by drawing software. Figure, the results are shown in Figure 6.

 Different genomic regions have different biological functions, and the methylation distribution in the form of heat maps helps to further understand the thiolated graphical features of these regions.

 The results are shown in Figure 7: where each graph represents the heat map distribution characteristics of a genomic region, where n represents the number of CpGs involved in the analysis (each CpG has a coverage depth greater than 4X on a single strand). The horizontal axis represents the number of CpGs per 200 bp window (ie, CpG density), and the vertical axis represents the C 1⁄2 average thiolation level of CpG. The black thin line in the figure indicates the median of the average thiolation level of c bases in this window at a specific CpG density. From the shallow to the deep, the B region indicates the number of CpGs at a specific thiolation level and CpG density. The upper A area represents the distribution of CpG density, mapped to the horizontal axis, and the right C area represents the distribution of the thiolation level and is mapped to the vertical axis.

The transcription of genes is regulated by different regulatory elements, dividing all coding gene sequences into seven different regions of transcriptional elements, here! ^ Statistics on the level of thiolation of different transcriptional element regions. The distribution of DNA thiolation levels in different functional regions helps to understand the role of DNA thiolation modifications in different regions from the genome-wide level. The statistical results are shown in Fig. 8: the horizontal axis represents 7 different transcription element regions, and the vertical axis represents the average thiolation level of each point in a specific region, and each red dot is a bin (person For the average thiolization level of a defined length interval, the curve represents the average 5-bin shift, and the dashed line between a and b is the TSS (transcription start site).

 7. Differential methylation region (DMR) analysis

 The methylation difference region was determined by the thiolation level of the two samples as follows: First, a CG site was identified as a seed site, and two samples were obtained to support thiolation and thiolization at the site, respectively. Four values were used to check whether the thiolation level of the two samples was significant at the site by a chi-square test; if the difference was significant, a CG site was extended backward, and two samples were supported at these two sites.曱 化 和 和 和 和 和 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个 四个Area. This site is then used as a seed site to cycle the above process, resulting in a differential thiolated region.

 After finding the differential thiolated regions, the corresponding genes, namely the thiolated differential genes, are found through these regions, and the GO clustering function analysis is performed on these genes to analyze whether the differentially related genes have obvious functional clustering. That is, whether to perform certain types of functions in a targeted manner. Example results:

 Table 1, data output statistics for two samples

Table 2, comparison results statistics

表 3 , 不同类型曱基化 C威基的平均甲基化水平

Table 3, average methylation levels of different types of thiolated C-Wiki

Table 4, the readings that reflect the singularization information

Table 5, Number and proportion of different distribution types thiolated C

The level of thiolation of various types of C in different gene regions

实施例二， 基于 RRBS测序的生物信息分析方法 样本： d液样本 DNA

Example 2: Samples of biological information analysis methods based on RRBS sequencing: d liquid sample DNA

 Have a relationship

 First, data preprocessing

1. Filter the sequencing sequence; In this example, the liquid sample DNA numbered YH is subjected to heavy subhydro acid hydrogenation sequencing or reduced representative heavy acid sulfite sequencing in the Shenzhen Hua λ ^ Institute, wherein the sequencing technology is Illumina GA sequencing technology. . Sequencing yielded 94.5 mega-sequences, and the sequencing sequence was then filtered to remove unqualified sequences, resulting in 83.62 mega-qualified sequences. The unqualified sequence includes: The number of bases whose sequencing quality is lower than the preset threshold (66 in this example) exceeds 50% of the number of bases minus the entire sequence. The number of subtractions (in the present example, the number of N in the Illumina GA sequencing results) in the sequence is more than 10% of the entire sequence; in addition to the sample linker sequence, the exogenous sequence introduced by other experiments; The presence of an exogenous sequence in the sequence is considered a non-conforming sequence.

 2. Determine the methylation information represented by the above filtered sequencing sequence

The reference genomic sequence of the species used in this embodiment is the human genome hgl8, from the database UCSC. The filtered sequencing sequence and the species reference sequence are pretreated as follows: i) the sequence obtained by converting all C (cytosine) tt in the respective sequencing sequences into T (thymine) bases is defined as Fql; The sequence obtained by converting all G (guanine) bases in the complementary sequence of each sequencing sequence into A (adenine) ^^ is defined as fq2, and the schematic is as follows: fql ΘΘΤ bond βΆΘΘΑΤΤΜΤ霞(3⁄4'ΤΪ职3⁄43⁄4$ TT 3⁄43⁄4S Beach STAGI^ ββ T TeC3⁄4 CCTC Job TO嚷"GCT _: Cft3⁄43⁄4TCG: The simplified surface is shown in the figure, C ^m represents the thiolated C (cytosine)

 Ii) converting all C (cytosine) 1⁄2 in the reference genomic sequence into T (thymine) to generate a sequence defined as refl; converting all G (guanine) subunits in the reference genome sequence to A (gland)嘌呤) The generated sequence is defined as ref2;

 Iii) aligning each sequence fql and fq2 with reference sequences refl and ref2, respectively, and obtaining positional information of each sequence fql and fq2 on the transformed reference genomic sequence, which is retained on the transformed reference genome sequence. The unique sequences fql and fq2 are as follows:

_ef2 T iv )将一个序列 fql和 fq2在转化后的参考基因组序列上的位置信息参考序 列在相应位置处的序列信息相结合， 确定该序列 fql和 fq2的甲基化信息： 如 果在该序列 fql上和在参考序列相应位置上都是 C (胞嘧啶）碱基， 则表示在 该序列 fql的该位置被曱基化； 如果在该序列 fq2上和在参考序列相应位置都 是 G 则表示在该序列 fq2的该位置被曱基化示意图如下：

_r

_e f2 T Iv) combining the sequence information of a sequence fql and fq2 on the transformed reference genome sequence with the sequence information at the corresponding position to determine the methylation information of the sequences fql and fq2: if and in the sequence fql If the C (cytosine) base is at the corresponding position of the reference sequence, it means that the position of the sequence fql is thiolated; if both the sequence fq2 and the corresponding position in the reference sequence are G, the sequence is represented. The position of fq2 is singulated as follows:

, 权面 i^eG^TceT:!^卿: 雄听3⁄4?3⁄4:3⁄433⁄43⁄4 3⁄43⁄4|0<≤^3⁄43 ^: 3⁄4 listen AT:e 2: τ^Αε^'Θ σϊ^^ΐ^Ε^^^ Βόδ^^τ, ίρί^^^τΐ^^漏;^ [" Xiang: In the figure, * and # indicate undensified, and C ^m denotes thiolated C (cytosine) 1⁄2. After the above treatment, the sample YH obtained 63.45 megabytes with 3⁄4 ^ columns, covering 6.2% of the CG sites of the entire genome.

 Since RRBS sequencing is performed by BS sequencing of the digested fragments, we need to filter the obtained effective sequences by cleavage sites, and then obtain different cleavage points according to different alcohols. Here, Mspl is used to identify CCGG sites. Enzyme. Both ends of the effective sequence are judged whether there is a restriction site for Mspl. If one of the ends has the alcohol cleavage site, the effective sequence is considered to be a normal restriction fragment, which is used for subsequent analysis.

 3. Identify the thiolation site on the reference genome sequence

For each locus on the reference genome, the number of sequencing sequences supporting the thiolation at this site, mC, and the number of sequencing sequences that do not support thiolation, are counted. According to the binomial distribution, the following inequality is used to determine whether the site is a sulfhydryl group. Site: nmC * ()p ^k (l - v ~ ^k < 0.01 * mC; where n is the total number of sequencing sequences at this site, p is the error rate (0.01 in this case), k is the supporting thiol group The expected number of reads, mC is the number of sequencing sequences that support thiolation, and nmC is the number of sequencing sequences that do not support thiolation;

The expected value k of the number of thiolated reads supported by each site on the reference genome sequence is obtained by the above formula. For a certain site, if the number of sequencing sequences supporting thiolation is greater than the expected value of the site, the site is a thiol group. The resulting site preferably forms a result file *.cout including the following information: Chr-a b +/- cdefgh First column Chr-a: chromosome number,

 The second column b: the position of the C site,

 The third column +/-: positive and negative chain information,

 The fourth column c: the type of C site,

 The fifth column d: the specific ^,

 Column 6 e: average alignment position of the sequencing sequence,

 Column 7 f: the number of sequencing sequences that are ligated at the corresponding position by the first sequencing sequence or the second sequencing sequence covering the C site,

 Column 8 g: the number of sequencing sequences covering the first or the second sequencing sequence of the C site that are not thiolated at the corresponding position,

 Ninth column h: Number of minimum expected sequencing sequences for thiolation of the C site.

The results are as follows:

 The YH methylation CG site of the sample was 1.84 megabytes, accounting for 52.93% of the covered CG sites; the thiolated CHG site was 0.013 megagrams, accounting for 0.49% of the total CHG sites; the thiolated CHH site was 0.036. Mega, accounting for 2.55% of the total CHH locus.

 Second, biological information analysis

 Data output

 The sequencing data of the DNA of the two jk liquid samples were generated by counting the libraries of different enzyme fragments, and the statistics included the length of the sequencing fragments, the number of sequencing fragments, and the total amount of samples. See Table 7 for specific results.

 2. Sequencing fragment alignment

 The sequenced fragments are aligned to the genome of the species, and the specific comparison method is described.

 According to step 2 of the pretreatment: the determination of the position of the sequence on the genome. After the data comparison is completed, the comparison results of the two human samples are counted, and the statistics include the Raw reads, the Raw Data, the Mapped reads, and the comparisons. See Table 8 for specific results of Mapped bases, Mapped rate, and Enzyme cutting rate.

 3. Insert length statistics

Compare the length of the reads on the pair to determine whether the comparison results match the length of the database. Calculate the percentage of reads in the total number of reads for different insert lengths. The statistical results are shown in Figure 9. 4. Coverage of thiolated sites

 The C site can be divided into three types, namely CG, CHG, CHH, where H stands for non-G. From the *.cout file we can get whether the site is a thiolated site, and the site that supports the thiolized reads greater than the expected value is the thiolation site. The ratio of the mC locus to the total locus under different coverage depths (the number of times the locus is covered by the reads) is counted. The statistical result is shown in Fig. 10. The horizontal axis represents greater than or equal to the coverage depth, and the vertical axis represents the locus at the corresponding depth. proportion. At the same time, the coverage of the sites in different coverage depths in each region of the genome is also counted, and the coverage is the proportion of the number of sites covered in each region of the genome. The statistical results are shown in Fig. 11, the horizontal axis ^ the depth of coverage, and the vertical axis represents the coverage of the site at the corresponding depth. Finally, the average coverage depth of CG loci in each genomic region was counted. See Table 9 for the results. In addition, the average coverage depth of CG loci in each region of the statistical genome

5. Coverage of CpG sites

 RRBS was digested to obtain fragments and then subjected to BS sequencing. Here, the theoretical fragment (40-220 bp) of the enzyme digestion was counted in the number of CG sites and their proportions in each chromosome and genome, and then the two samples were actually counted. The coverage of the obtained data at the theoretical site is shown in Table 10 and Figure 12.

 6. Differential thiolized region (DMR) analysis

 The thiolation difference region was determined by the thiolation level of the two samples as follows:

 First, a CG locus is determined as a seed site, and two samples are supported at the site to support the thiolation and the four values that do not support thiolation, and then the chi-square level of the two samples is tested by the chi-square test. Whether the difference in the locus is significant; if the difference is significant, a CG locus is extended backward, and two samples are supported at these two loci to support the thiolation and the four values that do not support thiolation, and then the chi-square test is performed. Such a push cycle, when encountering a point where the difference is no longer significant, <more stop, the front significant area is the differential thiolated area. This site is then used as a seed site to cycle the above process to obtain a differential thiolated region.

 After finding the differential thiolated regions, the corresponding genes, namely methylation differential genes, are found through these regions, and the GO clustering function analysis is performed on these genes to analyze whether the differentially related genes have obvious functional clustering. That is, whether to perform certain types of functions in a targeted manner.

The RRBS is subjected to restriction enzyme digestion to obtain fragments and then subjected to BS sequencing, and the results include, but are not limited to, the above. RRBS can also be analyzed as appropriate for the content of bioinformatics analysis of BS sequencing. Example results:

 Table 7, Data output statistics for two samples

Table 8, statistics of comparison results

Table 9. Average coverage depth of CG loci in each genomic region

 CG average 3⁄4 degree (X) gene component

 LB-NB9 LB-NB10

Upstream of CG Island 2k 18.04 17.61

CG Island 20.92 19.89

Downstream of CG Island 2k 17.68 17.39 Upstream 2k 18.7 17.33

5-untranslated region 21.17 19.95 coding sequence 19.98 20.67 intron 18.02 17.55

3-untranslated region 17.69 17.73 downstream 2k 20.11 19.85 promoter 19.94 18.73 ii due 20.78 20.98 tandem repeat 18.46 17.59

Alu repeat sequence 13.49 11.3 ribosomal RNA 982.78 726.59 transfer RNA 20.45 18.81 0, coverage of CpG sites

Claims

Claim  CLAIMS 1. A method for detecting DNA thiolation by heavy bisulfite hydrogenation sequencing or reduced representative heavy acid salt sequencing, the method comprising the steps of:  1) subjecting the sample to bisulfite sequencing or reduced representative bisulfite sequencing, and preferably filtering the sequencing sequence to remove unqualified sequences;  2) Determine the thiolation site represented by the above sequencing sequence, the steps are as follows  i) a first sequencing sequence obtained by converting all C (cytosine) in the respective sequencing sequences into T (thymine) bases; converting all G (guanine) bases in the complementary sequences of the respective sequencing sequences into A (adenine) obtained second sequencing sequence;  ϋ ) converting all C (cytosine) 1⁄2 in the reference genome sequence of the corresponding species into the first reference sequence generated by the Τ (thymidine) subtraction; converting all G (guanine) bases in the reference genome sequence of the corresponding species a second reference sequence generated for Α (adenine);  Iii) aligning each of the first sequencing sequence and the second sequencing sequence with the first reference sequence and the second reference sequence, respectively, and obtaining each of the first sequencing sequence and the second sequencing sequence in the first reference sequence and the second Selecting, in the position on the sequence, each of the first sequencing sequence and the second sequencing sequence having unique alignments on the first reference sequence and the second reference sequence;  Iv) determining the methylation information of the first sequencing sequence and the second sequencing sequence by using the sequence information of the first sequencing sequence and the second sequencing sequence aligned by the second sequencing sequence, as follows: If in the first sequencing sequence and The corresponding positions of the aligned reference genomic sequences are C (cytosine) bases, indicating that they are thiolated at this position of the first sequencing sequence; if at the corresponding position of the second sequencing sequence and the aligned reference genome sequence All of which are G (guanine) bases, indicating that they are thiolated at this position of the second sequencing sequence; 3) for each C site on the reference genomic sequence, the number of sequencing sequences that are thiolated at the corresponding position of the first or second sequencing sequence covering the site and the number of sequencing sequences that are not thiolated, The minimum number of expected sequencing sequences for thiolation of the site is determined according to the following inequality, ie, the minimum k satisfying the following formula: nmC * ( )p ^k (tp ~ ^k < 0.01 * mC; wherein  mC is the number of sequencing sequences in which the first sequencing sequence or the second sequencing sequence covering the C site is thiolated at the corresponding position, nmC is the first sequencing sequence covering the C site or the second sequencing sequence is not thiol at the corresponding position Number of sequencing sequences,  n is the total number of sequencing sequences at this site,  p is 4 Wu rate,  If the number of sequencing sequences that are methylated at the corresponding position by the first or second sequencing sequence covering the C site is greater than the minimum number of desired sequencing sequences at the site, the C site is the site of thiolation .  2. The method of claim 1, wherein the step 1) -2) is carried out by replacing the sample with an unmodified control DNA such as lamda DNA, and the p value is usually determined according to the conversion ratio of C to T minus, p is usually It is 0.001-0.01, such as 0.01.  3. The method of claim 1 or 2, further comprising forming a result file comprising: a chromosome number, a position of the C site, positive and negative strand information, a type of C site, a specific base, and a coverage of the C site The number of the first sequencing sequence or the second sequencing sequence of the point, the number of sequencing sequences that are thiolated at the corresponding position by the first sequencing sequence or the second sequencing sequence covering the C site, and the first number covering the C site The number of sequencing sequences that the sequencing sequence or the second sequencing sequence is not thiolated at the corresponding position and/or the minimum number of desired sequencing sequences that are thiolated by the C site.  4. A method of bioinformatics analysis based on bisulfite sequencing or reduced representative bisulfite sequencing, the method comprising detecting DNA thiolation by the method of any of claims 1-3, Bioinformatics analysis is then performed, which analyzes one or more selected from the group consisting of: data production, sequencing of fragments, coverage of thiolated sites, thiolation sites The level of thiolation, genome-wide thiolation data distribution trends, genome-wide DNA thiolation maps, differential thiolated regions, insert length statistics, and coverage of CpG sites.  5. The method of claim 4, wherein the coverage analysis of the thiolation site is performed as follows: C tt sites are classified into three types, namely CG, CHG, CHH, wherein H represents a non-G base; The result file obtains whether the site is a thiolation site, and the site supporting the thiolated reads larger than the expected value is the thiolation site; the ratio of the mC site to the total site at different coverage depths is also calculated; For sites with different coverage depths on the genome, ί3⁄4Α is counted, and the coverage is the proportion of the number of sites on the whole genome. 6. The method of claim 4, wherein the thiolation level analysis of the thiolation site is performed as follows: To understand the overall characteristics of the genomic thiolation map and the average methylation level of the whole genome, each methyl group is counted. The average thiolation level of C bases; statistically reflects the effective coverage of reads of thiolated information. The calculation of coverage is based on three aspects: chromosome, gene interval and gene function After the components are excluded from the ambiguous reads, the effective coverage calculation is performed based on the read-supported singly-based reads and the reads data that does not support methylation, and the statistical results are obtained.  7. The method of claim 4, wherein said genome-wide thiolation data distribution trend analysis is performed as follows: To obtain a genome-wide thiolation distribution trend, statistical thiolation C 3⁄4 CG, CHG and The number and distribution ratio of CHH, ie the percentage of each type of thiolated C to the total thiolation C; the distribution of the C thiolation level of each type (CG, CHG and CHH) thiolation; Gene intervals and gene functional elements count C CG, CHG and each type of distribution in each region The thiolation level of CHH); or the distribution of bases near the thiolated C at non-CG sites, and the probability of occurrence of different thiolation patterns (CHG and CHH).  8. The method of claim 4, wherein said genome-wide DNA thiolation pattern analysis is performed as follows: Delineating the distribution of methylated C bases from the chromosomal level, and counting the density of thiolation at the chromosome of 10 kb on the chromosome The distribution of the upper; at the same time, the density distribution of different types of thiolized C 1⁄2 ( CG, CHG and CHH ) in each window of 100 kb is counted, and different genomic regions have different biological functions, which are expressed in the form of heat maps. The basic distribution helps to further understand the thiolated graphical features of these regions. The transcription of genes is regulated by different regulatory elements, dividing all coding gene sequences into seven different transcriptional element regions, here! The statistics on the thiolation levels of different transcriptional element regions were made.  9. The method of claim 4, said differential thiolated region (DMR) analysis being performed as follows: determining the thiolation difference region by the thiolation level of the two samples, as follows: First determining a CG site As a seed site, two samples were obtained at the site to support thiolation and do not support thiolated four numbers, and then the chi-square test was used to check whether the thiolation level of the two samples was significant at the site; If the difference is significant, then a CG site is extended backwards, and two samples are supported at these two sites to support thiolation and do not support thiolated four numbers, and then perform chi-square test; When the point where the difference is no longer significant, it stops, and the area that is prominent in the front is the difference thiolized area. This site is used as a seed site to circulate the above process, thereby obtaining differential methylation regions. After finding the differential thiolated regions, the corresponding genes are found through these regions, namely the thiolated differential genes, and These genes were analyzed by GO clustering function to analyze whether the differentially related genes have obvious functional clustering, that is, whether to specifically control the exercise of certain functions.  10. The method of claim 4, wherein the statistical analysis of the length of the insert is performed as follows: comparing the length of the reads on the pair to determine whether the comparison result is consistent with the length of the database, and counting the length of the different inserts. The percentage of the number. 11. The method of claim 4, wherein the coverage analysis of the CpG sites is performed as follows: In RRBS, the fragment was obtained by enzyme digestion and then subjected to BS sequencing. Here, the theoretical fragment (40-220 bp) of the enzyme digestion was counted in the number of CG sites and their proportions in each chromosome and genome, and then the statistics were counted. The coverage of the actual data obtained from a sample at the theoretical site.  12. A device for bioinformatics analysis based on heavy bisulfite sequencing or reduced representative heavy sulphate acidation sequencing, the device comprising bisulfite sequencing and/or reduced representation a bisulfite sequencing module for performing bisulfite sequencing or reduced representative bisulfite sequencing of the sample, and preferably filtering the sequencing sequence to remove unqualified sequences; data output analysis module, Statistics are made on the output of the sequencing data, including the conversion rate, the length of the sequencing fragment, the number of sequencing fragments, and the total amount of data of the sample; the sequencing fragment alignment analysis module is used to compare the fragments obtained by sequencing to the reference. On the genome, and compare the results, the statistical content includes the initial fragment data volume, the number of comparison fragments, the amount of comparison data, the average alignment ratio, the average coverage depth; the coverage analysis module of the thiolated site It is used to calculate the ratio of the number of C-sites of the thiolated sequence of 3⁄4Α to the total number of sites at different depths of coverage; A horizontal analysis module for counting the average thiolation level of each methylated C tt; a genome-wide thiolation data distribution trend analysis module for obtaining a genome-wide thiolation distribution trend, statistical thiolation (中CG , the number and distribution ratio of CHG and CHH, where H represents non-G ^tt; a genome-wide DNA thiolation map analysis module for describing the distribution of methylated C bases from the chromosomal level; differential thiolation A region analysis module is used to determine the thiolation difference region by the methylation level of the two samples; the insert length statistical analysis module, and the length of the sequencing sequence on the comparison pair to determine whether the alignment result conforms to the database fragment Length; CpG site coverage analysis module, used to count the number of CG sites and their proportions covered in each chromosome and genome, and then statistically report the coverage of the actual data on the theoretical site.