Comparative analysis of genome-scale, base-resolution DNA methylation profiles across 580 animal species.
Methylation of cytosines is the prototypic epigenetic modification of the DNA. It has been implicated in various regulatory mechanisms throughout the animal kingdom and particularly in vertebrates. We mapped DNA methylation in 580 animal species (535 vertebrates, 45 invertebrates), resulting in 2443 genome-scale DNA methylation profiles of multiple organs. Reference-genome independent analysis of this comprehensive dataset quantified the association of DNA methylation with the underlying genomic DNA sequence throughout vertebrate evolution. We observed a broadly conserved link with two major transitions – once in the first vertebrates and again with the emergence of reptiles. Cross-species com-parisons focusing on individual organs supported a deeply conserved association of DNA methylation with tissue type, and cross-mapping analysis of DNA methylation at gene promoters revealed evolu-tionary changes for orthologous genes with conserved DNA methylation patterns. In summary, this study establishes a large resource of vertebrate and invertebrate DNA methylomes, it showcases the power of reference-free epigenome analysis in species for which no reference genomes are available, and it contributes an epigenetic perspective to the study of vertebrate evolution.
More supplementary data and interactive data exploration are available on the supplementary website.
This repository provides a set of scripts used to analyze the data. Below please find step-by-step instructions on how to reproduce the analysis in the paper. There are two categories of scripts: those that have to be run in HPC mode on each species in parallel, and those used for summarizing and visualizing results by generating summary tables and figure panels. The provided HPC scripts are written for a SLURM scheduling system and would have to be adjusted for your own computational environment and scheduler.
| Panel | Script |
|---|---|
| 1A | Overview, no script |
| 1B | 02.3_stats_detail.R |
| 1C | 02.6_create_ITOL.R |
| 1D | 02.3_stats_detail.R |
| 1E | 02.5_stats_methylation.R |
| 1F | 03.1_predict_average_meth_levels.R & 03.2_kmer_stability.R |
| 1G | 04.1_phylogeny_reconstruction.R |
| 1H | 03.3_phyloglm.R |
| 1I | 02.7_speciesMeta.R |
| 1J | 02.7_speciesMeta.R |
| 1K | 02.7_speciesMeta.R |
| 2A | Overview, no script |
| 2B | 05.2_summarize_within_sp_prediction.R |
| 2C | 05.2_summarize_within_sp_prediction.R |
| 2D | 05.2_summarize_within_sp_prediction.R |
| 2E | 05.3_feature_weight_analysis.R |
| 2F | 05.3_feature_weight_analysis.R |
| 3A | 05.5_summary_cross_species.R |
| 3B | 05.4_predict_crossSp.R |
| 3C | 05.5_summary_cross_species.R |
| 3D | 06.1_inverted.R |
| 3E | 06.4_explore_inverted_final_WHH.R |
| 3F | 06.4_explore_inverted_final_WHH.R |
| 4A | 08.2_tissue_correlation.R |
| 4B | Overview, no script |
| 4C | 07.5_ame_selex.R |
| 4D | 07.8_followup_analysis.ipynb & cytoscape |
| 4E | 07.8_followup_analysis.ipynb & cytoscape |
| 5A | validation/01.5.3_crossMapping_plot.R |
| 5B | validation/01.5.3_crossMapping_plot.R |
| 5C | validation/01.5.3_crossMapping_plot.R |
| 5D | validation/01.5.3_crossMapping_plot.R |
| 5E | validation/01.5.3_crossMapping_plot.R |
For computational analysis, we use R version 3.6.1 with all the corresponding packages. The SVM module requires an independent block with R version 4.1.
We provide a full list of packages installed in each of the subsets with their versions in the envs/ folder as csv files. For your convinience, we also provide the mirrors of conda environments as yaml files, generated with conda env export. Please note that not all packages were installed within conda, so additionally you would still have to install missing R packages from the csv tables if needed for the corresponding part of the analysis.
The project folder contains the following folders: -src -meta -raw_data (can be downloaded from SRA SRP357738) -results_pipeline -results_analysis
The results_pipeline folder should contain the per-species output of the RefFreeDMA. Most relevant files (deduced reference genome and methylation profiles) are shared and could be downloaded from the GEO repository under accession number GSE195869 or from the supplementary website.
The results_analysis folder contains the output of most of the scripts below. We also have provided the most relevant summary tables for download on the supplementary website.
00.0_init.R is the script, to initialize all further scripts, it will standardize settings, make annotations available, and set paths to data and scripts. A few absolute paths need to be adjusted and two environmental variables are required (they have to be visible by both Rscript and bash):
CODEBASE (location of this github repository)
LOGDIR (location for pipeline log and tmp files)
For data processing you would have to use RefFreeDMA and bisulfiteBlast. The 00.1 script will prepare the run folders and shell scripts containing the commands to run RefFeeDMA. For setting up RefFreeDMA please follow the instructions in the corresponding repository.
Rscript 00.1_prepareRefFreeDMA.R Following scripts will generate some additional data metrices, that are used in most subsequent analysis. Here we create the big summary annotation table and other initial stats:
Rscript 01.1_collectStats.R
Rscript 01.2_taxonomicAnnotation.R. #generates the stats annot table
Rscript 01.3_quality_stratification.RCalcualating the frequency of CGs for subsequent analysis:
Rscript 01.4_dedRef_CpG_count.R ##takes a while
Calculating the frequency of kmers, using the MEME tool: !!Importantly, the path in 01.5.sbatch has to be adjusted
sbatch sbatch/01.5.1.sbatch ##prefilteing
sbatch sbatch/01.5.2.sbatch
Rscript 01.5.3_kmer_count_summary.R ##summary for further useGeneration of the phylogenetic tree for iTOL, it's optimization for future analysis and exporting the taxonomimc order of species for use in subsequent analysis:
Rscript 01.6_create_ITOL.R
Rscript 01.7_species_order.RCalculating the CpG island properties:
sbatch sbatch/01.8.sbatch Coverage statistics and sequence features:
Rscript 01.9_coverage_analysis.R
Rscript 01.10_seq_features.RThe following scripts create summary statistics and generate summary plots, giving us an overview of the data:
Rscript 02.1_stats_uc.R
Rscript 02.2_stats_general.R
Rscript 02.3_stats_detail.R
Rscript 02.4_stats_location.R
Rscript 02.5_stats_methylation.R
Rscript 02.7_speicesMeta.RSaving the supplementary tables and the stats for the paper:
Rscript 02.6_save_annot.RStep-by-step exploration of features, that explain mean DNA methylation level per species:
Rscript 03.1_predict_average_meth_levels.R ##creating a table of features & prediction values for each group (takes a couple of hours on a large dataset)
Rscript 03.2_kmer3_stability.R ##stepwise selection on 3mers
Rscript 03.3_phyloglm.R ##phylogenetic modelsUsing kmer frequencies to create a hierarchical clustering of species, based on the sequence composition of the consensus genomes:
Rscript 04.1_phylogeny_reconstruction.RWith-in species DNA methylation prediction:
sbatch sbatch/05.1.sbatch ##repeated on k 1:10, k = 3, k = 4, k = 7
Rscript 05.2_summarize_within_sp_prediction.R
Rscript 05.3_feature_weights_analysis.RAdditional check with more stringent thresholds:
sbatch sbatch/05.1.1.sbatch
Rscript 05.2.1_summary_alt.RCross-species prediction:
sbatch sbatch/05.4.sbatch
Rscript 05.5_summary_cross_species_prediction.RIdentifying the species, that we classify as "inverted":
Rscript 06.1_inverted.RConfirm the observed effect with additional experiments on subsets of data (by tissue, individual, and random sub-sampling):
sbatch sbatch/06.2.1.sbatch <candidate species>
sbatch sbatch/06.2.2.sbatch <candidate species>
sbatch sbatch/06.2.3.sbatch <candidate species>
Rscript 06.2_verify_inv_summary.RAnalysing the frequencies of kmers/repeats specifically in highly/lowly methylated sequences: !!first two python scripts require additional specification of the data_dir path
python 06.3.1_to_fasta_for_kmer_count.py
sbatch 06.3.2.sbatch #runs on all the species, where classes contain inverted
RScript 06.3_freq_kmers.RAnd summarizing all the knowledge, oriented on the WHH species
Rscript 06.4_explore_inverted_final_WHH.R
As a preparation, the JASPAR matrix of all the TF profiles has to be downloaded: http://jaspar.genereg.net/downloads/ We used the non-redundant vertebrate TFBS set (additionally download the JASPAR annotation table) Calculate and report frequencies of each TFBS in consensus reference of each species:
sbatch sbatch/07.1.sbatch
sbatch sbatch/07.2.sbatch
Rscript 07.3_TFBS_summary.RIdentification of differentially methylated fragemnts between heart and liver tissues and motif enrichment analysis: (make sure you have ame installed)
sbatch sbatch/07.4.1.sbatch
sbatch sbatch/07.4.2.sbatch
Rscript 07.4.ame_summary.RFollow-up on enriched TFBS: For integrating the methylation preference for TF binding, we used the data from Yin et. al, 2017.
Rscript 07.5_ame_selex.R
Rscript 07.6_GO_enrichment.R
Rscript 07.7_GRN.RInteractively: 07.8_followup_analysis.ipynb
Multidimensional scaling analysis for selected species:
Rscript 08.1_tissue_MDS.RComparison of DNA methylation variability due to tissue and individual:
Rscript 08.2_tissue_correlation.R
Rscript 08.3_tissue_distance.RCreate correlation plots between liver and heart DNA methylattion levels for selected species:
08.4_tissue_diffMeth.Rdetermining which genomes to process, cleaning them up and concatting:
Rscript 01.1_find_mapping_match.R ## finding pairs of genomes to match
Rscript 01.1.2_axolotl.R ## rearrange and remove small chromosomes of the genome of the axolotl
sbatch 01.2_concat_genome.sbatch ## running genome concatThe download procedure is described downstream in 02.1_insilico_digest_prepare_genomes_and_tracks.ipynb. Generating the config files and running the cross-mapping:
Rscript 01.3_save_configs.Rafterwards run the job-subbmitting bash script same way as for original RefFreeDMA.
Summarise the results of crossMapping:
Rscript 01.4_crossMapping_stats.R
sbatch sbatch/01.5_cross-mapping_analysis.sbatch
Rscript 01.5.2_crossMappingAnalysis_summary.R
Rscript 01.5.3_crossMapping_plot.RIn interactive jupyer notebook run the preparation notebook for genomes: 02.1_insilico_digest_prepare_genomes.R
sbatch 02.2_insilicoDigest.sbatch
Rscript 02.3_insilicoDigest_summary.RWGBS data download is described at 03.1_download_WGBS.ipynb:
Rscript 03.2.process_WGBS_data.R
Rscript 03.2.2_compare_mean_meth.R
sbatch 03.3_WGBS_fragment_CpGmeth.sbatch ##preparing the datasbatch 03.4.sbatch
sbatch 03.5.sbatch
sbatch 03.6.sbatch
Rscript 03.4.1_WGBS_prediction_summary.R
Rscript 03.5_summary.R
Rscript 03.6_summary.R