Dinh/Dinh 2014/NOTES/2014-1-23
Jump to navigation
Jump to search
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
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