Transcriptome sequencing-based screening method for differentially expressed genes in muscovy ducks and application thereof

Abstract

The application discloses a kind of based on transcriptome sequencing's muscovy duck difference gene screening method and application, belong to gene detection technical field, eukaryotic transcriptome sequencing technology is used to the analysis of muscovy duck breast muscle sample, realizes difference gene screening by data quality control processing, genome alignment and expression amount standardization, difference gene preliminary screening, multidimensional verification and core difference gene determination.This application has formulated strict and adapted transcriptome sequencing data quality control standard for muscovy duck breast muscle sample, clear original data fuzzy base N%, GC content, Q20 / Q30 proportion and clean data Clean Reads / Clean Data proportion such as quantitative threshold value, while ensuring that data is unbiased and reaches saturation threshold by multidimensional verification such as gene coverage uniformity, saturation analysis.

Description

Technical Field

[0001] This invention relates to the field of gene detection technology, specifically to a method and application for screening differential genes in Muscovy ducks based on transcriptome sequencing. Background Technology

[0002] Muscovy ducks are an important livestock and poultry breed. Research on the genetic improvement and molecular regulation mechanism of their pectoral muscle-related traits is an important direction in the field of livestock and poultry breeding. Transcriptome sequencing technology, as a core means to analyze gene expression regulation and discover functional genes, has key application value in Muscovy duck molecular breeding research.

[0003] Eukaryotic transcriptome sequencing with reference is a routine technique for studying gene expression characteristics in species. The standard procedure involves filtering the raw sequencing data to obtain high-quality, clean data, then aligning this clean data to a reference genome for subsequent studies such as gene expression level calculation, differential expression analysis, and functional enrichment analysis. This technology has already been applied in livestock and poultry research. Data alignment with the reference genome can be performed using HISAT2 software, and precise gene expression level statistics can be achieved using HTSeq software combined with FPKM normalization. Differentially expressed genes are screened using Deseq2 software with |log2FoldChange| > 1 and P-value < 0.05 as criteria. The biological functions and signaling pathways involved in these differentially expressed genes are then analyzed using GO, KEGG enrichment analysis, and GSEA analysis.

[0004] Meanwhile, transcriptome sequencing analysis can be further used for research such as transcription factor prediction, transcript structure analysis, and gene variation detection. Transcription factor families can be identified by comparison with the AnimalTFDB database, transcript splicing and new transcript mining can be completed using StringTie software, differential alternative splicing events such as skipped exons and alternative splicing sites can be identified using rMATS software, and gene variations such as SNPs and InDel can be detected and functionally annotated using Varscan program combined with ANNOVAR annotation. This provides technical support for the discovery of key genes, molecular markers and regulatory elements related to livestock and poultry traits.

[0005] However, in transcriptome studies related to Muscovy duck breast muscles, systematic analysis of gene expression characteristics, differentially expressed gene profiles, key functional genes, and regulatory pathways for specific Muscovy duck samples still needs further development. Based on standardized eukaryotic reference transcriptome sequencing analysis procedures, the discovery of differentially expressed genes, core regulatory pathways, and molecular markers related to Muscovy duck breast muscle traits is of great practical significance for the selection and genetic improvement of superior Muscovy duck breeds and is also an urgent problem to be solved in current Muscovy duck molecular breeding research.

[0006] Based on this, the present invention designs a method and application for screening differential genes in Muscovy ducks based on transcriptome sequencing to solve the above problems. Summary of the Invention

[0007] To address the aforementioned shortcomings of existing technologies, this invention provides a method and application for screening differentially expressed genes in Muscovy ducks based on transcriptome sequencing. The method employs eukaryotic transcriptome sequencing technology to analyze Muscovy duck breast muscle samples. Differentially expressed gene screening is achieved through data quality control processing, genome alignment and expression level standardization, initial screening of differentially expressed genes, multi-dimensional verification, and identification of core differentially expressed genes. Specifically, the method includes:

[0008] The original data shows that the percentage of fuzzy bases (N%) is 0.0081-0.0366%, the GC content is 48.10-49.39%, the percentage of Q30 bases is ≥95.0%, and the percentage of Q20 bases is ≥98.0%.

[0009] Clean data refers to Clean Reads accounting for ≥98.0% and Clean Data accounting for ≥97.8%. Clean Reads refers to the number of high-quality sequencing sequences after filtering, and Clean Data refers to the total number of bases in the high-quality sequencing sequences after filtering.

[0010] The differential gene screening threshold is |log2FoldChange|>1 and P-value<0.05.

[0011] Furthermore, the method includes the following steps:

[0012] S1: Quality control processing of raw transcriptome sequencing data from Muscovy duck breast muscle samples;

[0013] The raw sequencing data in FASTQ format was filtered using Fastp software or equivalent transcriptome data filtering software to remove sequences with adapters at the 3′ end and reads with an average quality score lower than Q20, thus obtaining clean data.

[0014] The clean data were sequentially subjected to base quality distribution detection, base content distribution detection, average quality distribution of reads detection, gene coverage uniformity analysis and saturation analysis. The clean data were confirmed to be unbiased and the sequencing data volume reached the saturation threshold by the gene coverage uniformity curve showing no single peak bias and the relative error of gene expression level in the saturation analysis being ≤5%.

[0015] S2: Genome alignment and gene expression level normalization analysis of clean data;

[0016] Using HISAT2 software or a functionally equivalent genome alignment software, clean data that passed quality control were aligned to the Muscovy duck reference genome. Default parameters were used for non-strand-specific libraries, while the --rna-strandness RF parameter was used for the first strand library and the --rna-strandness FR parameter was used for the second strand library. The overall alignment rate between the clean data and the reference genome was 59.11%–67.69%.

[0017] The original gene expression level was counted using the Union statistical scheme of HTSeq software or a functionally equivalent expression level statistical method. When constructing a chain-specific library, it is necessary to determine whether the Read and Feature directions in the annotation are consistent. The original gene expression level was converted into the number of fragments per thousand bases of transcript per million aligned fragments using the FPKM method to eliminate the influence of gene length and sequencing depth.

[0018] Pearson correlation coefficient was used to analyze the correlation of gene expression levels among samples, and samples with correlation coefficients in the range of 0.8-1.0 were screened. Principal component analysis (PCA) was used to verify the effectiveness of sample grouping and clustering, ensuring that samples in the same group clustered together and samples in different groups were clearly separated.

[0019] S3: Initial screening of differentially expressed genes in Muscovy ducks;

[0020] Deseq2 software or functionally equivalent differential expression analysis software was used to perform inter-group differential expression analysis on the original expression levels of genes that passed S2 validation. With |log2FoldChange|>1 and P-value<0.05 as the screening thresholds, differentially expressed genes in Muscovy ducks were obtained. The total number of differentially expressed genes in different comparison groups was 298-400, of which 151-208 were upregulated and 147-246 were downregulated.

[0021] Collect common and unique genes among different comparison groups;

[0022] S4: Multidimensional verification of differential genes in Muscovy ducks;

[0023] Protein-protein network interaction analysis was performed on differentially expressed genes using the STRING database, and protein-protein interaction pairs with a score > 0.95 were selected. If the STRING database did not have Muscovy duck PPI information, chicken, turkey, or Muscovy duck species with close kinship were selected. The BLAST algorithm was used to compare the Muscovy duck protein sequence with the protein sequence of the selected species, and homologous protein interaction relationships with similarity ≥ 80% were selected and mapped back to Muscovy duck as Muscovy duck protein-protein interaction pairs with a score > 0.95.

[0024] ClusterProfiler software or equivalent enrichment analysis software were used to perform GO enrichment analysis and KEGG enrichment analysis. Significantly enriched GO terms and KEGG pathways were screened with a P-value < 0.05. The GO enrichment analysis analyzed the functions of differentially expressed genes from three dimensions: molecular function, cellular components, and biological processes. The KEGG enrichment analysis completed KO annotation and metabolic pathway annotation.

[0025] GSEA analysis was conducted, using |NES|>1, NOM p-val<0.05, and FDR q-val<0.25 as thresholds to screen genes with insignificant differential expression but important biological significance.

[0026] S5: Identification of core differentially expressed genes in Muscovy ducks;

[0027] To further characterize the regulatory features of differentially expressed genes, at least one of the following analyses will be performed:

[0028] By comparing with the AnimalTFDB database, we can identify transcription factors and their families in differentially expressed genes.

[0029] The rMATS software was used to identify five types of alternative splicing events: exon skipping, variable 5′ splicing site, variable 3′ splicing site, mutually exclusive exons, and retained introns. Differential alternative splicing events were screened based on IncLevelDiff≥0.1 and FDR<0.05.

[0030] SNP and InDel detection of differentially expressed genes were performed using Varscan. The filtering criteria were Q > 20, number of reads for covered sites > 8, number of reads for supporting mutation sites > 2, and P-value < 0.01. Functional annotation of variant sites was performed using ANNOVAR software. The annotation scope included intron regions, intergenic regions, exon regions, 5′UTR regions, 3′UTR regions, splice regions, upstream regions, downstream regions, and ncRNA exon regions.

[0031] Based on the differential screening results of S3, the functional enrichment and protein interaction results of S4, and the above analysis information, genes that simultaneously meet the criteria of significant differential expression, high-confidence protein interaction relationship, and enrichment in functional pathways related to Muscovy duck pectoral muscle traits are prioritized as core differential genes of Muscovy duck.

[0032] Genes that are significantly differentially expressed and enriched in key functional pathways but for which transcription factor regulation or gene variation has not been detected can also be included in the core differentially expressed gene set, taking into account their biological function and their core position in the pathway.

[0033] Furthermore, the HISAT2 software described in S2 simultaneously calculates the proportion of multiple alignment sequences and the proportion of unique alignment sequences, and the distribution of alignment regions in the CDS region, introns, and intergenic regions conforms to the characteristics of eukaryotic transcriptome sequencing.

[0034] Furthermore, the rMATS software described in S5 uses a junction counts quantitative method, through a formula... The expression level of alternative splicing is calculated, where I is the number of transcript reads containing the alternative splicing region, LI is the effective length of the corresponding transcript, S is the number of transcript reads that skip the alternative splicing region, and LS is its effective length.

[0035] Furthermore, the core pathways enriched by KEGG described in S4 include the PPAR signaling pathway, the adrenergic signaling pathway, the estrogen signaling pathway, the cholesterol metabolism pathway, and the fatty acid degradation pathway.

[0036] Furthermore, the saturation analysis described in S1 uses RSeQC software to resample clean data at ratios of 5%, 10%, 15%...100%. The expression levels of all genes are calculated 20 times at each sampling ratio. The relative error is calculated by comparing the gene expression levels at different sampling ratios with the actual expression levels at 100% sampling. All genes are divided into four groups, Q1-Q4, according to their expression levels. When the relative error of all four gene groups is ≤5%, the data volume is determined to have reached the saturation threshold.

[0037] Furthermore, the PCA principal component analysis described in S2 retains the first two principal components, whose cumulative contribution rate is ≥70%.

[0038] Furthermore, the protein network interaction analysis results described in S4 can be plotted using Cytoscape software or a functionally equivalent visualization tool to present the gene interaction network topology.

[0039] An application of the Muscovy duck differential gene screening method based on transcriptome sequencing is disclosed, which is applied to Muscovy duck molecular breeding, genetic improvement of Muscovy duck breast muscle traits and mining of Muscovy duck functional genes. The application is based on the core differential genes of Muscovy ducks screened by the method.

[0040] Compared with the prior art, the beneficial effects of this invention are as follows:

[0041] 1. This invention establishes stringent and suitable quality control standards for transcriptome sequencing data from Muscovy duck breast muscle samples. It specifies quantitative thresholds for raw data, including the percentage of ambiguous bases (N%), GC content, Q20 / Q30 ratio, and the ratio of clean reads to clean data. Simultaneously, it ensures data unbiasedness and reaches saturation thresholds through multi-dimensional verification, such as gene coverage uniformity and saturation analysis. The analysis workflow standardizes parameters and methods for genome alignment, expression level standardization, and differential screening, unifying the differential gene screening threshold for |log2FoldChange|>1 and P-value<0.05, significantly improving the accuracy of differential gene screening and the reliability of experimental results.

[0042] 2. This invention breaks through the limitations of traditional transcriptome research, which relies solely on differential expression analysis and simple functional enrichment to screen genes. It utilizes protein network interaction analysis (PPI), GO / KEGG enrichment analysis, and GSEA analysis to perform multi-dimensional functional verification of differentially expressed genes. Simultaneously, it combines transcription factor mining, differential alternative splicing analysis, and SNP / InDel variant detection with multi-level regulatory feature analysis to establish a standard for identifying core differentially expressed genes that are significantly different, have high-confidence interactions, are enriched in core pathways, and have clearly defined regulatory characteristics. This allows for the precise identification of core differentially expressed genes related to traits such as Muscovy duck pectoral muscle development and lipid metabolism, effectively excluding interference from biologically insignificant differentially expressed genes, and providing a precise candidate gene set for Muscovy duck functional gene mining. Attached Figure Description

[0043] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the accompanying drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are merely some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without any creative effort.

[0044] Figure 1 Line plot of base quality distribution for clean transcriptome data of Muscovy duck breast muscle samples;

[0045] Figure 2 A graph showing gene saturation analysis at different expression levels in the clean transcriptome data of Muscovy duck breast muscle samples.

[0046] Figure 3 Genome alignment of clean transcriptome data from Muscovy duck breast muscle samples; Uniform distribution of gene coverage.

[0047] Figure 4 Principal component analysis scatter plot of transcriptome gene expression levels in Muscovy duck breast muscle samples;

[0048] Figure 5 Volcano plot of differentially expressed genes in transcriptome among Muscovy duck pectoral muscle sample groups;

[0049] Figure 6 Clustering heatmap of differential gene expression patterns in transcriptome of Muscovy duck pectoral muscle samples;

[0050] Figure 7 Bar chart for KEGG pathway enrichment analysis of differentially expressed genes in transcriptome of Muscovy duck breast muscle samples;

[0051] Figure 8 This is a topology diagram of the differentially expressed protein-protein interaction (PPI) network in the transcriptome of Muscovy duck breast muscle samples. Detailed Implementation

[0052] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0053] Example 1: This example provides a method for screening differentially expressed genes in Muscovy ducks based on transcriptome sequencing. Eukaryotic transcriptome sequencing with reference was used to analyze Muscovy duck breast muscle samples. Differentially expressed genes were screened through standardized data processing, genome alignment, expression level analysis, and multi-dimensional validation. Specifically, the method includes:

[0054] The original data showed that the percentage of fuzzy bases (N%) ranged from 0.0081% to 0.0366%, and the GC content ranged from 48.10% to 49.39%.

[0055] Clean Reads accounted for ≥98.0% and Clean Data accounted for ≥97.8%;

[0056] The threshold for differential gene screening was |log2FoldChange|>1 and P-value<0.05, which specifically included data quality control processing, genome alignment and expression level standardization, initial screening of differential genes, and multi-dimensional verification steps.

[0057] The method includes the following steps:

[0058] S1: Quality control processing of raw transcriptome sequencing data from Muscovy duck breast muscle samples;

[0059] We selected pectoral muscle tissue samples from healthy Muscovy ducks, extracted RNA, constructed libraries, and then performed eukaryotic transcriptome sequencing. The sequencing platform's built-in software converted the generated image files into raw sequencing data in FASTQ format. The raw data included core statistical information such as sequencing data volume, Q20 base percentage, Q30 base percentage, N% percentage of fuzzy bases, and GC content.

[0060] Key quality indicators of the raw sequencing data were pre-validated, including: the percentage of ambiguous bases (N%) controlled within the range of 0.0081-0.0366%, the GC content controlled within the range of 48.10-49.39%, the percentage of Q30 bases ≥95.0%, and the percentage of Q20 bases ≥98.0%.

[0061] A standardized filtering process is performed on the raw sequencing data using Fastp software or a functionally equivalent transcriptome data filtering software, specifically including:

[0062] (1) Remove the sequence carrying the sequencing adapter at the 3′ end;

[0063] (2) Remove Reads with an average quality score below Q20 (i.e., sequences with a high proportion of bases and a Phred score < 20), and obtain high-quality clean data after filtering.

[0064] Subsequently, a multi-dimensional quality verification of the cleaned data was conducted, with the following specific steps:

[0065] (1) Detect the Q value distribution of each base position (5′→3′ direction) in the clean data to confirm that the sequencing error rate conforms to the inherent characteristics of high-throughput sequencing technology, that is, the error rate gradually increases with the increase of sequencing sequence length, and the high error rate caused by the incomplete binding of random primers to RNA template during the library construction of the first 6 base genes is a normal fluctuation, with no abnormal sudden high error rate region.

[0066] (2) Detect the coverage distribution of sequencing reads from the 5′ end to the 3′ end of the gene (divided according to the base length ratio of 0-100%). The coverage curve verifies that there is no obvious single peak bias, that is, the reads are uniformly distributed in different regions of the gene (5′ end, middle region, 3′ end) without significant bias coverage, ensuring the accuracy of gene expression level calculation.

[0067] (3) Saturation analysis was performed using RSeQC software. The specific process is as follows:

[0068] The clean data was resampled at gradient ratios of 5%, 10%, 15%...100%. The expression levels of all genes were calculated 20 times at each sampling ratio. The gene expression levels at different sampling ratios were compared with the actual expression levels under 100% sampling (i.e., complete clean data) and the relative error was calculated.

[0069] When the relative error of all genes (divided into four groups Q1-Q4 according to their expression levels, where Q1 represents the lowest 25% of expression and Q4 represents the highest 25% of expression) is ≤5%, the sequencing data volume is considered to have reached the saturation threshold. At this point, the data volume is sufficient to support the accurate calculation of gene expression levels, and there is no need to increase the sequencing depth.

[0070] After the above multi-dimensional verification, the cleaning data must meet the following quality standards:

[0071] Clean Reads account for ≥98.0%, and Clean Data accounts for ≥97.8%, ensuring the integrity and reliability of cleaning data.

[0072] S2: Genome alignment and gene expression level normalization analysis of clean data;

[0073] HISAT2 software or a functionally equivalent genome alignment software was used to align clean data with the Muscovy duck reference genome using the improved BWT algorithm. This algorithm has the advantages of fast alignment speed and low resource consumption. Default parameters were used for non-strand-specific libraries, while parameters were precisely specified according to the type of strand-specific libraries—the --rna-strandness RF parameter was used for the first strand library, and the --rna-strandness FR parameter was used for the second strand library to ensure alignment specificity.

[0074] After alignment is completed, the total alignment rate (Total Mapped), which is the proportion of sequences aligned with the reference genome to the total number of clean data sequences, should be controlled within the range of 59.11-67.69%.

[0075] The Union statistical scheme of HTSeq software or a functionally equivalent expression level statistical method is used to count the Read Count value aligned to each gene as the original expression level. If the library is constructed for strand specificity, it is necessary to additionally determine whether the Read direction is consistent with the Feature direction in the annotation to ensure statistical accuracy.

[0076] Since the Read Count value is directly related to gene length and sequencing depth, it needs to be standardized using the FPKM method. Specifically, the original gene expression level is converted into the number of fragments per thousand bases of transcript per million aligned fragments, which completely eliminates the interference of gene length and sequencing depth and enables direct and comparable analysis of expression levels between different genes and different samples.

[0077] The specific steps for validating the sample validity are as follows:

[0078] (1) Correlation analysis between samples: Pearson correlation coefficient was used to assess the similarity of gene expression levels between samples, and samples with correlation coefficients in the range of 0.8-1.0 were screened.

[0079] (2) Principal component analysis (PCA) was performed on the standardized expression data to reduce the dimensionality of the samples and to verify the clustering effect of the samples, ensuring that samples in the same group were clustered together and samples in different groups were clearly separated.

[0080] S3: Initial screening of differentially expressed genes in Muscovy ducks;

[0081] The Deseq2 software or a functionally equivalent differential expression analysis software was used to perform inter-group differential expression analysis on the raw expression levels (Raw Counts) of genes obtained in step S2. The normalization algorithm built into the DESeq2 software automatically eliminates the influence of sequencing depth and library composition, without the need to pre-input FPKM values. During the analysis, |log2FoldChange|>1 was used as the fold change threshold and P-value<0.05 was used as the statistical significance threshold. The dual thresholds were used to screen for differentially expressed genes in Muscovy ducks. The core statistical indicator after screening was that the total number of differentially expressed genes in different comparison groups was 298-400, of which the number of upregulated differentially expressed genes was 151-208 and the number of downregulated differentially expressed genes was 147-246.

[0082] At the same time, a volcano plot was drawn for visual verification. The volcano plot used log2 (FoldChange) as the horizontal axis and -log10 (pvalue) as the vertical axis. Significantly different regions were defined by double dashed lines, and the top 10 genes that were upregulated or downregulated were marked.

[0083] Using the Pheatmap or ComplexHeatmap software packages in R, bidirectional clustering analysis was performed on the differential gene unions and samples of all comparison groups. During the clustering process, the Euclidean method was used to calculate the distances between genes and samples, and the hierarchical clustering longest distance method (Complete Linkage) was used for clustering.

[0084] The clustering results are presented as a heatmap, with the horizontal axis representing differentially expressed genes and the vertical axis representing samples. The color intensity reflects the relative level of gene expression (red indicates high expression, and blue indicates low expression).

[0085] The left side of the heatmap is the sample clustering tree, where samples with similar expression patterns are clustered together; the right side is the gene clustering tree, where genes with similar functions or expression patterns are grouped together.

[0086] Based on the differential expression analysis results, common and unique genes among different comparison groups were counted.

[0087] S4: Multidimensional verification of differential genes in Muscovy ducks;

[0088] Based on the STRING database, we can predict the protein-protein interaction relationships of differentially expressed genes. This database comprehensively includes protein-protein interaction information from experimental verification, data mining, and homology prediction, which can meet the needs of Muscovy duck differential gene interaction analysis.

[0089] When Muscovy duck PPI information is already included in the database, PPI interaction pairs containing differentially expressed genes and with an interaction confidence score > 0.95 are directly selected. If Muscovy duck-specific PPI information is not available in the database, poultry species closely related to Muscovy ducks, such as chickens and turkeys, are selected. The BLAST algorithm is used to perform homology comparison between Muscovy duck protein sequences and the protein sequences of the selected species. Homologous protein interaction relationships with sequence similarity ≥ 80% are selected and mapped back to Muscovy ducks as high-confidence Muscovy duck protein interaction pairs with a score > 0.95. If the PPI network size obtained after mapping is too large or too small, the network structure can be optimized by adjusting the score threshold. Finally, a result file containing all target gene interaction relationships is generated. This can be visualized using Cytoscape software or the Paisenno Gene Cloud Platform to present the gene interaction network topology and identify core interacting genes.

[0090] Differential gene function enrichment analysis was conducted using clusterProfiler software or functionally equivalent enrichment analysis software based on the Gene Ontology (GO) database. The GO database is an internationally standardized gene function classification system that covers three core dimensions: molecular function, cellular component, and biological process. Its basic unit is the Term (a unique identifier consisting of "GO:" followed by 7 numbers). The Term of various Ontologies forms a directed acyclic topology structure through is_a, part_of, and regulate relationships, which can systematically describe the attribute characteristics of genes and gene products.

[0091] During the analysis, all annotated genes obtained from sequencing were used as the background gene set. The number of differentially expressed genes in each GO Term and the number of background genes corresponding to that Term were counted. The P-value was calculated using the hypergeometric distribution method to assess whether the GO Term was significantly enriched in the list of differentially expressed genes (i.e., differentially expressed genes tend to be concentrated in this Term relative to the background gene set). Significantly enriched GO Term was screened with a P-value < 0.05 as the threshold to clarify the main biological functions performed by differentially expressed genes.

[0092] The enrichment results can be visualized using a bar chart with -log10 (P-value) as the horizontal axis, showing the degree of enrichment of significantly enriched GO terms.

[0093] Differential gene metabolic pathway enrichment analysis was conducted using clusterProfiler software or equivalent enrichment analysis software, based on the Kyoto Encyclopedia of Genes and Genomes.

[0094] KEGG annotations contain two core elements:

[0095] First, there is KO (KEGG Ortholog) annotation, which enables cross-species annotation of molecular networks;

[0096] Second, KEGG Pathway annotation constructs a network of intermolecular interactions and responses within a species, which can clarify the metabolic and signal transduction pathways involved by differentially expressed genes; the analysis workflow is as follows:

[0097] Using all sequenced annotated genes as the background gene set, the number of differentially expressed genes in each KEGG Pathway and the number of background genes corresponding to that Pathway were counted. The P-value was calculated using the hypergeometric distribution method, and KEGG Pathways with significant enrichment were screened with a P-value < 0.05. The core metabolic and signaling pathways involved by the differentially expressed genes (such as the PPAR signaling pathway, adrenergic signaling pathway, estrogen signaling pathway, cholesterol metabolism pathway, etc.) were identified. The enrichment results were visualized using a bar chart with -log10(P-value) as the x-axis to show the enrichment degree of each Pathway.

[0098] S5: Identification of core differentially expressed genes in Muscovy ducks;

[0099] To further characterize the regulatory features of differentially expressed genes and identify Muscovy duck differentially expressed genes with core biological functions, one or more of the following analyses can be performed: transcription factor analysis, differential alternative splicing analysis, and SNP and InDel variant detection. The core differentially expressed genes in Muscovy ducks can then be determined by combining the results of the above analyses. The specific procedure is as follows:

[0100] Using the differentially expressed genes of Muscovy ducks obtained from step S3 as the research object, their gene sequences were compared with the AnimalTFDB animal transcription factor database for homology. Gene sequences belonging to transcription factors were screened out from the differentially expressed genes through sequence matching. At the same time, the family type of each transcription factor (such as zf-C2H2, Homeobox, bHLH, CSRNP_N, etc.) was determined.

[0101] The rMATS software was used to identify and screen differentially expressed genes for differentially expressed alternative splicing events. The software defaults to Junction Counts quantification, using only Reads crossing splice sites for quantitative analysis, and the results were obtained using the formula... Calculate the expression level of alternative splicing, where I is the number of transcript reads containing the alternative splicing region, LI is the effective length of the corresponding transcript, S is the number of transcript reads that skip the alternative splicing region, and LS is its effective length;

[0102] The software can accurately identify five types of alternative splicing events: exon skipping (SE), variable 5′ splice site (A5SS), variable 3′ splice site (A3SS), mutually exclusive exons (MXE), and retained introns (RI). It uses IncLevelDiff≥0.1 and FDR<0.05 as dual screening criteria to screen for significantly differentially expressed genes.

[0103] Varscan was used to detect single nucleotide polymorphisms (SNPs) and insertion / deletion mutations (InDels) in differentially expressed genes. To ensure the reliability of the variant sites, a four-fold filtering standard was set:

[0104] (1) The Q value of the base at the variant site is >20;

[0105] (2) The number of reads covering this site is >8;

[0106] (3) The number of Reads supporting the mutation site is >2;

[0107] (4) The P-value of the mutation site is <0.01. After filtering, the effective SNPs and InDel variant sites in the differentially expressed genes are obtained.

[0108] The above-mentioned effective variant sites were functionally annotated using ANNOVAR software. The annotation scope covered intron regions, intergenic regions, exon regions, 5′UTR regions, 3′UTR regions, splice regions, upstream regions, downstream regions, and ncRNA exon regions.

[0109] Based on the differential expression screening results of S3 (|log2FoldChange|>1, P-value<0.05), and integrating the protein network interaction (PPI Score>0.95), GO enrichment analysis (P-value<0.05), KEGG enrichment analysis (P-value<0.05), and GSEA analysis (|NES|>1, NOM p-val<0.05, FDR q-val<0.25) results of S4, a multi-dimensional comprehensive screening was conducted: genes that simultaneously meet the criteria of significant differential expression, high-confidence protein interaction relationship, and enrichment in functional pathways related to Muscovy duck pectoral muscle traits were prioritized as core differentially expressed genes in Muscovy ducks;

[0110] Genes that show significant differential expression and are enriched in key functional pathways but for which transcription factor regulation or gene variation has not been detected can also be included in the Muscovy duck core differential gene set, based on their biological function and core position in the pathway, and finally identify Muscovy duck core differential genes with important breeding value and molecular regulatory significance.

[0111] Example 2: Application of the transcriptome sequencing-based differential gene screening method for Muscovy ducks as described in Example 1 in Muscovy duck molecular breeding, genetic improvement of Muscovy duck breast muscle traits, mining of Muscovy duck functional genes, and regulation of Muscovy duck gene expression.

[0112] Example sample: Composed of 30 Muscovy duck breast muscle samples;

[0113] Table 1. Basic Information of Muscovy Duck Breast Muscat Muscovy Muscat Muscat Pectoral Mu

[0114] It should be noted that the biological duplicate numbers for G1-56, G2-56, G3-56, G4-56, and G5-56 in Table 1 correspond to 1, 2, 3, 4, and 5, respectively; the biological duplicate numbers for the remaining control groups (G1-70 to G5-70, G1-120 to G5-120) and experimental groups (M1-56 to M5-56, M1-70 to M5-70, M1-120 to M5-120) are sequentially assigned according to the suffix of the sample numbers (e.g., M1-56 corresponds to duplicate 1, M2-56 corresponds to duplicate 2, and so on).

[0115] The proportion of Q30 in sample G5-56 was 94.86%, slightly lower than the set threshold (≥95.0%). However, the percentage of fuzzy bases N% (0.036577%), GC content (48.74%), and Clean Reads percentage (98.25%) in this sample all met the standards. Furthermore, the data quality was verified to be qualified and without abnormal bias by saturation analysis (relative error ≤5%) and sample correlation analysis (correlation coefficient with the same group of samples ≥0.9959). Therefore, it can be included in subsequent analysis.

[0116] Table 2. Sample data of Muscovy duck breast muscle: ,

[0117] It should be noted that Table 2 is based on eukaryotic reference transcriptome sequencing and normalized analysis, as detailed below:

[0118] (1) The raw data indicators (number of raw reads, Q30 percentage, GC content, and fuzzy base N%) were obtained by converting the sequencing image files into FASTQ format raw data using the software provided by the sequencing platform, and then statistically analyzed by the Fastp software quality control module.

[0119] (2) The clean data indicators (Clean Reads percentage, Clean Data percentage) were obtained by filtering the 3′ end connector sequence and Reads with an average quality score lower than Q20 using Fastp software, and then calculating the clean data volume / original data volume × 100%.

[0120] (3) The 70-day and 120-day sample data and the 56-day sample data were sequenced and analyzed using the same procedure, and the data range was consistent with the quality control level of the 56-day sample.

[0121] Characterization Example 1: Clean data from 30 Muscovy duck breast muscle samples after S1 quality control, i.e., the sequencing sequences obtained after filtering 3′ adapter sequences and reads with an average quality score lower than Q20 using Fastp software;

[0122] FastQC v0.11.9 was used, following Illumina 1.9 encoding rules. Base quality value (Q value) and sequencing error rate were analyzed. The conversion formula is:

[0123] in, This represents the sequencing error rate of the base (i.e., the probability of incorrect base identification). The value range is 0-40, corresponding to an error rate range of 100-0.0001%;

[0124] The FASTQ format clean data of 30 samples were imported into the analysis software. The software automatically identified the base sequence and corresponding quality label of each Read. According to the 5′→3′ direction of the bases in the Read, the bases of each Read were grouped equally at positions (position 1 to position 151), forming a total of 151 independent position groups. For each position group (e.g., position 1, position 2... position 151), the Q value of the bases at that position in all sample clean data was calculated, and the arithmetic mean was calculated (the average Q value was calculated only once for each position group, without repeated sampling or weighting).

[0125] like Figure 1 The base quality distribution map shown shows that the average Q value of the full read length (positions 1-151) is distributed in the range of 38-40, with no positions where the Q value is <30. The corresponding sequencing error rate is ≤0.000158%, which is much higher than the conventional standards of Q20≥98% and Q30≥95% for transcriptome sequencing.

[0126] The 5′ end (positions 1-6) showed no low-quality initiation regions, with an average Q value ≥ 38. The average Q value in the middle region (positions 7-146) remained stable at around 40, without any sudden high error rate fluctuations. This demonstrates that the sequencing instrument was highly stable, the reagents were of uniform quality, and the samples were uncontaminated. The average Q value at the 3′ end (positions 147-151) decreased slightly from 40 to around 38, which is consistent with the inherent characteristic of a slight increase in the terminal error rate due to chemical reagent consumption in high-throughput sequencing, and there was no abnormal decrease. These characterization results demonstrate that the sequencing data of the 30 Muscovy duck breast muscle samples were of excellent quality, without base quality bias, and could effectively eliminate the interference of low-quality sequences on subsequent analysis.

[0127] Characterization Example 2: Clean data from 30 Muscovy duck breast muscle samples, identical to Characterization Example 1, passed S1 quality control. The analysis software was RSeQC v4.0.0. Data were sorted from low to high according to gene FPKM-normalized expression levels and divided into 4 groups:

[0128] Q1 expression level was at least 25%;

[0129] Q2 expression level 25-50%;

[0130] Q3 expression level 50-75%;

[0131] Q4 expression level was the highest at 25%;

[0132] The saturation characteristics of genes at different expression levels were analyzed separately. The gene expression level was calculated 20 times for each sampling ratio, and the average value was used for error analysis.

[0133] Gene FPKM-normalized expression data were extracted from 30 samples. Low-expression genes with FPKM=0 were filtered out (to avoid interference from non-expressed genes in error calculation), and the effective expression gene set (a total of 14,691 genes) was retained.

[0134] Using RSeQC software, clean data were randomly resampled at gradient ratios of 5%, 10%, 15%...100% (each ratio was independently repeated 20 times). The expression level corresponding to the 100% sampling ratio was considered the actual expression level (true value). For each resampling ratio (e.g., 5%, 10%...95%), the average expression level of the four gene groups Q1-Q4 was calculated and compared with the actual expression level at 100% sampling. The relative error was calculated using the following formula:

[0135]

[0136] The relative error for each scale is the arithmetic mean of 20 repeated calculations.

[0137] like Figure 2 The saturation analysis plot shows that the relative errors of the four gene groups (Q1-Q4) gradually decrease with the increase of the resampling ratio. When the sampling ratio is ≥50%, the error decreasing trend slows down significantly. When the sampling ratio is ≥80%, the error is basically stable (fluctuation range ≤0.5%). Among all resampling ratios, the relative errors of the four gene groups Q1-Q4 are all ≤5%. Among them, Q4 (high expression gene) has the best stability when the sampling ratio is ≥30%, with a relative error of ≤2%. Q1 (low expression gene) has a relative error of ≤4% when the sampling ratio is ≥60%. Even the gene with the lowest expression level can accurately capture expression features with the existing data. When sampling at 100%, the relative errors of the four gene groups are all ≤1%, further verifying that the data is fully saturated.

[0138] The characterization results demonstrate that the sequencing data from 30 Muscovy duck breast muscle samples has reached the saturation threshold, and the existing data is sufficient to support the accurate quantification of all genes (especially low-expression genes) without the need to increase the sequencing depth.

[0139] Characterization Example 3: Thirty Muscovy duck breast muscle samples, identical to those in Characterization Examples 1 and 2, were analyzed using BAM format files after S2 genome alignment. Additionally, a gene structure annotation file (GTF format, containing gene CDS, intron, UTR, and other location information) of the Muscovy duck reference genome was prepared. The analysis software was RSeQC v4.0.0. The genome was divided into 100 equal parts according to the base length along the 5′→3′ direction (i.e., each part represents 1% of the total gene length), with the x-axis labeled from 0 to 100% (0% for the 5′ end and 100% for the 3′ end).

[0140] The BAM files of 30 samples were quality filtered to retain only unique alignment sequences (multiple alignment sequences were removed to avoid duplicate counting). The BAM files were sorted (by chromosome position) and indexed using SAMtools software. Protein-coding genes were extracted from the GTF annotation files (ncRNA, pseudogenes, etc. were excluded), and genes with expression levels FPKM ≥ 0.1 were retained (to avoid coverage noise interference from low-expression genes). Finally, 12,436 effective genes were obtained.

[0141] Each effective gene was divided into 100 segments along the 5′→3′ direction (segment length = total gene length × 1%). For each sample, the number of reads that aligned to each gene among the 100 segments was counted (only effective reads that cross segments were counted, and reads that cross segments were not counted repeatedly). The coverage of the 100 segments of each gene was normalized (i.e., the coverage of each segment ÷ the total coverage of the 100 segments of the gene) to eliminate the numerical bias of coverage caused by the difference in expression levels of different genes, so that the coverage range of all genes was uniformly 0-1.

[0142] The normalized coverage of 30 samples was summarized, and the average coverage was calculated according to the gene equivalence positions (1-100). That is, for the first part (1% of the 5′ end region of all genes), the arithmetic mean of the normalized coverage of 30 samples was calculated; the second part to the 100th part were calculated in the same way, and finally the average coverage data of 100 positions were obtained.

[0143] like Figure 3 The gene coverage uniformity plot shown has a generally stable distribution of peaks, without obvious single-peak bias or local spikes or drops. The normalized average coverage fluctuates between 0.75 and 0.85, and there is no biased feature where the 5′ or 3′ end coverage is significantly higher than that of the middle region.

[0144] The average coverage of the 5′ end of the gene (0-10% region) was 0.78-0.82, with no sudden drop in coverage at the start end (excluding 5′ sequence loss due to RNA degradation); the average coverage of the middle region of the gene (10-90% region) was stable at 0.80-0.85, with a fluctuation range of ≤0.03, proving that reads uniformly cover the core coding region; the average coverage of the 3′ end of the gene (90-100% region) was 0.76-0.80, with a slight decrease but no significant drop;

[0145] The characterization results demonstrate that the genome alignment of S2 is effective. The uniform coverage of Reads on genes ensures that the Read Count value calculated by HTSeq software can accurately reflect the gene expression level, avoiding the problem of overestimation of short gene expression due to 5′ end bias coverage and underestimation of long gene expression due to 3′ end bias coverage.

[0146] Characterization Example 4: Gene expression data of 30 Muscovy duck breast muscle samples, the same as those in Characterization Examples 1-3, after S2 FPKM normalization. At the same time, low-expression genes with FPKM=0 need to be removed (to avoid noise interference from genes without expression). Finally, 12436 effective expression genes are retained.

[0147] The expression matrix was standardized (to eliminate the magnitude difference in expression levels of different genes and ensure the objectivity of the analysis), and the first two principal components (PC1 and PC2) were retained, requiring that their cumulative contribution rate be ≥70% (which can explain most of the data differences).

[0148] To avoid interference from genes specifically expressed in a single sample, at least three genes with FPKM ≥ 0.1 were retained from all samples. An expression matrix of 12436 × 30 was constructed with genes as rows and samples as columns. PCA was performed on the standardized expression matrix, and the contribution rate of each principal component and the score of each sample in each principal component were extracted and plotted (PC1 and PC2 values, see Table 3).

[0149] Table 3:

[0150] As shown in Table 3, all samples in group G had negative PC1 values ​​(ranging from -28.4657 to -20.8433) and negative PC2 values ​​(ranging from -11.6088 to -7.5412), corresponding to... Figure 4 Left-side clustering; in group M, all samples have negative PC1 values ​​but smaller absolute values ​​(-23.0756 to -11.7853), and all have positive PC2 values ​​(9.8725-13.6842), corresponding to... Figure 4 The right-side clustering shows no overlap or mixing between the two groups.

[0151] like Figure 4 The PCA analysis diagram and Table 3 above demonstrate that the sample pretreatment and grouping design of S2 are effective. The differences in expression levels between samples are mainly determined by the experimental grouping, rather than individual differences, experimental errors, or sample contamination.

[0152] Characterization Example 5: Core results of S3 differential expression analysis of 30 Muscovy duck breast muscle samples identical to those in Characterization Examples 1-4, including Gene_ID, log2FoldChange (fold change in expression), and P-value (statistical significance) for all genes. Differential gene screening thresholds |log2FoldChange|>1 (expression difference threshold) and P-value<0.05 (statistical significance threshold) are used. Genes that meet both thresholds are considered significantly differentially expressed genes (further divided into upregulated genes log2FoldChange>1; downregulated genes log2FoldChange<-1). Genes that do not meet these thresholds are considered not significantly differentially expressed genes.

[0153] The core columns Gene_ID, log2FoldChange, and P-value were extracted from the differential analysis results table. Genes with empty P-values ​​or meaningless log2FoldChange were removed, ultimately retaining 14,989 valid genes. Genes were then categorized according to screening thresholds.

[0154] Upregulated genes log2FoldChange>1 and P-value<0.05 (151-208 genes, slightly different depending on the comparison group);

[0155] Downregulated genes log2FoldChange < -1 and P-value < 0.05 (147-246 genes);

[0156] Genes with no significant differences did not meet the above dual thresholds (14,691).

[0157] Table 4 shows the key information on differentially regulated genes (Top 10 genes upregulated and downregulated):

[0158] like Figure 5 The volcano plot shown has scatter points clustered in a volcano-like pattern. Significantly upregulated genes (red) cluster on the right side of the horizontal axis (log2FoldChange>1), significantly downregulated genes (blue) cluster on the left side of the horizontal axis (log2FoldChange<-1), and genes with no significant difference (gray) cluster in the middle region of the horizontal axis (-1≤log2FoldChange≤1). The boundaries are clear. Taking the G_56_vs_M_56 comparison group as an example, there are 151 significantly upregulated genes, 147 significantly downregulated genes, and 14691 genes with no significant difference. The dual-threshold screening effect is clear. The top 10 upregulated and downregulated genes (sorted by P-value) are all marked next to the corresponding scatter points and are all distributed in the top region of the volcano plot (higher vertical axis values, i.e., smaller P-values, stronger significance). This characterization result ( Figure 5 (Table 4) demonstrates that the differential gene screening of S3 is effective. The dual threshold can accurately capture genes with significant differences in expression between groups and exclude interference from genes with no biological significance.

[0159] Characterization Example 6: The set of significantly differentially expressed genes in Muscovy ducks selected by S3 (|log2FoldChange|>1 and P-value<0.05, totaling 298-400 genes, with slight differences due to different comparison groups), and the FPKM-normalized expression data of these genes in 30 Muscovy duck breast muscle samples (the influence of gene length and sequencing depth has been eliminated). The expression levels of differentially expressed genes were Z-score normalized (to eliminate the difference in expression magnitude between samples and facilitate cross-sample comparison of expression patterns).

[0160] Euclidean distance was used to calculate the similarity of gene expression patterns, and gene clusters were divided by hierarchical clustering (Complete Linkage, longest distance method). The rationality of clustering was verified by grouping samples (G group / M group) using the same distance and clustering method. A red-blue two-color mapping was used, where red represents high gene expression, blue represents low gene expression, and gray represents no expression or extremely low expression.

[0161] Significantly differentially expressed genes were screened from the S3 results, and a gene × sample expression level matrix was constructed (rows = differentially expressed genes, columns = 30 samples, such as G1_56, M1_56, etc.). The expression level of each gene in the matrix was standardized using the following formula: ,

[0162] in, This represents the original expression level of gene i in sample j. The average expression level of gene i across all samples. Let be the standard deviation of gene i's expression level across all samples. After standardization, the expression level distribution of each gene has a mean of 0 and a standard deviation of 1, which facilitates cross-sample comparison of expression patterns.

[0163] Based on the standardized expression matrix, the Euclidean distance between all genes is calculated, and a gene clustering tree is generated through hierarchical clustering (CompleteLinkage) to group genes with similar expression patterns into the same cluster;

[0164] Based on the same matrix, calculate the Euclidean distance between all samples, generate a sample clustering tree using the same clustering method, and verify whether samples in the same group are clustered and whether samples in different groups are separated.

[0165] The normalized expression levels of differentially expressed genes in representative samples are shown in Table 5. Table 5:

[0166] It should be noted that the table above is a core data example, showing the normalized expression levels of differentially expressed genes in representative samples. The complete data can be reproduced from the differentially expressed gene results in S3 using the method described in this characterization example.

[0167] like Figure 6 The differential gene clustering heatmap shown shows that the sample clustering tree on the left side of the heatmap indicates that all samples (15) in group G are tightly clustered into one large branch, and all samples (15) in group M are tightly clustered into another large branch. There is no overlap or mixing between the two groups, which proves that the expression patterns of samples in the same group are highly similar, and the expression patterns of different groups are significantly separated.

[0168] The gene clustering tree on the right divides 298-400 differentially expressed genes into 6-8 core gene clusters. The expression patterns of genes within each cluster are highly consistent (e.g., cluster 1 is highly expressed in group G and lowly expressed in group M, while cluster 2 is the opposite), reflecting the synergistic expression characteristics of functionally related genes.

[0169] The heatmap color distribution shows obvious group specificity. The gene clusters in the G group sample area are mainly red (high expression), while the corresponding area in the M group sample area turns blue (low expression), and vice versa. This intuitively shows that the expression pattern of differentially expressed genes is highly correlated with the group.

[0170] The characterization results ( Figure 6 (Table 5) demonstrates that the differentially expressed genes screened in S3 have clear biological significance, and their expression patterns are highly correlated with the experimental groups rather than fluctuating randomly. At the same time, the gene clustering characteristics provide a set of co-expressed genes for the subsequent functional enrichment analysis in S4 (such as the KEGG pathway and PPI interaction), ensuring the specificity and reliability of functional validation.

[0171] Characterization Example 7: Through KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analysis, the core metabolic pathways and signal transduction pathways involved by significantly differentially expressed genes in Muscovy ducks were identified, verifying that the differentially expressed genes screened by S3 have clear biological functions, rather than random fluctuations;

[0172] The significantly differentially expressed genes in Muscovy ducks selected in step S3 (|log2FoldChange|>1 and P-value<0.05, a total of 298-400 genes) are obtained from the KEGG annotation information (KO number, pathway attribution) of the Muscovy duck reference genome or mapped from the KEGG database of closely related species (such as chickens and turkeys) through homology comparison (BLAST).

[0173] Using P-value < 0.05 as the criterion for significant enrichment, the hypergeometric distribution test was used to evaluate the significance of differentially enriched genes in each KEGG pathway. A bar chart was used to display the top enriched pathways and show the specific enrichment location of genes in the pathway.

[0174] Gene_ID of differentially expressed genes was extracted and compared with the KEGG database (if Muscovy ducks do not have specific annotations, the KO numbers of closely related species with sequence similarity ≥80% were selected). The KO (KEGG Ortholog) annotation and the pathway information corresponding to each differentially expressed gene were obtained. All annotated genes obtained from sequencing (a total of 14,989) were used as the statistical benchmark for the hypergeometric distribution test.

[0175] The hypergeometric distribution test is used to calculate the p-value for each KEGG Pathway, using the following formula:

[0176] Where N is the total number of background genes, M is the number of background genes annotated to this pathway, n is the total number of differentially expressed genes, and k is the number of differentially expressed genes annotated to this pathway;

[0177] Sorting by P-value from smallest to largest, we screened for significantly enriched pathways with P-value < 0.05, prioritizing pathways related to Muscovy duck breast muscle development and lipid metabolism (such as the PPAR signaling pathway and cholesterol metabolism pathway).

[0178] The top 5 significantly enriched KEGG pathway data are shown in Table 6.

[0179] Table 6:

[0180] It should be noted that the table above is a core data example, showing the key statistics of 5 core pathways selected from the Top 10 significantly enriched pathways. The complete data can be reproduced from the differential gene annotation results using the method described in this characterization example.

[0181] like Figure 7 The KEGG enrichment analysis diagram shown indicates that the top 10 enriched pathways are all highly related to Muscovy duck breast muscle development and lipid metabolism (such as the PPAR signaling pathway and cholesterol metabolism pathway), with P-values ​​all < 1 × 10⁻⁶. -3 The enrichment was significant; the enrichment factor was >0.4, indicating that the enrichment of differentially expressed genes in these pathways was significantly higher than the background level.

[0182] The characterization results ( Figure 7 (Table 6) demonstrates that the differentially expressed genes screened in S3 have clear biological functions and are significantly enriched in core pathways related to Muscovy duck pectoral muscle development and lipid metabolism, rather than exhibiting random fluctuations. At the same time, the enrichment results provide a functional screening basis for the subsequent identification of core genes in S5, prioritizing differentially expressed genes enriched in key pathways as core candidates, thus ensuring the biological effectiveness and application value of the entire differential gene screening method.

[0183] Characterization Example 8: Through protein-protein interaction (PPI) network analysis, the interaction relationships between proteins encoded by significantly differentially expressed genes in Muscovy ducks were explored, verifying that these differentially expressed genes do not exist in isolation, but form a functionally coordinated regulatory network.

[0184] The significantly differentially expressed genes in Muscovy ducks selected by S3 (|log2FoldChange|>1 and P-value<0.05, totaling 298-400 genes) were obtained using general PPI databases such as STRING v11.5 and BioGRID. If there was no specific interaction data for Muscovy ducks, it was obtained by mapping from the interaction networks of closely related species (such as chickens and turkeys) through homology alignment (BLAST, sequence similarity ≥80%). Only high-confidence interaction relationships with a confidence score ≥0.7 were retained (excluding random interaction interference). Protein nodes with a degree value (number of connections) ≥10 were regarded as regulatory hubs. Nodes represent proteins encoded by differentially expressed genes, edges represent interaction relationships, node size is positively correlated with degree value, and color represents functional classification.

[0185] Protein sequences of differentially expressed genes were extracted and compared with the PPI database to obtain interaction pairs (protein A - protein B) and their corresponding confidence levels. Interactions with confidence levels <0.7 were removed, and the high-confidence interaction network was retained. A PPI network was constructed using proteins encoded by differentially expressed genes as nodes and high-confidence interaction relationships as edges. A module detection algorithm (such as MCODE) was used to discover tightly connected functional modules (each module has dense protein interactions and participates in the same biological process) from the network. The number of connections (degree value) of each node was counted, and core nodes with a degree value ≥10 were selected as regulatory hubs.

[0186] The node data is shown in Table 7;

[0187] Table 7:

[0188] It should be noted that the table above is a core data example, showing the key statistics of the Top 5 core nodes. The complete data can be reproduced from the differential gene interaction results using the method described in this representation example.

[0189] like Figure 8 The PPI interaction network diagram shown is composed of 3-4 core functional modules (such as lipid metabolism module, myofibril development module, and signal regulation module). There are dense protein interactions within the modules and a small number of connections between the modules, reflecting functional synergy. The top 10 core nodes with the highest degree values ​​(such as PPARγ and MYOD1) are all regulatory genes that are highly related to Muscovy duck pectoral muscle development and lipid metabolism. They connect a large number of differentially expressed genes and are the regulatory hubs of the network. The confidence of all edges is ≥0.7 and the confidence of interactions between core nodes is ≥0.9, ensuring the reliability of the interaction relationships.

[0190] The characterization results ( Figure 8 (Table 7) demonstrates that the differentially expressed genes screened in S3 have functional synergy and are not isolated genes with random fluctuations, but rather form an interaction network with core regulatory genes as the hub, synergistically participating in key biological processes such as pectoral muscle development and lipid metabolism. At the same time, the core nodes provide a priority screening basis for the subsequent identification of core differentially expressed genes in S5, ensuring the biological effectiveness and application value of the entire differential gene screening method.

[0191] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions will not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for screening differentially expressed genes in Muscovy ducks based on transcriptome sequencing, characterized in that, Eukaryotic transcriptome sequencing technology was used to analyze Muscovy duck breast muscle samples. Differential gene screening was achieved through data quality control, genome alignment and expression level standardization, initial screening of differentially expressed genes, multi-dimensional validation, and identification of core differentially expressed genes. Specifically, this included: The original data shows that the percentage of fuzzy bases (N%) is 0.0081-0.0366%, the GC content is 48.10-49.39%, the percentage of Q30 bases is ≥95.0%, and the percentage of Q20 bases is ≥98.0%. Clean data refers to Clean Reads accounting for ≥98.0% and Clean Data accounting for ≥97.8%. Clean Reads refers to the number of high-quality sequencing sequences after filtering, and Clean Data refers to the total number of bases in the high-quality sequencing sequences after filtering. The differential gene screening threshold is |log2FoldChange|>1 and P-value<0.

05.

2. The method for screening differentially expressed genes in Muscovy ducks based on transcriptome sequencing according to claim 1, characterized in that, The method includes the following steps: S1: Quality control processing of raw transcriptome sequencing data from Muscovy duck breast muscle samples; The raw sequencing data in FASTQ format was filtered using Fastp software or equivalent transcriptome data filtering software to remove sequences with adapters at the 3′ end and reads with an average quality score lower than Q20, thus obtaining clean data. The clean data were sequentially subjected to base quality distribution detection, base content distribution detection, average quality distribution of reads detection, gene coverage uniformity analysis and saturation analysis. The clean data were confirmed to be unbiased and the sequencing data volume reached the saturation threshold by the gene coverage uniformity curve showing no single peak bias and the relative error of gene expression level in the saturation analysis being ≤5%. S2: Genome alignment and gene expression level normalization analysis of clean data; Using HISAT2 software or a functionally equivalent genome alignment software, clean data that passed quality control were aligned to the Muscovy duck reference genome. Default parameters were used for non-strand-specific libraries, while the --rna-strandness RF parameter was used for the first-strand library and the --rna-strandness FR parameter was used for the second-strand library. The overall alignment rate between the clean data and the reference genome was 59.11%–67.69%. The original gene expression level was counted using the Union statistical scheme of HTSeq software or a functionally equivalent expression level statistical method. When constructing a chain-specific library, it is necessary to determine whether the Read and Feature directions in the annotation are consistent. The original gene expression level was converted into the number of fragments per thousand bases of transcript per million aligned fragments using the FPKM method to eliminate the influence of gene length and sequencing depth. Pearson correlation coefficient was used to analyze the correlation of gene expression levels among samples, and samples with correlation coefficients in the range of 0.8-1.0 were screened. Principal component analysis (PCA) was used to verify the effectiveness of sample grouping and clustering, ensuring that samples in the same group clustered together and samples in different groups were clearly separated. S3: Initial screening of differentially expressed genes in Muscovy ducks; Deseq2 software or functionally equivalent differential expression analysis software was used to perform inter-group differential expression analysis on the original expression levels of genes that passed S2 validation. With |log2FoldChange|>1 and P-value<0.05 as the screening thresholds, differentially expressed genes in Muscovy ducks were obtained. The total number of differentially expressed genes in different comparison groups was 298-400, of which 151-208 were upregulated and 147-246 were downregulated. Collect common and unique genes among different comparison groups; S4: Multidimensional verification of differential genes in Muscovy ducks; Protein-protein network interaction analysis was performed on differentially expressed genes using the STRING database, and protein-protein interaction pairs with a score > 0.95 were selected. If the STRING database did not have Muscovy duck PPI information, chicken, turkey, or Muscovy duck species with close kinship were selected. The BLAST algorithm was used to compare the Muscovy duck protein sequence with the protein sequence of the selected species, and homologous protein interaction relationships with similarity ≥ 80% were selected and mapped back to Muscovy duck as Muscovy duck protein-protein interaction pairs with a score > 0.

95. ClusterProfiler software or equivalent enrichment analysis software were used to perform GO enrichment analysis and KEGG enrichment analysis. Significantly enriched GO terms and KEGG pathways were screened with a P-value < 0.

05. The GO enrichment analysis analyzed the functions of differentially expressed genes from three dimensions: molecular function, cellular components, and biological processes. The KEGG enrichment analysis completed KO annotation and metabolic pathway annotation. GSEA analysis was conducted, using |NES|>1, NOM p-val<0.05, and FDR q-val<0.25 as thresholds to screen genes with insignificant differential expression but important biological significance. S5: Identification of core differentially expressed genes in Muscovy ducks; To further characterize the regulatory features of differentially expressed genes, at least one of the following analyses will be performed: By comparing with the AnimalTFDB database, we can identify transcription factors and their families in differentially expressed genes. The rMATS software was used to identify five types of alternative splicing events: exon skipping, variable 5′ splicing site, variable 3′ splicing site, mutually exclusive exons, and retained introns. Differential alternative splicing events were screened based on IncLevelDiff≥0.1 and FDR<0.

05. SNP and InDel detection of differentially expressed genes were performed using Varscan. The filtering criteria were Q > 20, number of reads for covered sites > 8, number of reads for supporting mutation sites > 2, and P-value < 0.

01. Functional annotation of variant sites was performed using ANNOVAR software. The annotation scope included intron regions, intergenic regions, exon regions, 5′UTR regions, 3′UTR regions, splice regions, upstream regions, downstream regions, and ncRNA exon regions. Based on the differential screening results of S3, the functional enrichment and protein interaction results of S4, and the above analysis information, genes that simultaneously meet the criteria of significant differential expression, high-confidence protein interaction relationship, and enrichment in functional pathways related to Muscovy duck pectoral muscle traits are prioritized as core differential genes of Muscovy duck. Genes that are significantly differentially expressed and enriched in key functional pathways but for which transcription factor regulation or gene variation has not been detected can also be included in the core differentially expressed gene set, taking into account their biological function and their core position in the pathway.

3. The method for screening differentially expressed genes in Muscovy ducks based on transcriptome sequencing according to claim 2, characterized in that, The HISAT2 software described in S2 simultaneously counts the proportion of multiple alignment sequences and the proportion of unique alignment sequences, and the distribution of alignment regions in the CDS region, introns, and intergenic regions conforms to the characteristics of eukaryotic transcriptome sequencing.

4. The method for screening differentially expressed genes in Muscovy ducks based on transcriptome sequencing according to claim 2, characterized in that, The rMATS software described in S5 uses the Junction Counts quantitative method, through a formula. The expression level of alternative splicing is calculated, where I is the number of transcript reads containing the alternative splicing region, LI is the effective length of the corresponding transcript, S is the number of transcript reads that skip the alternative splicing region, and LS is its effective length.

5. The method for screening differentially expressed genes in Muscovy ducks based on transcriptome sequencing according to claim 2, characterized in that, The core pathways for KEGG enrichment described in S4 include the PPAR signaling pathway, the adrenergic signaling pathway, the estrogen signaling pathway, the cholesterol metabolism pathway, and the fatty acid degradation pathway.

6. The method for screening differentially expressed genes in Muscovy ducks based on transcriptome sequencing according to claim 2, characterized in that, The saturation analysis described in S1 uses RSeQC software to resample clean data at ratios of 5%, 10%, 15%...100%. The expression levels of all genes are calculated 20 times at each sampling ratio. The relative error is calculated by comparing the gene expression levels at different sampling ratios with the actual expression levels at 100% sampling. All genes are divided into four groups, Q1-Q4, according to their expression levels. When the relative error of all four gene groups is ≤5%, the data volume is considered to have reached the saturation threshold.

7. The method for screening differentially expressed genes in Muscovy ducks based on transcriptome sequencing according to claim 2, characterized in that, The PCA principal component analysis described in S2 retains the first two principal components, whose cumulative contribution rate is ≥70%.

8. The application of the Muscovy duck differential gene screening method according to any one of claims 1-7 in Muscovy duck molecular breeding, genetic improvement of Muscovy duck breast muscle traits, and mining of Muscovy duck functional genes.

Images