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

From ZhangLabWiki
Jump to navigation Jump to search
>Dinh
>Dinh
Line 105: Line 105:
| align="center" style="background:#f0f0f0;"|'''T98G'''
| align="center" style="background:#f0f0f0;"|'''T98G'''
|-
|-
| 0||0||8||21||15||32||18
| 0||2||19||19||24||24||21
|-
|-
| 0.1||6||50||82||86||94||92
| 0.1||8||91||103||121||109||123
|-
|-
| 0.15||10||58||91||97||107||102
| 0.15||14||108||122||135||132||146
|-
|-
| 0.2||11||66||98||107||124||110
| 0.2||38||128||143||152||145||160
|-
|-
| Total Reads Analyzed||177||177||178||176||173||174
| Total Reads Analyzed||159||174||171||174||174||175
|}
|}


Line 119: Line 119:
* 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.  


* Since 1-6% of blood reads were identified as non-blood, this means that we can only recognize cancer when cancer sample is mixed at >> 6%.
* Since up to 1% of blood reads were identified as non-blood (MaxScoreFilter=0), this means that we can only recognize cancer when cancer sample is mixed at >> 1%.


* How to reject more blood reads and improve cancer detection specificity?
* How to reject more blood reads and improve cancer detection specificity?

Revision as of 01:38, 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 reads with 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 ~200 reads from 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 ~200 of another blood sample to compare (Indx8).
  • Round 1 of <=200 random reads analyzed
MaxScoreFilter BLOOD BE2C U87MG BXPC3 PANC1 T98G
0 1 11 20 16 19 21
0.1 6 83 105 105 115 129
0.15 12 97 128 126 129 141
0.2 39 123 150 143 141 152
Total Reads Analyzed 169 174 180 176 173 178
  • Round 2 of <= 200 random reads analyzed
MaxScoreFilter BLOOD BE2C U87MG BXPC3 PANC1 T98G
0 2 19 19 24 24 21
0.1 8 91 103 121 109 123
0.15 14 108 122 135 132 146
0.2 38 128 143 152 145 160
Total Reads Analyzed 159 174 171 174 174 175
  • Cancer sample reads were identified from less than 200 total cancer reads.
  • 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.
  • Since up to 1% of blood reads were identified as non-blood (MaxScoreFilter=0), this means that we can only recognize cancer when cancer sample is mixed at >> 1%.
  • How to reject more blood reads and improve cancer detection specificity?
 (1) Use cancer haplotype information to reject non-cancer reads.
 (2) Investigate the locus on blood reads that were identified as non-blood, are they overlapping SNPs? 
     We need to remove locus that contains CpG-SNPs.
 (3) Try larger haplotype sizes with 5-7 CpGs.

Calculate the r2 between pairs of CpGs