Dinh/Dinh 2014/NOTES/2014-1-23: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Dinh
(Created page with "== Haplotype Analysis == === SAM record parser === * In order to perform haplotype analysis, I needed to write a SAM record parser so that I can analyze each read and determin...")
 
>Dinh
Line 20: Line 20:
  sam parser:  
  sam parser:  
  chr9 100610863 100610864 '9/49' 49 + 100610863 100610864 30,210,0
  chr9 100610863 100610864 '9/49' 49 + 100610863 100610864 30,210,0
* Fixed an error with extractMethyl.pl, the Regex for indels can't be used to remove them.
  Regex for indels: \+[0-9]+[ACGTNacgtn]+' \-[0-9]+[ACGTNacgtn]+'
  Example:
  ...TTTT.T.T..TTT-1GTT-1GT-1GTTTTT-1GTTTTT-1GT-1GTT-1GTTTTTTTTTTTTTT..T
                  ---------------------------------------------------
                  ***  *** ***    ***    *** ***  ***
  The positions marked with dashes at the bottom were removed using the regex for indels.
  However, only the positions marked with stars should be removed.
* Also, samtools pileup always convert 64 to 33 (even when the base quality is already in 33).
* Re-generate BED files from extractMethyl & calculated the correlation between the two
=== Count locus alleles ===
=== Count locus alleles ===
* Count the different alleles produced by four consecutive CpGs on the same read. There are 16 different possibilities.
* Count the different alleles produced by four consecutive CpGs on the same read. There are 16 different possibilities.

Revision as of 20:15, 24 January 2014

Haplotype Analysis

SAM record parser

  • In order to perform haplotype analysis, I needed to write a SAM record parser so that I can analyze each read and determine the CG positions and their methylation status. This works like samtools pileup but since the current version of samtools mpileup requires a bam file as input and a reference fasta file as input, it is probably more efficient to skip using samtools mpileup.
  • The SAM parser needs to identify the position of each read base, and the mismatch positions. Since there are INDELS and H/S clippings, it needs to integrate both the CIGAR and MD:Z fields.
  • Script for finding mismatch positions in SAM file: File:DdSAMparse.txt
  • Script for counting the CpG locus alleles (also makes BED): File:CgLocusAlleles.txt
  • To compare the parsing results, use sam parser to generate a BED format file for CpG positions & compare with results obtained from samtools pileup -> methylFreq -> BED
NP-GP1V6-Ind1-16-Dec26-Ind10_S10.BED.txt : cor = 0.9999383
NP-GP1V6-Ind1-16-Dec26-Ind11_S11.BED.txt : cor = 0.9999235
NP-GP1V6-Ind1-16-Dec26-Ind12_S12.BED.txt : cor = 0.9999451
NP-GP1V6-Ind1-16-Dec26-Ind13_S13.BED.txt : cor = 0.9999285
NP-GP1V6-Ind1-16-Dec26-Ind14_S14.BED.txt : cor = 0.9998338
NP-GP1V6-Ind1-16-Dec26-Ind15_S15.BED.txt : cor = 0.9999977
NP-GP1V6-Ind1-16-Dec26-Ind16_S16.BED.txt : cor = 0.9999913
NP-GP1V6-Ind1-16-Dec26-Ind9_S9.BED.txt   : cor = 0.9999657
  • There are some positions where in the SAM record parsing produced slightly different results from samtools pileup.
An example of non-concordance call:
samtools pileup: 
chr9	100610863	100610864	'9/19'	19	+	100610863	100610864	90,150,0
sam parser: 
chr9	100610863	100610864	'9/49'	49	+	100610863	100610864	30,210,0
  • Fixed an error with extractMethyl.pl, the Regex for indels can't be used to remove them.
 Regex for indels: \+[0-9]+[ACGTNacgtn]+' \-[0-9]+[ACGTNacgtn]+' 
 Example:
 ...TTTT.T.T..TTT-1GTT-1GT-1GTTTTT-1GTTTTT-1GT-1GTT-1GTTTTTTTTTTTTTT..T
                 ---------------------------------------------------
                 ***  *** ***     ***     *** ***  *** 
 The positions marked with dashes at the bottom were removed using the regex for indels.
 However, only the positions marked with stars should be removed.
  • Also, samtools pileup always convert 64 to 33 (even when the base quality is already in 33).
  • Re-generate BED files from extractMethyl & calculated the correlation between the two

Count locus alleles

  • Count the different alleles produced by four consecutive CpGs on the same read. There are 16 different possibilities.
  • Examples of the locus count:
chr1:1475064:1475086:1475088:1475095    2UMMM,1MMUM,28MMMM,     31
chr1:1475086:1475088:1475095:1475100    1MUMM,30MMMM,   31
chr1:1475095:1475100:1475108:1475130    32MMMM, 32
chr1:1475108:1475130:1475141:1475143    1MMUM,1MMMU,30MMMM,     32
chr1:1475130:1475141:1475143:1475162    1MMUM,1MMMU,29MMMM,     31
chr1:1475141:1475143:1475162:1475164    1MUMM,1MMUM,29MMMM,     31

Calculate the r2 between pairs of CpGs