Computational protocol: Genome wide DNA methylation profiles reveal novel candidate genes associated with meat quality at different age stages in hens

Similar protocols

Protocol publication

[…] The raw reads produced by the Illumina HiSeq in FastQ format were first preprocessed through in-house Perl scripts. We filtered out reads that contained adapters, N (unknown bases) and those in which over 50% of the sequence exhibited a low quality value (PHRED score ≤5). At the same time, the Q20, Q30, and GC content of the data were calculated. The remaining reads that passed the filters were called clean reads, and all of the subsequent analyses were based on these.The methylation data were aligned to the reference genome Gallus gallus-4.0 by Bismark software (version 0.12.5). The results of the sequencing and the reference genome were transformational (C to T and G to A). The transformation of the sequencing results and the genome were a pairwise comparison, which was then indexed using bowtie2. The best one of the four parallel alignment results was selected as the final alignment result. The proportion of the number of aligned reads in the total number of reads was regarded as the mapping rate. The same reads that aligned to the same regions of the genome were regarded as duplicated reads. The sequencing depth and coverage were summarized using the deduplicated reads.In the second-generation sequencing, the sequencing depth of each locus was not the same, and the proportion of methylation to unmethylation was different at different sites. However, the probability of methylation at each site obeyed the Binomial distribution. The Binomial Distribution, or the Bernoulli Experiment, which repeats n times, is the probability distribution of a discrete random variable with a wide range of applications. Based on the methylation data of the Bismarck software, the methylated cytosine and unmethylated cytosine frequencies at each site were tested for Binomial distributions in order to identify whether the site was a true methylation site. To find the exact methylation site, a set of thresholds were set during the analysis process: (1) Sequencing with a depth greater than or equal to 5 and (2) a q-value less than or equal to 0.05.For the identified methylation sites, the methylation level (ML) was calculated as follows (mC and umC represent the number of methylated C and unmethylated C, respectively): According to previous studies, due to the influence of the bisulfite conversion rate, the level of methylation needed to be corrected. ML_corrected represents the corrected methylation level, r represents the bisulfite non-conversion rate, and the corrected ML was estimated as: [...] To determine the differential methylation regions between the two age groups, the present study used swDMR software (, and a sliding-window approach (Select every 1000 bp as a window, with 100 bp as the step length) was performed to identify the differential methylation regions in the genome. Those regions, with corrected P-values less than 0.05, a read coverage greater than 5, fold-change values greater than 2 and a false discovery rate (FDR) less than 0.05, were regarded as differential methylation regions (DMRs). When the region where a DMR and a specific gene function element overlapped, the corresponding gene was selected as the DMR-related gene, namely, a differentially methylated gene (DMG). [...] For the WGBS validation, 250 ng of the pooled genomic DNA (diluted to 250 ng/μL) from each group was treated with bisulfite sodium using the MethylEdge® Bisulfite Conversion System (Promega Research, Wisconsin, USA) according to the manufacturer’s instructions. The bisulfite-treated DNA was used for PCR. Two hypermethylated regions and one hypomethylated region were selected randomly. Online MethPrimer software ( was used to design the bisulfite sequencing PCR (BSP) primers (). The BSP reaction was performed in 50 μL containing 100 ng genomic DNA, 1 μL of each primer (20 μM), 0.25 μL of TaKaRa EpiTaq HS (5 U/μL), 5 μL of 10× EpiTaq PCR Buffer (TaKaRa, Dalian, China), 5 μL of MgCl2 (25 mM), and 6 μL of a dNTP mixture. The PCR was performed using the following program: 94 °C for 5 min; 40 cycles of 94 °C for 1 min; 60 °C for 1 min; 72 °C for 30 s; and ending with incubation at 72 °C for 10 min. The PCR products were detected by 2% agarose gel electrophoresis and cloned into the pMD18-T vector (TaKaRa). Ten clones were selected for each gene and were subsequently sequenced by BGI (Wuhan, Hubei, China). All of the sequences were analysed by the online software QUMA ( All the BSP experiments were repeated at least three times. […]

Pipeline specifications

Software tools Bismark, Bowtie2, swDMR
Application BS-seq analysis
Organisms Gallus gallus