Dinh/Dinh 2014/NOTES/2014-9-15: Difference between revisions
Jump to navigation
Jump to search
>Dinh mNo edit summary |
>Dinh mNo edit summary |
||
Line 23: | Line 23: | ||
MixingRatio:0.265 | 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: | * 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}' | samtools depth XXX.CGI.sorted.bam | awk '{sum+=$3} END { print "Sum = ", sum, " Average = ",sum/NR}' | ||
Line 41: | Line 43: | ||
* 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). | * 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 | 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 |
Revision as of 20:13, 15 September 2014
MethylPurify
- Tumor deconvolution paper: http://genomebiology.com/2014/15/8/419
- Python package: https://pypi.python.org/pypi/MethylPurify/2.0-20140819
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