Dinh/Dinh 2014/NOTES/2014-1-23: Difference between revisions
Jump to navigation
Jump to search
>Dinh |
>Dinh |
||
Line 94: | Line 94: | ||
|} | |} | ||
* 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 | * 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. | ||
* | * Next, I decreased the number of reads to <=200 | ||
{| class="wikitable" cellpadding=2 | |||
| align="center" style="background:#f0f0f0;"|'''MaxScoreFilter''' | |||
| align="center" style="background:#f0f0f0;"|'''BLOOD''' | |||
| align="center" style="background:#f0f0f0;"|'''BE2C''' | |||
| align="center" style="background:#f0f0f0;"|'''U87MG''' | |||
| align="center" style="background:#f0f0f0;"|'''BXPC3''' | |||
| align="center" style="background:#f0f0f0;"|'''PANC1''' | |||
| align="center" style="background:#f0f0f0;"|'''T98G''' | |||
|- | |||
| 0||3||12||12||15||14||15 | |||
|- | |||
| 0.1||7||41||66||79||96||74 | |||
|- | |||
| 0.15||7||47||73||88||106||88 | |||
|- | |||
| 0.2||7||57||88||106||113||99 | |||
|- | |||
| Total Reads Analyzed||172||180||172||175||168||175 | |||
|} | |||
* Cancer sample reads were identified from less than 200 total cancer reads. | |||
* Since 1-4% of blood reads were identified as non-blood, we need the level of cancer to be greater than >>1.1% to make a cancer call. We can also score the identified non-blood read by comparing them to cancer haplotypes to increase confidence. | |||
=== Calculate the r2 between pairs of CpGs === | === Calculate the r2 between pairs of CpGs === | ||
* | * |
Revision as of 00:11, 25 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.
- Next, I decreased the number of reads to <=200
MaxScoreFilter | BLOOD | BE2C | U87MG | BXPC3 | PANC1 | T98G |
0 | 3 | 12 | 12 | 15 | 14 | 15 |
0.1 | 7 | 41 | 66 | 79 | 96 | 74 |
0.15 | 7 | 47 | 73 | 88 | 106 | 88 |
0.2 | 7 | 57 | 88 | 106 | 113 | 99 |
Total Reads Analyzed | 172 | 180 | 172 | 175 | 168 | 175 |
- Cancer sample reads were identified from less than 200 total cancer reads.
- Since 1-4% of blood reads were identified as non-blood, we need the level of cancer to be greater than >>1.1% to make a cancer call. We can also score the identified non-blood read by comparing them to cancer haplotypes to increase confidence.