Dinh/Dinh 2014/NOTES/2014-9-15: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Dinh
mNo edit summary
>Dinh
Line 33: Line 33:
  SRR1232303.200105      0      chr1    30557958        255    90M    *      0      0      ATTCCTTTATTTTGGTGGTTGGGGTGTTTTAGGGTTTTTTTTTTAATTAATTTTTTTTTCGTTGTTTTTTTTTTTTTGTTTTATGTTTTT      @CCFFFFFFHHHHJICGHHIJJJJ:DFHIJHJJJHJJJJJJHFDCDDEECEEEDDDDDDDDDDDDDDDDDDDDDDDD38ACD3:(:ADDD      NM:i:3  ZS:Z:++
  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:-+
  SRR1232303.200112      16      chr1    30559194        255    90M    *      0      0      ATCCCTTCCCTTTCACAACCTACTCCGAAAACCGCTAACCTAAAAATACAGCAATCTCTACTATCCTAAAAAAACAAAATCCCAAGGAAT      DDDDDDDDBDDEDDDDDDDDCACAFFFHHGIJJIJIHHJIHJJJIIGJJJJJIJGIIJJIJJIHJIHJJJJJJJJJJHHHHHFFFFFCCC      NM:i:2  ZS:Z:-+
* BSMAP have extra ZS:Z:-+ field  
* BSMAP have extra ZS:Z:-+ field that's just the strand information.
* 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.
* 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.

Revision as of 17:35, 15 September 2014

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
  • 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 map PC-T-7_2 with BSMAP and try to run MethylPurify again and 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.
  • 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.