Dinh/Dinh 2014/NOTES/2014-9-15

From ZhangLabWiki
Revision as of 17:54, 18 September 2014 by >Dinh
Jump to navigation Jump to search

MethylPurify

Test run with WGBS data and BSMAP aligner

  • I tried running MethylPurify on our RRBS MONOD bam files, but could never get the program to complete running/obtain informative bins and mixing ratio.
  • Two major issues with MethylPurify: requires enough reads to get minimum depth of coverage over informative bins and may need to realign data with BSMAP.
  • Downloaded WGBS sequencing data from the paper: SRR1232303 (BioSample: SAMN02725471; lung adenocarcinoma tissue, male), containing ~20 million PE reads x 90 bp each.
  • Downsampled to get ~25% of only reads 1 to get 5,077,371 x 90 bp single end reads to perform testing.
  • Performed mapping with BSMAP as that aligner is listed by MethylPurify to get the necessary bam input file.
  • BSMAP mapping command:
 /bsmap-2.74/bsmap -u -s 12 -v 0.04 -p 4 -a downsampled.SRR1232303_1.fastq -d hg19.fa -o downsampled.SRR1232303_1.fastq.sam
  • BSMAP took 36,221 seconds (10.06 hours) to complete mapping for ~5 million SE x 90 bp reads with 4 CPU threads.
  • Performed deconvolution with MethylPurify:
 python2.7 /home/ddiep/Downloads/MethylPurify-2.0-20140819/methylpurify/bin/MethylPurify -f downsampled.SRR1232303_1.fastq.bam 
      -g hg19.fa -i /home/ddiep/Downloads/MethylPurify-2.0/methylpurify/db/CGI_hg19_slop1000.bed -c 10 -s 50 -b 300
  • MethylPurify outputs:
 ...
 The predicted mixing ratio is:  0.265
 The number of informative bins is:  29.0
 The number of informative bins is less than 400, too less for the prediction of mixing ratio, continue
 INFO  @ Mon, 15 Sep 2014 09:33:27: Running get_methyl_profile:
 MixingRatio:0.265
  • CGI_hg19_slop1000.bed was only found in the MethylPurify-2.0.tar.gz package, but not the latest package. I couldn't find the genome.sh script that was necessary to generate this, so I need to figure out how to generate this file in the future because I am not 100% confident this is the appropriate file to use.

Testing with MONOD RRBS bam files from BisReadMapper

  • Since WGBS and RRBS may cover different CGIs, the CGI selected regions by MethylPurify have this depth of coverage:
 samtools depth XXX.CGI.sorted.bam  |  awk '{sum+=$3} END { print "Sum = ", sum, " Average = ",sum/NR}' 
 >> Sum =  41941637  Average =  2.30644
  • Our RRBS data have even better coverage at CGIs, the CGI regions selected by MethylPurify for PC-T-7_2 have this depth of coverage:
 samtools depth PC-T-7_2.merged.CGI.sorted.bam  |  awk '{sum+=$3} END { print "Sum = ", sum, " Average = ",sum/NR}'
 >> Sum =  393749988  Average =  27.9027
  • This could means that MethylPurify failed when I ran on our RRBS bam files because it couldn't understand the alignments from bisReadMapper, but requires the BSMAP alignment format (may have special instructions for MethylPurify). I will figure out the differences in BSMAP and BisReadMapper bam files.
  • BSMAP sam file:
SRR1232303.200105       0       chr1    30557958        255     90M     *       0       0       ATTCCTTTATTTTGGTGGTTGGGGTGTTTTAGGGTTTTTTTTTTAATTAATTTTTTTTTCGTTGTTTTTTTTTTTTTGTTTTATGTTTTT      @CCFFFFFFHHHHJICGHHIJJJJ:DFHIJHJJJHJJJJJJHFDCDDEECEEEDDDDDDDDDDDDDDDDDDDDDDDD38ACD3:(:ADDD      NM:i:3  ZS:Z:++
SRR1232303.200112       16      chr1    30559194        255     90M     *       0       0       ATCCCTTCCCTTTCACAACCTACTCCGAAAACCGCTAACCTAAAAATACAGCAATCTCTACTATCCTAAAAAAACAAAATCCCAAGGAAT      DDDDDDDDBDDEDDDDDDDDCACAFFFHHGIJJIJIHHJIHJJJIIGJJJJJIJGIIJJIJJIHJIHJJJJJJJJJJHHHHHFFFFFCCC      NM:i:2  ZS:Z:-+
  • BSMAP have extra ZS:Z field that's just the strand information. MethylPurify would fail similarly as for BisReadMapper bam files, after I removed the ZS:Z field from the BSMAP bam file. I may just need to generate the ZS:Z field to get MethylPurify to work on BisReadMapper bam files.
  • Encoding for ZS:Z field is as follows:
       ++: forward strand of Watson of reference (BSW)
       +-: reverse strand of Watson of reference (BSWC)
       -+: forward strand of Crick of reference (BSC)  
       --: reverse strand of Crick of reference (BSCC)
  • Make MethylPurify think that BisReadMapper bam files are from BSMAP: Add in NM:i:0 and ZS:Z:-+ (for Crick) and ZS:Z:++ (for Watson).
 samtools view $f | awk '{if($2 == 0) print $0"\tNM:i:0\tZS:Z:++"; if($2 == 16) print $0"\tNM:i:0\tZS:Z:-+"}' | samtools view -bSt /media/2TB_storeA/BisRef/Hg19_plusLamba/hg19_lambda.fa.fai - > mod.$f
  • Run MethylPurify:
 python2.7 /home/ddiep/Downloads/MethylPurify-2.0-20140819/methylpurify/bin/MethylPurify -f ../mod.PC-T-7_2.merged.bam -g hg19.fa -i /home/ddiep/Downloads/MethylPurify-2.0/methylpurify/db/CGI_hg19_slop1000.bed -c 10 -s 50 -b 300
  • MethylPurify outputs:
 ...
 The predicted mixing ratio is:  0.175
 The number of informative bins is:  304.0
 The number of informative bins is less than 400, too less for the prediction of mixing ratio, continue
 INFO  @ Mon, 15 Sep 2014 13:09:12: Running get_methyl_profile:
 MixingRatio:0.175
  • Results from MethylPurify:
 File:MixRatio.MethylPurify test1.png
 File:InformativeBins.MethylPurify test1.png

MethylPurify run on merged duplicates tumor samples

  • Results from MethylPurify:
 File:MixRatio.MethylPurify testMergedTumors.png
 File:InformativeBins.MethylPurify testMergedTumors.png
  • I merged all of the informative bins for colon cancer, lung cancer, and pancreatic cancers together into one BED file.
 File:Colon cancer methylpurify informativeCGIs merged.txt
 File:Lung cancer methylpurify informativeCGIs merged.txt
 File:Pancreatic cancer methylpurify informativeCGIs merged.txt
  • I uploaded the BED files to GREAT & used default settings to get the biological processes that are enriched. This may be because the regions are all CGIs, so not certain whether the signal is real or not.
  • Results from MethylPurify:
 File:Colon cancer informative CGIs biologicalprocess.png
 File:Lung cancer informative CGIs biologicalprocess.png
 File:Pancreatic cancer informative CGIs biologicalprocess.png