Percentile Time Series Analysis | Feature Methylation Rate Analysis | Whole-genome Methylated Fraction Analysis | Whole-genome Methylation Rate Analysis | Phasing Analysis | Installation | Notice | Citation | Reference
This software analyzes the methylated fraction of GATC sites in BigWig files generated by the Snakemake pipeline snakemakeMethylFrac. The default resource files are provided to analyze the data in human MCF7 and MCF10A cells, but the users could also use their only annotation files to analyze in other organisms. The resource files for annotations are in the resources directory, and the metadata of the resource files are stored in the config/resource.toml file. The users could modify the config/resource.toml file to use their own annotation files. The config/example_sample.toml file is the metadata of the specific samples. Currently, our softwares used t2t annotation for most of the analysis (centromeric elements, ChromHMM, whole-genome maps), and the hg38 annotation is only used for the analysis of general genetic features (e.g., gene body, promoter, enhancer, etc.) and the phasing plots relative to genes and CTCF sites. The users could modify the config/example_sample.toml file to analyze their own samples. The software has the following functions:
Usage | Description | Options
analzyeMethylFrac featPercentile -c {path/to/sample.toml} -r [full/calc/plot]
This subcommand calculates the percentiles (5%, 15%, 25%, 50%, 75%, 85%, 95%) in the provided features in config/resource.toml file. The provided resource files will process the following features: General Features in Hg38 (Gencode v43 annotation) including cpgislands, enhancers, silencers, tRNA genes, origin of replications, protein coding genes, centromeres and whole-genome data (details in the Notice section); Detailed centromeirc elements in T2T v1.1 assembly; Histone Marks and ChromHMM states. The results are stored in CSV files within the {figure_dir}/source directory ({figure_dir} is specificed in the config toml file).
This subcommand also plots the the percentiles in Time Series according to the sample sheet (sample/example_sample_sheet.tsv), which is the same sample sheet in the snakemake workflow. The plots are stored in separate folders in the {figure_dir}. Furthermore, median plots comparing different features are also generated in the {figure_dir}/median.
-c --configfile: The path to the sample toml file. The example sample toml file is config/example_sample.toml.
-r --run [full/calc/plot]: The run mode of the subcommand. The full mode will calculate the percentiles and plot the time series. The calc mode will only calculate the percentiles. The plot mode will only plot the time series. The plot must be used after the calc mode.
--help: Print the help message.
Optional arguments:
-R --resourcefile: The path to the resource toml file. The resource toml file for the provided resource files is config/resource.toml.
Usage | Description | Options
analzyeMethylFrac featRate -c {path/to/sample.toml} -r [full/calc/plot]
This subcommand calculates the methylation rate in the provided features in config/resource.toml file, which are listed in Percentile Time Series Analysis section. This subcommand calculated the methylation rate in each feature using the median methylated fractions in each feature, which are calculated in the featPercentile subcommand. The results are stored in CSV files within the {figure_dir}/source directory ({figure_dir} is specificed in the config toml file). The relative methylation rate is calculated by dividing the methylation rate using the whole-genome methylation rate (rate calculated by the median methylated fraction).
This subcommand also plots the methylation rate in each feature. The plots are stored in {figure_dir}/rate.
-c --configfile: The path to the sample toml file. The example sample toml file is config/example_sample.toml.
-r --run [full/calc/plot]: The run mode of the subcommand. The full mode will calculate the methylation rate. The calc mode will only calculate the methylation rate. The plot mode will only plot the methylation rate. The plot must be used after the calc mode.
--help: Print the help message.
Usage | Description | Options
analzyeMethylFrac genomeFrac -c {path/to/sample.toml} -r [full/calc/plot]
This subcommand calculates the average methylated fraction in the whole genome in 100 kb bins. The results are stored in CSV files within the {figure_dir}/source directory ({figure_dir} is specificed in the config toml file).
This subcommand also plots the whole genome methylation fraction map of each sample. The plots are stored in {figure_dir}/frac_map. The centoremeres are highlighted in red blocks in the plots.
-c --configfile: The path to the sample toml file. The example sample toml file is config/example_sample.toml.
-r --run [full/calc/plot]: The run mode of the subcommand. The full mode will calculate the methylation rate. The calc mode will only calculate the methylation rate. The plot mode will only plot the methylation rate. The plot must be used after the calc mode.
--help: Print the help message.
Optional arguments:
-R --resourcefile: The path to the resource toml file. The resource toml file for the provided resource files is config/resource.toml.
Usage | Description | Options
analzyeMethylFrac genomeRate -c {path/to/sample.toml} -r [full/calc/plot]
This subcommand calculates the methylation rate in the whole genome in 100 kb bins. The results are stored in CSV files within the {figure_dir}/source directory ({figure_dir} is specificed in the config toml file). The methylation rate is calculated using the median methylated fractions in each 100 kb bin, which are calculated in the genomeFrac subcommand. The relative rate is calculated by dividing the methylation rate using the whole-genome methylation rate in Feature Methylation Rate Analysis section.
This subcommand also plots the whole genome methylation rate map of each sample. The plots are stored in {figure_dir}/rate_map. The centoremeres are highlighted in red blocks in the plots.
-c --configfile: The path to the sample toml file. The example sample toml file is config/example_sample.toml.
-r --run [full/calc/plot]: The run mode of the subcommand. The full mode will calculate the methylation rate. The calc mode will only calculate the methylation rate. The plot mode will only plot the methylation rate. The plot must be used after the calc mode.
--help: Print the help message.
Optional arguments:
-R --resourcefile: The path to the resource toml file. The resource toml file for the provided resource files is config/resource.toml.
Usage | Description | Options
analzyeMethylFrac phasing -c {path/to/sample.toml} -r [full/calc/plot]
This subcommand calculates the phasing of Methylated Fraction relative to the TSS of protein-coding genes as well as that relative to the CTCF sites. The results are stored in CSV files within the {figure_dir}/source directory ({figure_dir} is specificed in the config toml file). The protein-coding genes are grouped by active and inactive NDRs based on the ATAC-seq data (details in Notice). The average methylated fraction in +/- 1 kb of the TSS and CTCF sites are calculated. The average methylated fractions are further smoothed in 21-bp windows using moving average.
This subcommand also plots the phasing of Methylated Fraction relative to the TSS of protein-coding genes as well as that relative to the CTCF sites. The plots are stored in {figure_dir}/phasing.
-c --configfile: The path to the sample toml file. The example sample toml file is config/example_sample.toml.
-r --run [full/calc/plot]: The run mode of the subcommand. The full mode will calculate the methylation rate. The calc mode will only calculate the methylation rate. The plot mode will only plot the methylation rate. The plot must be used after the calc mode.
--help: Print the help message.
Optional arguments:
-R --resourcefile: The path to the resource toml file. The resource toml file for the provided resource files is config/resource.toml.
This software is written in Python 3.12. The required packages are in pyproject.toml file. The users could install the required packages using the following command using hatchling backend:
# optional: install and upgrade build setuptools and whell
pip install --upgrade build pip setuptools wheel
git clone https://github.com/zhuweix/methylFracAnalyzer.git
cd methylFracAnalyzer
# build the package
python -m build
# install the package
pip install .
The details of the resource files provided with the software are as follows:
Cpgislands, enhancers, silencers, tRNA genes, origin of replications and centromeres in Hg38 assemlby (gencode v43 annotation): resource/hg38_gencode_v43.other_feature.csv.gz
We used the protein-coding genes in Hg38 assemlby (gencode v43 annotation). For genes with multiple transcripts we chose the best transcript based on expression level. The expression level were calculated using the RNA-seq data in MCF7 and MCF10A cell lines (GEO accession: GSE201262 for MCF7 and GSE237066 for MCF10A)[1,2]. The normalized read count is calculated by Salmon [3] and the transcript with the highest expression level are chosen. For genes with no expression, the primary transcript is chosen. The expressed genes were grouped in 5 quantiles based on expression level. Furthermore, the genes are also grouped be active NDR and inactive NDR using ATAC-seq (GEO accession: GSE201262 for MCF7 and GSE152410 for MCF10A) [1,4]. The ATAC-seq data were normalized by the non-zero read coverage in each chromosome, and we used the average ATAC-seq data in 2 replicates of MCF7 and 4 replicates in MCF10A. We used the average normalized ATAC-seq signal in +/- 1 kb region relative to the TSS of the genes. The active genes in MCF7 and MCF10A are filtered based on ATAC-seq signal.
The list of used annotation files:
Genes in MCF7: resource/MCF7_RNAseq_ATAC.csv
Genes in MCF10A: resource/MCF10_RNAseq_ATAC.csv
The centromeric elements in T2T v1.1 assembly are obtained from the Altemose et al[5].
List of used annotation files:
Entire centromere location (all the elements excluding peri-centromeric regions): resource/centromere_location.csv
Location of each centromeric element: resource/centromere_element.csv
Location of active alpha satellites: resource/active_hor.csv
Location of CENPA regions: resource/cenpA_location.csv
We aligned the raw ChIP-seq data to T2T v1.1 assembly using chip-seq-pipeline2 and predicted optimal and conservative peaks for Histone Marks. We used the ChromHMM software to predict 15 ChromHMM states based on the aligned BAM files [6]. The states are renamed based on the combination the histone marks and enrichment in genomic locations in a similar fashion with Ernst et al [6]. The GEO accession number of ChIP-seq data is GSE85158 for MCF7 and MCF10A [7].
Optimal peaks for Histone Marks:
MCF7: resource/MCF7_HistoneMark.optimal.csv.gz
MCF10A: resource/MCF10_HistoneMark.optimal.csv.gz
Conservative peaks for Histone Marks:
MCF7: resource/MCF7_HistoneMark.conservative.csv.gz
MCF10A: resource/MCF10_HistoneMark.conservative.csv.gz
ChromHMM states:
MCF7: resource/MCF7_ChromHMM.t2t.csv.gz
MCf10A: resource/MCF10_ChromHMM.t2t.csv.gz
The MNase-seq data are used to plot the nucleosome phasing relative to CTCF sites and TSS of protein-coding genes. The MNase-seq data are available in GEO (TBD). The MNase-seq data were aligned to Hg38 genome and the single nucleosome fragments were filtered be size=120 to 180 bp. The nucleosome dyads were counted by the center of the fragments. For each gene transcript specified in the Protein-coding genes section, the nucleosome dyads in +/- 2010 bp of the TSS were counted, and the dyad count were normalized by the average dyad count in the flanking regions. The average dyad count of the active and inactive genes were further smoothed in 21-bp windows using moving average. The smoothed dyad count relative to CTCF sites were also calculated in the samilar way. MCF10A and MCF7 cells were similar cell lines and we currently also used the MNase data in MCF7 to show the nucleosome phasing in MCF10A cells as well.
List of used annotation files:
Gene: resource/MCF7_MNase.gene_average_phasing.csv CTCF: resource/MCF7_MNase.CTCF_average_phasing.csv
Please cite the following paper if you use this software: TBD
[1] Tian, C.. Impaired histone inheritance promotes tumor progression. Nature Communications 14, (2023).
[2] Dorgham, M. G., Elliott, B. A., Holley, C. L. & Mansfield, K. D.. m6A regulates breast cancer proliferation and migration through stage-dependent changes in Epithelial to Mesenchymal Transition gene expression. Frontiers in Oncology 13, (2023).
[3] Patro, R., Duggal, G., Love, M. I., Irizarry, R. A. & Kingsford, C.. Salmon provides fast and bias-aware quantification of transcript expression. Nature Methods 14, 417–419 (2017).
[4] Gross, S. M.. A multi-omic analysis of MCF10A cells provides a resource for integrative assessment of ligand-mediated molecular and phenotypic responses. Communications Biology 5, (2022).
[5] Altemose, N.. Complete genomic and epigenetic maps of human centromeres. Science 376, (2022).
[6] Ernst, J. & Kellis, M.. Chromatin-state discovery and genome annotation with ChromHMM. Nature Protocols 12, 2478–2492 (2017).
[7] Franco, H. L.. Enhancer transcription reveals subtype-specific gene expression programs controlling breast cancer pathogenesis. Genome Research 28, 159–170 (2018).