Kun:LabNotes/CpgSeq/2008-5-23

From ZhangLabWiki
Jump to navigation Jump to search

Informatics

Read mapping

Should generate three files:

  1. Haplotype file: xxx.methylHap.txt
    1. One molecule per line;
    2. target_id, offset1, methylotype1, offset2, methylotype2...
  2. Methylation level file, average methylation level per site: xxx.methylFreq.txt
  3. Methylation LD statistics file, all LD statistics: xxx.methylLD.txt

Mock read generation

  • 50% methylation at every CpG site;
  • No correlation between adjacent CpG;
  • Read position is evenly distributed;
  • Can generate both single reads and pair-end reads;
  • Incorporate an error model for sequencing;
  • Should generate one raw read file and one methylHap file.
  Perl Script: Media:mockReadsGenV1.txt
  • I generated mock reads at 50X coverage (3,414,819 reads). 98% were mapped back to the target sequences with SOAP with the following parameters:
   -v 5 -e 10000 -c 0 -f 8 -r 0 -s 9
  • To find out whether the 2% unmappable reads were due to multiple mapping positions or too many C:T mismatches, I generated another set of mock reads from completely methylated sequences (so that there is no C:T mismatch). It turned out that 100% reads were successfully mapped to the targets, which means the mapping strategy is bias towards methylated sequence.
  • I generated a third set of mock reads from completely unmethylated sequences, and performed the mapping with the same parameters. Only 85% reads were mappable. In contrast, if I replace all CGs with TGs in the target sequences, 100% reads were mappable.
  • I think the best mapping strategy is to map the reads to the completely methylated targets, and to the completely unmethylated targets, then take the union of the mapping results.