Dinh/Dinh 2014/NOTES/2014-1-23: Difference between revisions
Jump to navigation
Jump to search
>Dinh |
>Dinh |
||
Line 1: | Line 1: | ||
== Haplotype Analysis == | == GP1 Haplotype Analysis == | ||
=== SAM record parser === | === 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. | * 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. |
Revision as of 23:54, 24 January 2014
GP1 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 again. They are improved slightly and even have perfect correlation for some data files.
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
- Generated a locus count table for each sample. The locus count table contains the haplotype information for each sample. Consider only locus that are covered by at least 10 reads.
- For any read, calculate the probability of seeing a particular allele at a locus in blood sample using the haplotype information for blood:
- Example 1:
M01186:73:000000000-A20K5:1:1101:15873:5723:1:N:0:3_AGCCTT:R chr4:172734347:172734362:172734393:172734395 UUMM 0 <--- This means that no blood sample read have this haplotype chr4:172734362:172734393:172734395:172734426 UMMM 0 chr4:172734393:172734395:172734426:172734431 MMMU 0 chr4:172734395:172734426:172734431:172734435 MMUM 0 chr4:172734426:172734431:172734435:172734442 MUMU 0 chr4:172734431:172734435:172734442:172734445 UMUM 0 chr4:172734435:172734442:172734445:172734453 MUMM 0 chr4:172734442:172734445:172734453:172734455 UMMM 0 chr4:172734445:172734453:172734455:172734459 MMMM 0 Average probability score is 0, meaning that this read have no chance of being originated from a blood sample. (This is a non-blood read).
- Example 2:
M01186:73:000000000-A20K5:1:2102:6845:17522:1:N:0:3_CCTTAG:F chr2:26407679:26407682:26407684:26407689, UUUU 0.966 chr2:26407682:26407684:26407689:26407691, UUUU 0.967 chr2:26407684:26407689:26407691:26407697, UUUU 0.975 chr2:26407689:26407691:26407697:26407713, UUUU 0.975 chr2:26407691:26407697:26407713:26407721, UUUU 0.979 chr2:26407697:26407713:26407721:26407723, UUUU 0.982 chr2:26407713:26407721:26407723:26407725, UUUU 0.980 chr2:26407721:26407723:26407725:26407728, UUUU 0.979 chr2:26407723:26407725:26407728:26407732, UUUU 0.971 chr2:26407725:26407728:26407732:26407741, UUUU 0.966 chr2:26407728:26407732:26407741:26407743, UUUU 0.968 The average probability is high, meaning that his read might have also originated from a blood sample. (This cannot be distinguished from blood).
- It is almost certain that 0 average scores are not from blood DNA, but it is harder to determine whether a 0.2 or even a 0.9 average score are not from blood. We may improve the confidence by increasing the haplotype size to include more than 4 consecutive CpG sites.
- To test specificity, I randomly sampled 1/1000 of each cancer BAM files (for read 1 of V4 capture ONLY), and calculate the probability score of finding each read in a blood sample (Indx6). As a control, I sampled 1/1000 of another blood sample to compare (Indx8).
MaxScoreFilter | BLOOD | BE2C | U87MG | BXPC3 | PANC1 | T98G |
0 | 0 | 45 | 44 | 108 | 52 | 71 |
0.1 | 2 | 171 | 201 | 457 | 202 | 200 |
0.15 | 3 | 194 | 225 | 515 | 215 | 221 |
0.2 | 3 | 224 | 246 | 552 | 243 | 254 |
Total Reads Analyzed | 267 | 576 | 475 | 899 | 365 | 402 |
- This analysis promisingly shows that we can positively identify many of the cancer sample reads in blood and reject most reads from blood samples in blood. It is also able to do this using fewer than 1000 cancer sample reads.
- Further check the performance using read 2 from V4, and using both reads from V6.
- Score reads based on their similarities to a cancer sample.