--- title: "12-Peve-WGBS-bismark.Rmd" author: "Zoe Dellaert" date: "2025-04-10" output: github_document: toc: true number_sections: true bibliography: references.bib --- ```{r setup, include=FALSE} library(knitr) knitr::opts_chunk$set( echo = TRUE, # Display code chunks eval = TRUE, # Evaluate code chunks warning = FALSE, # Hide warnings message = FALSE, # Hide messages comment = "" # Prevents appending '##' to beginning of lines in code output ) ``` ## This is the downstream methylation analysis of the WGBS data for *Porites evermanni* Reads were trimmed and QC'd in [this code](https://github.com/urol-e5/deep-dive-expression/blob/main/E-Peve/code/01.00-E-Peve-WGBS-trimming-cutadapt-FastQC-MultiQC.md) ## Important file locations: 1. [Trimmed WGBS Reads](https://gannet.fish.washington.edu/gitrepos/urol-e5/deep-dive-expression/E-Peve/output/01.00-E-Peve-WGBS-trimming-cutadapt-FastQC-MultiQC/) 2. [Bismark Genome](https://gannet.fish.washington.edu/gitrepos/urol-e5/deep-dive-expression/E-Peve/data/Bisulfite_Genome/) 3. [All Bismark output files (BAMs, .cov files, .bedgraph files)](https://gannet.fish.washington.edu/gitrepos/urol-e5/deep-dive-expression/E-Peve/output/12-Peve-WGBS/bismark_cutadapt/) 4. [Bismark and Qualimap MultiQC](https://gannet.fish.washington.edu/gitrepos/urol-e5/deep-dive-expression/E-Peve/output/12-Peve-WGBS/bismark_cutadapt/multiqc_report.html) and [base bismark report](https://gannet.fish.washington.edu/gitrepos/urol-e5/deep-dive-expression/E-Peve/output/12-Peve-WGBS/bismark_cutadapt/bismark_summary_report.html) ### Note: Most of this code is based on the [E5 Time Series Molecular](https://github.com/urol-e5/timeseries_molecular) code by Steven Roberts [here](https://github.com/urol-e5/timeseries_molecular/blob/main/D-Apul/code/15.5-Apul-bismark.qmd) ## Generate Bismark Bisulfite Genome ```{bash, eval=FALSE} #!/usr/bin/env bash #SBATCH --export=NONE #SBATCH --nodes=1 --ntasks-per-node=20 #SBATCH --mem=200GB #SBATCH -t 24:00:00 #SBATCH --mail-type=BEGIN,END,FAIL #email you when job starts, stops and/or fails #SBATCH --error=scripts/outs_errs/"%x_error.%j" #if your job fails, the error report will be put in this file #SBATCH --output=scripts/outs_errs/"%x_output.%j" #once your job is completed, any final job report comments will be put in this file # load modules needed module load uri/main module load Bismark/0.23.1-foss-2021b module load bowtie2/2.5.2 cd ../data bismark_genome_preparation --verbose --parallel 10 ./ ``` ### output: ```{bash, eval=FALSE} Using 10 threads for the top and bottom strand indexing processes each, so using 20 cores in total Writing bisulfite genomes out into a single MFA (multi FastA) file Bisulfite Genome Indexer version v0.23.1 (last modified: 27 Jan 2021) Step I - Prepare genome folders - completed Step II - Genome bisulfite conversions - completed Bismark Genome Preparation - Step III: Launching the Bowtie 2 indexer Preparing indexing of CT converted genome in /scratch3/workspace/zdellaert_uri_edu-deep_dive/deep-dive-expression/E-Peve/data/Bisulfite_Genome/CT_conversion/ Building a SMALL index Preparing indexing of GA converted genome in /scratch3/workspace/zdellaert_uri_edu-deep_dive/deep-dive-expression/E-Peve/data/Bisulfite_Genome/GA_conversion/ Building a SMALL index Renaming BS_GA.3.bt2.tmp to BS_GA.3.bt2 Renaming BS_GA.4.bt2.tmp to BS_GA.4.bt2 Renaming BS_GA.1.bt2.tmp to BS_GA.1.bt2 Renaming BS_GA.2.bt2.tmp to BS_GA.2.bt2 Renaming BS_GA.rev.1.bt2.tmp to BS_GA.rev.1.bt2 Renaming BS_GA.rev.2.bt2.tmp to BS_GA.rev.2.bt2 Renaming BS_CT.3.bt2.tmp to BS_CT.3.bt2 Renaming BS_CT.4.bt2.tmp to BS_CT.4.bt2 Renaming BS_CT.1.bt2.tmp to BS_CT.1.bt2 Renaming BS_CT.2.bt2.tmp to BS_CT.2.bt2 Renaming BS_CT.rev.1.bt2.tmp to BS_CT.rev.1.bt2 Renaming BS_CT.rev.2.bt2.tmp to BS_CT.rev.2.bt2 ``` ### Compress and generate md5 ```{bash, eval=FALSE} cd ../data tar -czvf Bisulfite_Genome.tar.gz Bisulfite_Genome md5sum Bisulfite_Genome.tar.gz | tee Bisulfite_Genome.tar.gz.md5 ``` ```{bash, eval=FALSE} 885c84c233ca313f23b8670903a18db6 Bisulfite_Genome.tar.gz ``` ### Output file location: [Bismark Genome](https://gannet.fish.washington.edu/gitrepos/urol-e5/deep-dive-expression/E-Peve/data/Bisulfite_Genome/) ## Test parameters ```{bash, eval=FALSE} #!/usr/bin/env bash #SBATCH --ntasks=1 --cpus-per-task=30 #split one task over multiple CPU #SBATCH --array=0-4 #for 5 samples #SBATCH --mem=100GB #SBATCH -t 00:30:00 #SBATCH --mail-type=END,FAIL,TIME_LIMIT_80 #email you when job stops and/or fails or is nearing its time limit #SBATCH --error=scripts/outs_errs/"%x_error.%j" #if your job fails, the error report will be put in this file #SBATCH --output=scripts/outs_errs/"%x_output.%j" #once your job is completed, any final job report comments will be put in this file # load modules needed module load uri/main module load Bismark/0.23.1-foss-2021b module load bowtie2/2.5.2 # Set directories and files reads_dir="../output/01.00-E-Peve-WGBS-trimming-cutadapt-FastQC-MultiQC/" genome_folder="../data/" output_dir="../output/12-Peve-WGBS/bismark_paramtest_cutadapt" checkpoint_file="${output_dir}/completed_samples.log" # make output directory mkdir -p ${output_dir} # Create the checkpoint file if it doesn't exist touch ${checkpoint_file} # Get the list of sample files and corresponding sample names files=(${reads_dir}*_R1_001.fastq.gz) file="${files[$SLURM_ARRAY_TASK_ID]}" sample_name=$(basename "$file" "_R1_001.fastq.gz") # Check if the sample has already been processed if grep -q "^${sample_name}$" ${checkpoint_file}; then echo "Sample ${sample_name} already processed. Skipping..." exit 0 fi # Define log files for stdout and stderr stdout_log="${output_dir}/${sample_name}_stdout.log" stderr_log="${output_dir}/${sample_name}_stderr.log" # Define the array of score_min parameters to test score_min_params=( "L,0,-0.4" "L,0,-0.6" "L,0,-0.8" "L,0,-1.0" "L,-1,-0.6" ) # Loop through each score_min parameter for score_min in "${score_min_params[@]}"; do echo "Running Bismark for sample ${sample_name} with score_min ${score_min}" # Create a subdirectory for this parameter param_output_dir="${output_dir}/${sample_name}_score_${score_min//,/}" mkdir -p ${param_output_dir} # Run Bismark alignment bismark \ -genome ${genome_folder} \ -p 8 \ -u 10000 \ -score_min ${score_min} \ --non_directional \ -1 ${reads_dir}${sample_name}_R1_001.fastq.gz \ -2 ${reads_dir}${sample_name}_R2_001.fastq.gz \ -o ${param_output_dir} \ --basename ${sample_name}_${score_min//,/} \ 2> "${param_output_dir}/${sample_name}-bismark_summary.txt" # Check if the command was successful if [ $? -eq 0 ]; then echo "Sample ${sample_name} with score_min ${score_min} processed successfully." else echo "Sample ${sample_name} with score_min ${score_min} failed. Check ${stderr_log} for details." fi done # Mark the sample as completed in the checkpoint file if [ $? -eq 0 ]; then echo ${sample_name} >> ${checkpoint_file} echo "All tests for sample ${sample_name} completed." else echo "Sample ${sample_name} encountered errors. Check logs for details." fi # Define directories summary_file="${output_dir}/parameter_comparison_summary.csv" # Initialize summary file echo "Sample,Score_Min,Alignment_Rate" > ${summary_file} # Loop through parameter output directories for dir in ${output_dir}/*_score_*; do if [ -d "$dir" ]; then # Extract sample name and score_min parameter from directory name sample_name=$(basename "$dir" | cut -d '_' -f1-3) score_min=$(basename "$dir" | grep -o "score_.*" | sed 's/score_//; s/_/,/g') # Locate the summary file summary_file_path="${dir}/${sample_name}_${score_min}_PE_report.txt" # Extract metrics mapping=$(grep "Mapping efficiency:" ${summary_file_path} | awk '{print "mapping efficiency ", $3}') # Append to the summary file echo "${sample_name},${score_min},${mapping}" >> ${summary_file} fi done ``` ### Results from parameter tests: | Sample | Score_Min | Alignment_Rate | |-----------------|-----------|----------------| | trimmed_421_S8 | L0-0.4 | 29.20% | | trimmed_421_S8 | L0-0.6 | 36.00% | | trimmed_421_S8 | L0-0.8 | 41.00% | | trimmed_421_S8 | L0-1.0 | 47.30% | | trimmed_421_S8 | L-1-0.6 | 36.30% | | trimmed_471_S6 | L0-0.4 | 35.90% | | trimmed_471_S6 | L0-0.6 | 42.80% | | trimmed_471_S6 | L0-0.8 | 48.40% | | trimmed_471_S6 | L0-1.0 | 54.70% | | trimmed_471_S6 | L-1-0.6 | 43.00% | | trimmed_487_S9 | L0-0.4 | 34.60% | | trimmed_487_S9 | L0-0.6 | 41.40% | | trimmed_487_S9 | L0-0.8 | 47.50% | | trimmed_487_S9 | L0-1.0 | 54.20% | | trimmed_487_S9 | L-1-0.6 | 41.80% | | trimmed_489_S10 | L0-0.4 | 32.50% | | trimmed_489_S10 | L0-0.6 | 40.30% | | trimmed_489_S10 | L0-0.8 | 47.00% | | trimmed_489_S10 | L0-1.0 | 54.20% | | trimmed_489_S10 | L-1-0.6 | 40.70% | | trimmed_491_S7 | L0-0.4 | 42.00% | | trimmed_491_S7 | L0-0.6 | 50.50% | | trimmed_491_S7 | L0-0.8 | 56.60% | | trimmed_491_S7 | L0-1.0 | 62.10% | | trimmed_491_S7 | L-1-0.6 | 50.80% | I ran the parameter testing before fixing the file sample names. For reference: - 421 = POR-82 - 471 = POR-76 - 487 = POR-71 - 489 = POR-79 - 491 = POR-73 ## Align to genome ```{bash, eval=FALSE} #!/usr/bin/env bash #SBATCH --ntasks=1 --cpus-per-task=48 #split one task over multiple CPU #SBATCH --array=0-4 #for 5 samples #SBATCH --mem=400GB #SBATCH -t 48:00:00 #SBATCH --mail-type=END,FAIL,TIME_LIMIT_80 #email you when job stops and/or fails or is nearing its time limit #SBATCH --error=scripts/outs_errs/"%x_error.%j" #if your job fails, the error report will be put in this file #SBATCH --output=scripts/outs_errs/"%x_output.%j" #once your job is completed, any final job report comments will be put in this file # load modules needed module load uri/main module load Bismark/0.23.1-foss-2021b module load bowtie2/2.5.2 # Set directories and files reads_dir="../output/01.00-E-Peve-WGBS-trimming-cutadapt-FastQC-MultiQC/" genome_folder="../data/" output_dir="../output/12-Peve-WGBS/bismark_cutadapt" checkpoint_file="${output_dir}/completed_samples.log" # make output directory mkdir -p ${output_dir} # Create the checkpoint file if it doesn't exist touch ${checkpoint_file} # Get the list of sample files and corresponding sample names files=(${reads_dir}*_R1_001.fastq.gz) file="${files[$SLURM_ARRAY_TASK_ID]}" sample_name=$(basename "$file" "_R1_001.fastq.gz") # Check if the sample has already been processed if grep -q "^${sample_name}$" ${checkpoint_file}; then echo "Sample ${sample_name} already processed. Skipping..." exit 0 fi # Define log files for stdout and stderr stdout_log="${output_dir}/${sample_name}_stdout.log" stderr_log="${output_dir}/${sample_name}_stderr.log" # Run Bismark alignment bismark \ -genome ${genome_folder} \ -p 48 \ -score_min L,0,-1.0 \ --non_directional \ -1 ${reads_dir}${sample_name}_R1_001.fastq.gz \ -2 ${reads_dir}${sample_name}_R2_001.fastq.gz \ -o ${output_dir} \ --basename ${sample_name} \ 2> "${output_dir}/${sample_name}-bismark_summary.txt" # Check if the command was successful if [ $? -eq 0 ]; then # Append the sample name to the checkpoint file echo ${sample_name} >> ${checkpoint_file} echo "Sample ${sample_name} processed successfully." else echo "Sample ${sample_name} failed. Check ${stderr_log} for details." fi # Define directories summary_file="${output_dir}/alignment_summary.csv" # Initialize summary file echo "Sample,Score_Min,Alignment_Rate" > ${summary_file} # Loop through parameter output directories for file in ${output_dir}/*_report.txt; do # Extract sample name and score_min parameter from directory name sample_name=$(basename "$file" | cut -d'_' -f1-3) score_min="L0-1.0" # Locate the summary file summary_file_path="${output_dir}/${sample_name}_PE_report.txt" # Extract metrics mapping=$(grep "Mapping efficiency:" ${summary_file_path} | awk '{gsub("%", "", $3); print $3}') # Append to the summary file echo "${sample_name},${score_min},${mapping}" >> ${summary_file} done ``` ### Output file location: [All Bismark output files](https://gannet.fish.washington.edu/gitrepos/urol-e5/deep-dive-expression/E-Peve/output/12-Peve-WGBS/bismark_cutadapt/) ## Post-alignment code is based once again on [Steven's code](https://github.com/urol-e5/timeseries_molecular/blob/main/D-Apul/code/15.5-Apul-bismark.qmd) ### Deduplication, Sorting, and methylation extraction & calling ```{bash, eval=FALSE} #!/usr/bin/env bash #SBATCH --export=NONE #SBATCH --nodes=1 --ntasks-per-node=24 #SBATCH --mem=250GB #SBATCH -t 24:00:00 #SBATCH --mail-type=BEGIN,END,FAIL #email you when job starts, stops and/or fails #SBATCH --error=scripts/outs_errs/"%x_error.%j" #if your job fails, the error report will be put in this file #SBATCH --output=scripts/outs_errs/"%x_output.%j" #once your job is completed, any final job report comments will be put in this file # load modules needed module load uri/main module load Bismark/0.23.1-foss-2021b module load parallel/20240822 # set directories bismark_dir="../output/12-Peve-WGBS/bismark_cutadapt/" genome_folder="../data/" ### deduplicate bams find ${bismark_dir}*.bam | \ xargs -n 1 basename | \ sed 's/^trimmed_//' | sed 's/_pe.bam$//' | \ parallel -j 8 deduplicate_bismark \ --bam \ --paired \ --output_dir ${bismark_dir} \ ${bismark_dir}trimmed_{}_pe.bam ### methylation extraction find ${bismark_dir}*deduplicated.bam | xargs -n 1 -I{} \ bismark_methylation_extractor --bedGraph --counts --comprehensive --merge_non_CpG \ --multicore 24 --buffer_size 75% --output ${bismark_dir} "{}" ### methylation call find ${bismark_dir}*deduplicated.bismark.cov.gz | \ xargs -n 1 basename | \ sed 's/^trimmed_//' | sed 's/_pe.deduplicated.bismark.cov.gz$//' | \ parallel -j 24 coverage2cytosine \ --genome_folder ${genome_folder} \ -o ${bismark_dir}{} \ --merge_CpG \ --zero_based \ ${bismark_dir}trimmed_{}_pe.deduplicated.bismark.cov.gz ### sort bams # change modules module purge module load samtools/1.19.2 find ${bismark_dir}*deduplicated.bam | \ xargs -n 1 basename | \ sed 's/^trimmed_//' | sed 's/_pe.deduplicated.bam$//' | \ xargs -I{} samtools \ sort --threads 24 \ ${bismark_dir}trimmed_{}_pe.deduplicated.bam \ -o ${bismark_dir}{}.sorted.bam ``` This took just over 6 hours with max memory used per node as 249.99GiB. ### View output ```{bash, eval=FALSE} head ${bismark_dir}*evidence.cov ``` ### Make summary reports ```{bash, eval=FALSE} #!/usr/bin/env bash #SBATCH --export=NONE #SBATCH --nodes=1 --ntasks-per-node=8 #SBATCH --mem=250GB #SBATCH -t 08:00:00 #SBATCH --mail-type=BEGIN,END,FAIL #email you when job starts, stops and/or fails #SBATCH --error=scripts/outs_errs/"%x_error.%j" #if your job fails, the error report will be put in this file #SBATCH --output=scripts/outs_errs/"%x_output.%j" #once your job is completed, any final job report comments will be put in this file module load uri/main module load Bismark/0.23.1-foss-2021b module load all/MultiQC/1.12-foss-2021b module load qualimap/2.2.1 cd ../output/12-Peve-WGBS/bismark_cutadapt/ bam2nuc --genome_folder ../../../data/ *_pe.deduplicated.bam mkdir -p qualimap/bamqc for bamFile in *sorted.bam; do prefix=$(basename $bamFile .bam) qualimap \ --java-mem-size=29491M \ bamqc \ \ -bam ${bamFile} \ \ -p non-strand-specific \ --collect-overlap-pairs \ -outdir qualimap/bamqc/${prefix} \ -nt 6 done bismark2report bismark2summary *pe.bam multiqc . ``` #### Output file location: [Bismark and Qualimap MultiQC](https://gannet.fish.washington.edu/gitrepos/urol-e5/deep-dive-expression/E-Peve/output/12-Peve-WGBS/bismark_cutadapt/multiqc_report.html) and [base bismark report](https://gannet.fish.washington.edu/gitrepos/urol-e5/deep-dive-expression/E-Peve/output/12-Peve-WGBS/bismark_cutadapt/bismark_summary_report.html) ### Sorting cov files and filtering for coverage and gene intersection ```{bash, eval=FALSE} # salloc -p cpu -c 8 --mem 16G # load modules needed module load bedtools2/2.31.1 # set directories and files root_dir="/scratch3/workspace/zdellaert_uri_edu-deep_dive/deep-dive-expression/E-Peve/" bismark_dir="${root_dir}/output/12-Peve-WGBS/bismark_cutadapt/" genome_folder="${root_dir}/data/" gtf_name="Porites_evermanni_validated" # Sort .cov files cd ${bismark_dir} for file in *merged_CpG_evidence.cov do sample=$(basename "${file}" .CpG_report.merged_CpG_evidence.cov) bedtools sort -i "${file}" \ > "${sample}"_sorted.cov done # Create bedgraphs for 5X coverage for file in *_sorted.cov do sample=$(basename "${file}" _sorted.cov) cat "${file}" | awk -F $'\t' 'BEGIN {OFS = FS} {if ($5+$6 >= 5) {print $1, $2, $3, $4}}' \ > "${sample}"_5x_sorted.bedgraph done # Create .tab files for 5x coverage for file in *_sorted.cov do sample=$(basename "${file}" _sorted.cov) cat "${file}" | awk -F $'\t' 'BEGIN {OFS = FS} {if ($5+$6 >= 5) {print $1, $2, $3, $4, $5, $6}}' \ > "${sample}"_5x_sorted.tab done # Intersection of all samples multiIntersectBed -i *_5x_sorted.tab > CpG.all.samps.5x_sorted.bed ## change number after == to your number of samples cat CpG.all.samps.5x_sorted.bed | awk '$4 ==5' > CpG.filt.all.samps.5x_sorted.bed # Intersection of all samples with gene bodies: awk '{if ($3 == "transcript") {print}}' "${genome_folder}/${gtf_name}.gtf" > "${genome_folder}/${gtf_name}_transcripts.gtf" for i in *5x_sorted.tab do intersectBed \ -wb \ -a ${i} \ -b "${genome_folder}/${gtf_name}_transcripts.gtf" \ > ${i}_gene done # Keep only loci intersecting with genes found in all samples for i in *_5x_sorted.tab_gene do intersectBed \ -a ${i} \ -b CpG.filt.all.samps.5x_sorted.bed \ > ${i}_CpG_5x_enrichment.bed done # Global methylation levels for file in *_sorted.cov; do sample=$(basename "$file" _sorted.cov) awk '{methylated+=$5; unmethylated+=$6} END {print "'$sample'", methylated/(methylated+unmethylated)}' "$file" done > global_methylation_levels.txt ``` ### Output file location: [All Bismark output files](https://gannet.fish.washington.edu/gitrepos/urol-e5/deep-dive-expression/E-Peve/output/12-Peve-WGBS/bismark_cutadapt/) ## Methylkit ```{r} # Load methylKit and other packages library("methylKit") library("tidyverse") library("parallel") file_list <- list.files("/scratch3/workspace/zdellaert_uri_edu-deep_dive_exp/deep-dive-expression/E-Peve/output/12-Peve-WGBS/bismark_cutadapt/",pattern = ".merged_CpG_evidence.cov$", full.names = TRUE, include.dirs = FALSE) sample <- gsub("_S\\d{1,2}.CpG_report.merged_CpG_evidence.cov", "", basename(file_list)) file_list <- list.files("/scratch3/workspace/zdellaert_uri_edu-deep_dive_exp/deep-dive-expression/E-Peve/output/12-Peve-WGBS/bismark_cutadapt/",pattern = "_pe.deduplicated.bismark.cov.gz", full.names = TRUE, include.dirs = FALSE) sample <- gsub("trimmed_", "", basename(file_list)) sample <- gsub("_S\\d{1,2}_pe.deduplicated.bismark.cov.gz", "", sample) file_list <- as.list(file_list) treatment <- c(0,0,0,0,0) sample file_list # methyl_list <- mclapply(seq_along(file_list), function(i) { # methRead( # file_list[[i]], # assembly = "Peve", # treatment = treatment[i], # sample.id = sample[[i]], # context = "CpG", # mincov = 10, # pipeline = "bismarkCoverage" # ) # }, mc.cores = 4) # # methylObj <- new("methylRawList", methyl_list) # methylObj@treatment <- c(0, 0, 0, 0,0) # # save(methylObj, file = "../output/12-Peve-WGBS/methylkit/MethylObj.RData") ``` ```{r, eval=FALSE} load("../output/12-Peve-WGBS/methylkit/MethylObj.RData") getMethylationStats(methylObj[[1]],plot=FALSE,both.strands=FALSE) getMethylationStats(methylObj[[1]], plot = TRUE,both.strands=FALSE) getCoverageStats(methylObj[[1]],plot=TRUE,both.strands=FALSE) ``` ```{r,eval=FALSE} filtered_methylObj=filterByCoverage(methylObj,lo.count=5,lo.perc=NULL, hi.count=NULL,hi.perc=99.9) filtered_methylObj_norm <- methylKit::normalizeCoverage(filtered_methylObj) meth_filter <- methylKit::unite(filtered_methylObj_norm,destrand=FALSE, min.per.group = 4L) save(meth_filter, file = "../output/12-Peve-WGBS/methylkit/MethylObj_filtered.RData") ``` ```{r} load("../output/12-Peve-WGBS/methylkit/MethylObj_filtered.RData") PCASamples(meth_filter) clusterSamples(meth_filter, dist = "correlation", method = "ward", plot = TRUE) getCorrelation(meth_filter) ``` ## Annotation- I want an intersection of methylated CpGs with the various regions in the genome: ```{r} library("genomationData") library("genomation") library("GenomicRanges") # convert methylation data to GRange object meth_GR <- as(meth_filter, "GRanges") #import malfomed gff gff_df <- read.delim("../data/Porites_evermanni_v1.annot.gff", header = FALSE, sep = "\t", comment.char = "#", col.names = c("seqid", "source", "type", "start", "end", "score", "strand", "phase", "attributes")) %>% mutate( ID = str_extract(attributes, "(?<=ID=)[^;]+"), Parent = str_extract(attributes, "(?<=Parent=)[^;]+") ) gff_gr <- GRanges( seqnames = gff_df$seqid, ranges = IRanges(start = gff_df$start, end = gff_df$end), strand = gff_df$strand, type = gff_df$type, source = gff_df$source, score = gff_df$score, phase = gff_df$phase, attributes = gff_df$attributes,ID = gff_df$ID, Parent = gff_df$Parent ) unique(gff_df$type) # Get only transcript entries transcripts <- gff_gr[gff_gr$type == "mRNA"] ``` ### Extract methylation count matrix for transcripts ```{r} #extract percent methylation and add this to GRanges object meth_matrix <- as.data.frame(percMethylation(meth_filter)) meth_GR$meth <- meth_matrix transcript_CpG <- findOverlaps(meth_GR, transcripts) # Create a data frame with CpG-to-transcript mapping df <- data.frame( cpg_ix = queryHits(transcript_CpG), transcript_id = transcripts$ID[subjectHits(transcript_CpG)]) # Merge with CpG methylation meth_df <- as.data.frame(meth_GR)[df$cpg_ix, ] meth_df$transcript_id <- df$transcript_id meth_long <- meth_df %>% select(starts_with("meth."), transcript_id) %>% rename_with(~gsub("meth\\.", "", .x)) %>% pivot_longer( cols = -transcript_id, names_to = "sample", values_to = "perc_meth" ) %>% filter(!is.na(perc_meth)) # Summarize: Mean % methylation per transcript × sample CpG_count_transcripts <- meth_long %>% group_by(transcript_id, sample) %>% summarise(mean_meth = mean(perc_meth, na.rm = TRUE), .groups = "drop") %>% pivot_wider(names_from = sample, values_from = mean_meth) write.csv(CpG_count_transcripts, "../output/12-Peve-WGBS/CpG_Transcript_CountMat.csv") ``` ### Extract additional features ```{r} #extract exons and introns exons = gff_gr[gff_gr$type == "CDS"] introns = GenomicRanges::setdiff(transcripts, exons) exons_ranges = GenomicRanges::reduce(split(exons, exons$Parent)) # extract UTRs from gff UTR = gff_gr[gff_gr$type == "UTR"] # Reduce CDS per transcript to get full CDS range cds_ranges <- GenomicRanges::reduce(split(exons, exons$Parent)) # one merged range per Parent cds_flat <- unlist(cds_ranges) cds_df <- data.frame( Parent = names(cds_flat), cds_start = start(cds_flat), cds_end = end(cds_flat), strand = as.character(strand(cds_flat)) ) # Make UTR dataframe with index utr_df <- data.frame( idx = seq_along(UTR), start = start(UTR), end = end(UTR), Parent = UTR$Parent, strand = as.character(strand(UTR)) ) # Join UTR with CDS info utr_df <- left_join(utr_df, cds_df, by = c("Parent","strand")) # Classify UTRs utr_df$utr_type <- case_when( utr_df$strand == "+" & utr_df$end < utr_df$cds_start ~ "5'UTR", utr_df$strand == "+" & utr_df$start > utr_df$cds_end ~ "3'UTR", utr_df$strand == "-" & utr_df$start > utr_df$cds_end ~ "5'UTR", utr_df$strand == "-" & utr_df$end < utr_df$cds_start ~ "3'UTR", TRUE ~ "Internal/Overlapping" ) # Add back to GRanges UTR$utr_type <- NA UTR$utr_type[utr_df$idx] <- utr_df$utr_type gff5UTR <- UTR[UTR$utr_type == "5'UTR"] gff3UTR <- UTR[UTR$utr_type == "3'UTR"] ``` #### Transposable elements: ```{r,eval=FALSE} # Import Transposable element file from repeat masker format #rmsk <- read.table("../data/apul.hifiasm.s55_pa.p_ctg.fa.k32.w100.z1000.ntLink.5rounds.fa.out", # skip = 3, fill = TRUE, stringsAsFactors = FALSE, quote = "") %>% drop_na() # Assign column names based on RepeatMasker .out format colnames(rmsk)[c(5, 6, 7, 10)] <- c("seqname", "start", "end", "repeat_name") # Create a GRanges object TE <- GRanges( seqnames = rmsk$seqname, ranges = IRanges(start = as.numeric(rmsk$start), end = as.numeric(rmsk$end)), strand = "*", repeat_name = rmsk$repeat_name ) TE$type <- "transposable_element" head(TE) ``` #### miRNA: ```{r} miRNAs = genomation::gffToGRanges("../output/05-Peve-sRNA-ShortStack_4.1.0/ShortStack_out/known_miRNAs.gff3") head(miRNAs) ``` #### lncRNA: ```{r} lncRNAs = genomation::gffToGRanges("../output/07-Peve-lncRNA-matrix/Peve-lncRNAs.gtf") head(lncRNAs) ``` ## Plotting annotation information ```{r} # Define intergenic = genome - all annotations # First combine all annotated regions all_annotated <- GenomicRanges::reduce(c(gff5UTR, gff3UTR, exons, introns)) # Define genome bounds from methylation object genome_range <- GenomicRanges::reduce(meth_GR) # Intergenic = genome - annotated intergenic <- GenomicRanges::setdiff(genome_range, all_annotated) ``` ```{r} # Priority list region_priority <- list( `5UTR` = gff5UTR, `3UTR` = gff3UTR, exon = exons, intron = introns, intergenic = intergenic ) # Start with all CpGs unassigned cpg_annot <- rep(NA, length(meth_GR)) # Keep track of which CpGs are still unassigned unassigned <- rep(TRUE, length(meth_GR)) # Assign in priority order for (region_name in names(region_priority)) { region <- region_priority[[region_name]] transcript_CpG <- findOverlaps(meth_GR[unassigned], region) # Map the indices back to full set full_indices <- which(unassigned)[queryHits(transcript_CpG)] cpg_annot[full_indices] <- region_name unassigned[full_indices] <- FALSE # Mark these as assigned } # Replace any remaining NA with "unassigned" or "intergenic" cpg_annot[is.na(cpg_annot)] <- "intergenic" ``` ```{r} df_annotated <- as.data.frame(meth_GR) df_annotated$region <- cpg_annot #set as factor df_annotated$region <- factor(df_annotated$region, levels = c("intergenic", "3UTR", "5UTR", "intron", "exon")) blue_palette <- c("intergenic" = "#c6dbef", "3UTR" = "#9ecae1","5UTR" = "#6baed6","intron"= "#3182bd","exon" = "#08519c" ) ggplot(df_annotated, aes(x = strand, fill = region)) + geom_bar(position = "fill",color="black") + theme_minimal() + labs(y = "% CpGs", fill = "Genomic Region") + theme(axis.title.x = element_blank(),axis.text.x = element_blank(), panel.grid.minor = element_blank(),panel.grid.major.x = element_blank()) + scale_fill_manual(values = blue_palette) ``` Percent meth ```{r} # Calculate average across all samples avg_meth <- rowMeans(meth_matrix, na.rm = TRUE) # Add to annotated data df_annotated$avg_meth <- avg_meth # Classify methylation status low_thresh <- 10 high_thresh <- 50 df_annotated$meth_status <- cut(df_annotated$avg_meth, breaks = c(0, low_thresh, high_thresh, 100), labels = c("lowly", "moderately", "highly")) ``` ```{r} ggplot(df_annotated, aes(x = region, fill = meth_status)) + geom_bar(position = "fill") + # stacked bar normalized to proportions theme_minimal() + labs( title = "Proportion of CpG Methylation Status by Region", x = "Genomic Region", y = "Proportion of CpGs", fill = "Methylation Status" ) + scale_fill_manual(values = c(lowly = "#3498db", moderately = "#bdc3c7", highly = "#e74c3c")) + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + scale_y_continuous(labels = scales::percent_format()) ``` ```{r} df_summary <- df_annotated %>% group_by(region) %>% summarise(mean_meth = mean(avg_meth, na.rm = TRUE)) df_summary ``` ```{r} # Add region info meth_long <- as.data.frame(meth_matrix) meth_long$region <- cpg_annot # Reshape: One row per CpG × sample meth_long_long <- meth_long %>% pivot_longer( cols = -region, names_to = "sample", values_to = "perc_meth" ) # Now summarize: mean methylation per sample × region meth_summary <- meth_long_long %>% group_by(sample, region) %>% summarise(mean_meth = mean(perc_meth, na.rm = TRUE), .groups = "drop") ggplot(meth_summary, aes(x = region, y = mean_meth)) + geom_boxplot(outlier.shape = NA) + geom_jitter(aes(color = sample), width = 0.1) + theme_minimal() + labs( title = "Sample-Wise Average Methylation by Genomic Region", x = "Genomic Region", y = "Mean % Methylation per Sample" ) + theme(axis.text.x = element_text(angle = 45, hjust = 1)) ```