A method for splitting and analyzing polyploid genome epigenetic data
By optimizing alignment parameters and screening strategies, the epigenetic data between homologous chromosomes can be accurately distinguished, solving the problem of inaccurate data alignment in traditional methods and achieving efficient splitting and analysis of homologous polyploid or genotypic genomes.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- WUHAN FRASERGEN CO LTD
- Filing Date
- 2025-09-30
- Publication Date
- 2026-07-17
AI Technical Summary
Traditional techniques struggle to accurately distinguish the epigenetic structure between homologous chromosomes in homologous polyploid or genotypic genomes, leading to data misalignment and low alignment rates.
After acquiring data using a high-throughput sequencing platform, the bowtie2 software was used to perform single-end alignment with the reference genome with the parameter -k'n'. Combined with specific screening rules and MACS software to identify peaks, polyploid subgenomes were split and independent reference genome indexes were constructed to improve the accuracy and efficiency of data alignment.
It enables accurate differentiation of epigenetic structures among homologous chromosomes in homologous polyploid or genotypic genomes, improves data alignment rate, and can also retrieve peaks lost in traditional methods, improving the data alignment rate by 3%.
Smart Images
Figure CN121354668B_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of bioinformatics technology, specifically relating to a method for splitting and analyzing epigenetic data of polyploid genomes. Background Technology
[0002] With the development of genome sequencing and assembly technology, more and more homologous polyploids and homologous chromosome typing genomes are being assembled. Analyzing the regulation of gene expression by epigenetic structural differences among homologous chromosomes of different species is of great significance in the study of heterosis and complex genetic diseases.
[0003] Traditional techniques such as ATAC-seq and CUT&TAG directly align reference genomes using alignment software like bowtie2 and BWA. Due to the influence of non-global alignment algorithms, they often fail to accurately distinguish signals between homologous chromosomes, resulting in data errors and low effective data rates. Therefore, there is an urgent need to provide a new method based on ATAC-seq and CUT&TAG techniques to improve the accurate differentiation of epigenetic structures between homologous chromosomes in homologous polyploid or genotypic genomes, and to increase the data alignment rate. Summary of the Invention
[0004] In view of this, the present invention provides a method for splitting and analyzing epigenetic data of polyploid genomes, which can accurately distinguish epigenetic data between homologous chromosomes, achieve accurate differentiation of epigenetic structures between homologous chromosomes of homologous polyploid or genotypic genomes, and improve the data comparison rate.
[0005] To achieve the above objectives, the present invention adopts the following technical solution:
[0006] A method for splitting and analyzing polyploid genome epigenetic data includes the following steps:
[0007] S1. Use a high-throughput sequencing platform to acquire data, remove data adapters and filter low-quality sequences. The data includes ATAC-seq data, ChIP-seq data, and CUT&TAG data.
[0008] S2. Use bowtie2 software with -k'n' parameter to align one end of the reference genome. Filter the alignment results. If only one of the multiple alignment sites in the Reads does not mismatch, record the read ID and the chromosome number of the matching reference genome.
[0009] S3. Record the paired-end matching Reads and the matching chromosome numbers, according to the following rules:
[0010] If both ends of the reads are marked, and the mismatched sites at both ends are on the same chromosome, record the read ID and chromosome number;
[0011] If a single-end read is marked and the other end has an ID that matches the same chromosome record, record the read ID and chromosome number.
[0012] S4. Divide the reference genome into multiple sets according to the polyploid subgenome or diploid genome, and construct a separate reference genome index for each set of genomes;
[0013] The polyploid subgenome includes homologous or heterologous tetraploids, homologous hexaploids, and homologous octaploids;
[0014] Among them, 2, 3, and 4 sets of reference genomes were constructed from homologous or heterologous tetraploids, homologous hexaploids, and homologous octoploids, respectively.
[0015] Two reference genomes were constructed for diploid genotyping.
[0016] S5. Extract the sequence from the record obtained in step S3 and assign it to the corresponding haploid genome according to the chromosome number;
[0017] S6. Align each sequence to the corresponding haploid genome according to the bowtie / bowtie2 / bwa pairwise alignment process. Finally, use MACS software to identify the peaks information of each genome, thus completing the polyploid genome epigenetic data splitting analysis.
[0018] Furthermore, in step S1, the connector removal and filtering of low-quality sequences includes the use of Trimmomatic, CutAdapt, SOAPnuke, and FastP quality control tools.
[0019] Furthermore, in step S2, n in parameter -k'n' represents the genomic ploidy or the number of alleles for the target split.
[0020] In some specific embodiments, specifically when the target to be analyzed is a diploid typing reference genome, n=2 in parameter -k'n';
[0021] When the target to be analyzed is a homologous tetraploid reference genome, n=2 in parameter -k'n';
[0022] When the target of analysis is a tetraploid genotyping reference genome, n=4 in parameter -k'n';
[0023] When the target to be analyzed is a homologous hexaploid reference genome, n=3 in parameter -k'n';
[0024] When the target to be analyzed is a homologous octoploid reference genome, n=4 in parameter -k'n'.
[0025] In some specific embodiments, preferably, the specific screening process in step S2 is as follows:
[0026] Parse the SAM file's 6th column, Sigar information, to calculate the alignment region length $len; sum the numbers before "M" and "D"; parse the MD:Z: field sequence alignment mismatch information $num; and count the total number of "ATCG" fields after "MD:Z:".
[0027] Calculate the mismatch rate: $sn_rate=$num / $len. A read has at most 2 matching results. A record has only one matching result $sn_rate=0 for reads.
[0028] In some specific embodiments, preferably, the reference genome index is constructed separately in step S4 using bowtie2-build, as shown in the following code:
[0029] bowtie2-build --threads 12 refA.fa refA;
[0030] bowtie2-build --threads 12 refB.fa refB.
[0031] In some specific embodiments, the identification of each set of genome peaks information in step S6 is as follows:
[0032] When the data is ATAC-seq data, use the MACS software with the parameters --shift -75 --extsize 150 to identify peaks;
[0033] When the data is ChIP-seq data, use the MACS software to identify peaks in either dual-ended (-f BAMPE, -f BEDPE) or single-ended (-f BED, -f BAM) mode.
[0034] When the data is CUT & TAG data, use the MACS software to identify peaks in either dual-ended (-f BAMPE, -f BEDPE) or single-ended (-f BED, -f BAM) mode.
[0035] Furthermore, MACS software also includes MACS2 and MACS3 software.
[0036] Compared with the prior art, the beneficial effects of the present invention are as follows:
[0037] This invention optimizes alignment parameters and incorporates subsequent sequence identification and screening strategies to accurately distinguish epigenetic data between homologous chromosomes, achieving accurate differentiation of epigenetic structures between homologous chromosomes in homologous polyploid or genotypic genomes. Peaks partially lost in traditional analysis methods can be retrieved using this approach, and the data alignment rate is improved by 3%. Attached Figure Description
[0038] Figure 1 This is a flowchart of the technology of the present invention.
[0039] Figure 2 This is a comparison chart of the results of Embodiment 1 and Comparative Example 1 of the present invention. Detailed Implementation
[0040] The present invention will be further described in detail below with reference to specific embodiments, so that those skilled in the art can more clearly understand the present invention. Unless otherwise specified, the technical means used in the following embodiments are all conventional means well known to those skilled in the art, and all reagents and consumables are commercially available products.
[0041] Example 1
[0042] This embodiment provides a method for splitting and analyzing epigenetic data of polyploid genomes. The specific steps are as follows (see detailed flowchart). Figure 1 ):
[0043] 1. Data Filtering
[0044] For the raw sequencing data of tetraploid plants: MS_1_ATAC.fq.gz; MS_2_ATAC.fq.gz, Trimmomatic was used to remove sequence adapters and filter low-quality sequences. The code is as follows:
[0045] java -Djava.io.tmpdir=. / Trimmomatic_tmp / -jar trimmomatic-0.39.jar PE-threads 8 -summary trim.stats -validatePairs MS_1.fq.gz MS_2.fq.gz;
[0046] ILLUMINACLIP:BGI-SEQ-PE.fa:2:30:10:8:true LEADING:3;
[0047] TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:50 -baseout MS_ATAC_clean.fastq.gz.
[0048] The resulting files are MS_ATAC_clean_1P.fastq.gz and MS_ATAC_clean_1P.fastq.gz.
[0049] 2. Reference genome index construction
[0050] Use bowtie2-build to build a reference genome index. The code is as follows:
[0051] bowtie2-build --threads 12 ref.fa ref.
[0052] 3. Data Comparison
[0053] The bowtie2 software was used to align one end of the reference genome with the -k2 parameter (meaning to retain the information of the two best-matching sites). The code is as follows:
[0054] bowtie2 -x ref -U MS_ATAC_clean_1P.fastq.gz -S out_ATAC_1P.sam -p 8 -k 2;
[0055] bowtie2 -x ref -U MS_ATAC_clean_1P.fastq.gz -S out_ATAC_2P.sam -p 8 -k 2.
[0056] 4. Compare the screening results
[0057] The comparison results are filtered and recorded according to the following rule: if only one of the maximum two alignment sites of a read does not mismatch, the read ID and the matching reference genome chromosome number are recorded. The specific method is as follows:
[0058] Parse the SAM file, column 6, Sigar information, calculate the alignment region length $len: sum the numbers before "M" and "D":
[0059] Parse the MD:Z: field sequence for mismatch information. $num: Count the total number of "ATCG" fields following "MD:Z:".
[0060] Calculate the mismatch rate: $sn_rate=$num / $len.
[0061] A read can have up to two alignment results. Reads with only one alignment result $sn_rate=0 are marked.
[0062] 5. Two-end matching results
[0063] The paired-end matching results are used to filter out the retained sequence IDs and chromosome numbers. The specific method is as follows:
[0064] A: Filter reads with labels at both ends, and then filter reads whose chromosome numbers match at the $sn_rate=0 site in the alignment results at both ends.
[0065] B: Filter reads with labels at one end and check the alignment results at the other end. If the chromosome number that matches the $sn_rate=0 position in the alignment results at the other end is retained.
[0066] Reads that satisfy condition A or B are retained sequences.
[0067] 6. Disassemble the reference genome
[0068] The reference genome was split into two parts based on the tetraploid subgenome information, and a separate index was built using bowtie2-build. The code for the index is as follows:
[0069] bowtie2-build --threads 12 refA.fa refA;
[0070] bowtie2-build --threads 12 refB.fa refB.
[0071] 7. Conventional method for paired-end reference genome alignment
[0072] Extract the corresponding subgenome sequences according to the reference genome numbers of the preserved reads: MS_1_ATAC_genomeA.fq.gz; MS_2_ATAC_genomeA.fq.gz and MS_1_ATAC_genomeB.fq.gz; MS_2_ATAC_genomeB.fq.gz.
[0073] The two sets of sequences were aligned with their corresponding reference genomes. The code is as follows:
[0074] bowtie2 -1 MS_1_ATAC_genomeA.fq.gz -2 MS_2_ATAC_genomeA.fq.gz -p 8 -xrefA -X 2000 --dovetail 2>Agenome.summary | samtools view -F 4 -u -| \samtools sort -@ 6 -o ATAC_genomeA.bam -;
[0075] bowtie2 -1 MS_1_ATAC_genomeB.fq.gz -2 MS_2_ATAC_genomeB.fq.gz -p 8 -xrefB -X 2000 --dovetail 2>Bgenome.summary | samtools view -F 4 -u -| \samtools sort -@ 6 -o ATAC_genomeB.bam -.
[0076] 8. Conventional methods for identifying peaks
[0077] Use MACS2 to identify subgenomic peak information. The code is as follows:
[0078] samtools sort --threads 8 -nATAC_genomeA.bam | bedtools bamtobed -bedpe -mate1 -i / dev / stdin | gzip -nc>ATAC_genomeA.bedpe.gz;
[0079] zcat -f ATAC_genomeA.bedpe.gz | awk '{printf $1"\t"$2"\t"$3"\tN\t1000\t"$9"\n", $4"\t"$5"\t"$6"\tN1000\t"$10"\n"}' | gzip -nc>ATAC_genomeA.bed.gz;
[0080] macs2 callpeak -f BED -t ATAC_genomeA.bed.gz -n sheep -g 659709883 --shift -75 --extsize 150 --nomodel -B --SPMR --keep-dup all --call-summits;
[0081] samtools sort --threads 8 -nATAC_genomeB.bam | bedtools bamtobed -bedpe -mate1 -i / dev / stdin | gzip -nc>ATAC_genomeB.bedpe.gz;
[0082] zcat -f ATAC_genomeB.bedpe.gz | awk '{printf $1"\t"$2"\t"$3"\tN\t1000\t"$9"\n", $4"\t"$5"\t"$6"\tN1000\t"$10"\n"}' | gzip -nc>ATAC_genomeB.bed.gz;
[0083] macs2 callpeak -f BED -t ATAC_genomeB.bed.gz -n sheep -g 662033183 --shift -75 --extsize 150 --nomodel -B --SPMR --keep-dup all --call-summits.
[0084] For details, please see [link / details]. Figure 2 .
[0085] Comparative Example 1
[0086] This comparative example uses the traditional ATAC-seq analysis method for data splitting, and the specific steps are shown in the flowchart (). Figure 1 The difference between this example and Example 1 is obvious. Specifically, steps 3, 4, 5, and 6 are omitted; the rest are the same. The final result is shown below. Figure 2 .
[0087] Furthermore, in order to understand the differences between the splitting results of the embodiments and the comparative examples, the specific data between the two were compared, and the specific comparison is shown in Table 1.
[0088] Table 1. Comparison details of data from each embodiment and comparative example.
[0089]
[0090] Combined Table 1 and Figure 2 It can be seen that this method, by optimizing alignment parameters and adding subsequent sequence identification and screening strategies, accurately distinguishes epigenetic data between homologous chromosomes, achieving accurate differentiation of epigenetic structures between homologous chromosomes in homologous polyploid or genotypic genomes. Figure 2 It can be seen that some of the missing peaks in Comparative Example 1 can be obtained in Example 1, and the data comparison rate is improved by 3%.
[0091] Unless otherwise specified, all raw materials used in this invention are existing substances that can be purchased directly from the market.
[0092] The above are merely preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.
Claims
1. A method for splitting and analyzing polyploid genome epigenetic data, characterized in that, Includes the following steps: S1. Use a high-throughput sequencing platform to acquire data, remove data adapters and filter low-quality sequences. The data includes ATAC-seq data, ChIP-seq data, and CUT&TAG data. S2. Use bowtie2 software with -k'n' parameter to align one end of the reference genome. Filter the alignment results. If only one of the multiple alignment sites in the Reads does not mismatch, record the read ID and the chromosome number of the matching reference genome. S3. Record the paired-end matching Reads and the matching chromosome numbers, according to the following rules: If both ends of the Reads are marked, and the mismatched sites at both ends are on the same chromosome, record the readsID and chromosome number; If a single-end read is marked and the other end has an ID that matches the same chromosome record, record the read ID and chromosome number; S4. Divide the reference genome into multiple sets according to the polyploid subgenome or diploid genome, and construct a separate reference genome index for each set of genomes; The polyploid subgenome includes homologous or heterologous tetraploids, homologous hexaploids, and homologous octaploids; Among them, 2, 3, and 4 sets of reference genomes were constructed from homologous or heterologous tetraploids, homologous hexaploids, and homologous octoploids, respectively. Two reference genomes were constructed for diploid genotyping. S5. Extract the sequence from the record obtained in step S3 and assign it to the corresponding haplotype genome according to the chromosome number; S6. Align each sequence to the corresponding haplotype genome according to the bowtie / bowtie2 / bwa pairwise alignment process. Finally, use MACS software to identify the peaks information of each genome, thus completing the polyploid genome epigenetic data splitting analysis.
2. The method according to claim 1, characterized in that, Step S1 involves removing and filtering low-quality sequences by using tools such as Trimmomatic, CutAdapt, SOAPnuke, and FastP.
3. The method according to claim 1, characterized in that, In step S2, n in parameter -k'n' represents the genomic ploidy or the number of alleles for the target split.
4. The method according to claim 3, characterized in that, When the target of analysis is a diploid typing reference genome, n=2 in parameter -k'n'; When the target to be analyzed is a homologous tetraploid reference genome, n=2 in parameter -k'n'; When the target of analysis is a tetraploid genotyping reference genome, n=4 in parameter -k'n'; When the target to be analyzed is a homologous hexaploid reference genome, n=3 in parameter -k'n'; When the target to be analyzed is a homologous octoploid reference genome, n=4 in parameter -k'n'.
5. The method according to claim 4, characterized in that, The specific filtering process in step S2 is as follows: Parse the SAM file's 6th column, Sigar information, to calculate the alignment region length $len; sum the numbers before "M" and "D"; parse the MD:Z: field sequence alignment mismatch information $num; and count the total number of "ATCG" fields after "MD:Z:". Calculate the mismatch rate: $sn_rate=$num / $len. A read has at most 2 matching results. A record has only one matching result $sn_rate=0 for reads.
6. The method according to claim 1, characterized in that, In step S4, the reference genome index is constructed separately using bowtie2-build, as shown in the code below: bowtie2-build --threads 12 refA.fa refA; bowtie2-build --threads 12 refB.fa refB.
7. The method according to claim 1, characterized in that, The specific steps for identifying peak information for each genome in step S6 are as follows: When the data is ATAC-seq data, use the MACS software with the parameters --shift -75 --extsize 150 to identify peaks; When the data is ChIP-seq data, use the MACS software to identify peaks in either dual-ended (-f BAMPE, -f BEDPE) or single-ended (-f BED, -f BAM) mode. When the data is CUT & TAG data, use the MACS software to identify peaks in either dual-ended (-f BAMPE, -f BEDPE) or single-ended (-f BED, -f BAM) mode.
8. The method according to claim 7, characterized in that, MACS software also includes MACS2 and MACS3 software.