/
Update - 4-14-15

Update - 4-14-15

Apr 14, 2015

Update - 4/14/15

  • First of all, I sorted all my data as needed with my shell script and samtools. That all worked nicely.
  • I spent most of the week struggling with methylKit, but have been successful running it on some test datasets (small samplings of my real data). It took me a while because I don't know much about R and the documentation for methylKit assumes a high degree of familiarity.
  • Reading Bismark .sam files into Methylkit requires them to be sorted by chromosome and position, as I said last week. You then need to run the following command: 
    read.bismark(location, sample.id, assembly,
         save.folder = NULL, save.context = c("CpG"),
         read.context = "CpG", nolap = FALSE, mincov = 10,
         minqual = 20, phred64 = FALSE, treatment)
    Some notes on this:
  1. The read.context argument will take CpG, CHG, and CHH.
  2. save.folder will save the methylation percentage calls in a folder of your choosing
  3. sample.id is a string that acts as an identifier for downstream analyses (differential methylation, etc.)
  4. assembly is a string that should reference the genome assembly (TAIR10), but I don't think it has any effect on the functioning of the program
  5. Some of those arguments are optional, but the documentation doesn't offer any explanation of which.
  6. I don't know enough about sequencing at the moment to know good values for some of them, but I'm working on that.
  • Here's an example of reading a Bismark .sam file:

    When running this for real it would be much more efficient to script it. This is just for documentation purposes.
  • Once .sam files have been read in they will need to be merged into one "methbase object": 
    methyl.obj=unite(sample_list, destrand=FALSE)
    - sample_list is a list of your methylation calls (different genetic backgrounds, same context)
  • After merging into one object the following command will identify differential (hyper) methylation: 
    > myDiff = calculateDiffMeth(methyl.obj, num.cores = 8)
> myDiff25p.hyper = get.methylDiff(myDiff, difference = 25,qvalue = 0.01, type = "hyper")
  • My goal of getting all this done, plus writing the python script was a bit too ambitious last week. So, for this coming week I am going to focus solely on python now that I have the basics down as far as methylKit. Viewing my methylation calls in CoGe is a good deliverable for the end of the project.

To-do:

  1. Bismark CoGe import script. Due 4/21/15