Comparing Crassostrea spp. methylation patterns

Guess I got to put together a poster for ASLO! My goal is to showcase some preliminary comparisons of methylation patterns in C. gigas and C. virginica gonad tissue, and DML in reponse to experimental ocean acidification.

Method comparison

Although the experiments are similar, there are some key differences in the experimental design that I need to consider. I made some summary tables to easily compare the experiments.

The first takeaway is that the C. gigas experimental duration was longer (7 weeks vs. 4 weeks), and potentially more extreme. Although the treatment pH conditions were the same for both species, pCO₂ was higher for C. gigas in both experimental conditions.

Table 1. Carbonate chemistry parameters for both experiments. The C. gigas experiment took place over 7 weeks and had ambient inflow, while the C. virginica experiment was 4 weeks long and had controlled, not ambient, inflow.

The C. gigas analysis had a smaller sample size. Pooled samples were used since I had lower DNA yield from histology blocks than from the tissue Katie and Alan sent. To compensate for the smaller sample size, I used more stringent analysis parameters (minimum sequencing depth used in downstream analyses and bismark alignment score). Another important note is that I only sequenced stage 2 female gonad for C. gigas, but sex and maturity are unknown for C. virginica.

Table 2. Sequencing and analysis parameters. For C. gigas, one pooled sample with two individuals of the same sex and reproductive sstage was created for each treatment. Alignment score refers to the score_min parameter used by bismark.

Methylation island analysis

I want to replicate this perl script from this paper to identify methylation islands (MI) in the C. gigas and C. virginica genomes. The goal is to create a MI track for each species that I can then compare with DML in response to ocean acidification.

My perl comprehension is pretty low…but from what I can tell, all I need to do is:

• Choose a window size: Probably going to start with 200 bp since that’s what the paper I’m replicating used, but I’ll try 300 bp too
• Define mCpG fraction: Again, starting with 0.02 since that’s what the authors used, but I want to play around with this since insect genomes are less methylated than C. virginica.
• Choose a step size: 50 bp (again from the paper)
• Provide a list of mCpG in the genome (column 1: scaffold, column 2: bp): The paper used Bis-class to identify mCpG, but I already identified CpG as methylated using a 50% methylation cut-off for C. virginica, and Claire and Mac have done this previously with C. gigas

I’m going to start with C. virginica. Since my BEDfile of methylated loci has three columns instead of two, I started by creating a new file with only two columns: chromosome and mCpG position in this Jupyter notebook.

awk '{print $1"\t"$2}' 2019-04-09-All-5x-CpG-Loci-Methylated.bed > 2019-04-09-All-5x-CpG-Loci-Methylated-Reduced.bed #Create new tab-delimited file with only chromosome and mCpG position

Then, I cloned the script to my repo and saved it here. Finally, I ran the script in my Jupyter notebook based on the README.md instructions! I ran several iterations of the script with different mCpG fractions.

#Run script with 200 bp window, 0.02 mCpG fraction, and 50 bp step size
./methyl_island_sliding_window.pl 200 0.02 50 2019-04-09-All-5x-CpG-Loci-Methylated-Reduced.bed > 2020-02-06-Methylation-Islands-200_0.02_50.tab

I tried this with various mCpG fractions and initial widow sizes and saved the output in this folder. I really should have found a way to generate all these files programmatically…but here we are!

Table 3. MI statistics for various parameters. In general, MI number decreases as mCpG fraction and window size increase.

Now that I have all of these different MI tracks, I need to visualize them in IGV to 1) make sure methylation islands make sense in general and 2) see which set of parameters best fits my data. I created a new IGV session and added in bedgraphs with the locations of all CpGs from the concatenated 5x data, as well as just the methylated CpGs, using gannet links so someone else could view my session easily. My tab-delimited files have an extra column that counts the number of mCpG in each MI. Unfortunately IGV wasn’t too happy about this extra column when I tried to import and view it. I used the following loop to create BEDfiles from my tab-delimited output:

%%bash
for f in *.tab
do
    awk '{print $1"\t"$2"\t"$3}' ${f} > ${f}.bed
done

I loaded all my BEDfiles in IGV! And looked at the MI tracks at various resolutions.

Figures 1-5. MI tracks generated using various parameters. Blue MI tracks are those from 200 bp windows, and purple MI tracks are those from 300 bp windows. The mCpG fraction increases down the screen.

After creating BEDfiles and putting them on gannet, I added them to IGV. Based on my IGV session, it seems like the 0.02 mCpG fraction actually does a pretty good job of capturing the methylation islands. I just don’t know if I should choose the 200 or 300 bp windows. I’m leaning towards 200 bp just because that’s the original method and I don’t really have a better way to choose. I posted this issue to get feedback.

Going forward

1. Generate C. gigas methylation islands
2. Charaterize C. gigas and C. virginica DML in relation to methylation islands
3. Compare general methylation landscapes
4. Complete a GO-MWU and DMG analysis with C. gigas
5. Compare C. gigas and C. virginica DML
6. Draft poster for ASLO and get feedback
7. Finalize ASLO poster and print!