Dinh/Dinh 2014/NOTES/2014-9-15
Jump to navigation
Jump to search
MethylPurify
- Tumor deconvolution paper: http://genomebiology.com/2014/15/8/419
- Python package: https://pypi.python.org/pypi/MethylPurify/2.0-20140819
- Installed MethylPurify on genome-miner. The scripts dont take much memory & finish relatively quickly.
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 GREAT:
File:Colon cancer informative CGIs biologicalprocess.png File:Lung cancer informative CGIs biologicalprocess.png File:Pancreatic cancer informative CGIs biologicalprocess.png
Overlap of informative CGIs with blood UMR/LMRs
- Blood UMRs and LMRs were identified using MethylSeekR
- I used bedtools intersect and required that 50% of the informative bins be overlapped by UMRs/LMRs
- Results:
' | colon | lung | pancreatic |
Number of informative CGIs | 2852 | 1248 | 1798 |
overlap with UMR | 2274 | 1144 | 1699 |
overlap with LMR | 12 | 14 | 14 |
%overlap with UMR | 79.73% | 91.67% | 94.49% |
%overlap with LMR | 0.42% | 1.12% | 0.78% |
% no overlap | 19.85% | 7.21% | 4.73% |
- Majority of informative CGIs are in UMRs.
- Since only a fraction of informative CGIs are overlapping with LMRs, I looked for region/functional enrichment using GREAT with the CGIs overlapping LMRs as test data and all the informative CGIs as background. No term was enriched.
- Since most of the informative CGIs are overlapping with UMRs, no term was enriched if I used all informative CGIs as background. So I used the whole genome was background. The GREAT results show developmental associated terms enriched in biological processes.
- GREAT results:
File:Colon cancer informative CGIs UMRs biologicalprocess.png File:Lung cancer informative CGIs UMRs biologicalprocess.png File:Pancreatic cancer informative CGIs UMRs biologicalprocess.png
Discussion
- Discussion with Dr. Zhang (9/12/14):
- Intersect informative bins with blood UMRs, what are those regions enriched for in GREAT? (DONE)
- How to tell whether the minor or major component is the tumor DNA ? It appears that these tumor references have varying amount of minor component.
- Ave. methyl colon cancer tumors: major = 0.7348, minor = 0.1602
- Ave. methyl lung cancer tumors: major = 0.6393, minor = 0.2574
- Ave. methyl pancreatic cancer tumors: major = 0.7623, minor = 0.1181
- Although majority of informative CGIs are overlapping with blood UMRs, the minor components are not always unmethylated, thus there may be other normal cells in the mixture.
- Use MethylPurify to generate cancer reference from primary tumors. So we can create a database of true cancer references (after substracting out normal cells). We can use the pure cancer reference for methylation haplotyping to detect cancer reads in blood.