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

From ZhangLabWiki
Jump to navigation Jump to search
>Dinh
>Dinh
mNo edit summary
Line 41: Line 41:
  chr1:1475141:1475143:1475162:1475164    1MUMM,1MMUM,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.
* 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.
=== Pearson's Chi-square test ===
=== Pearson's Chi-square test ===
* We need to filter for locus where blood is significantly different from cancer but not signficantly different from blood.
* We need to filter for locus where blood is significantly different from cancer but not signficantly different from blood.
Line 89: Line 90:
=== Identify cancer reads by comparing against blood haplotype profile ===
=== Identify cancer reads by comparing against blood haplotype profile ===
* For any read, calculate the probability of seeing a particular allele at a locus in blood sample using the haplotype information for blood:
* For any read, calculate the probability of seeing a particular allele at a locus in blood sample using the haplotype information for blood:
* Total probability is obtained by multiplying all of the found probabilities.
* Example 1:
* Example 1:
   M01186:73:000000000-A20K5:1:1101:15873:5723:1:N:0:3_AGCCTT:R
   M01186:73:000000000-A20K5:1:1101:15873:5723:1:N:0:3_AGCCTT:R
Line 100: Line 102:
   chr4:172734442:172734445:172734453:172734455  UMMM  0
   chr4:172734442:172734445:172734453:172734455  UMMM  0
   chr4:172734445:172734453:172734455:172734459 MMMM  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).
   Total 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:
* Example 2:
Line 115: Line 117:
   chr2:26407725:26407728:26407732:26407741, UUUU 0.966
   chr2:26407725:26407728:26407732:26407741, UUUU 0.966
   chr2:26407728:26407732:26407741:26407743, UUUU 0.968
   chr2:26407728:26407732:26407741:26407743, UUUU 0.968
   The average probability is high, meaning that this read might have also originated from a blood sample. (This cannot be distinguished from blood).
   The total probability is high, meaning that this 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.
* 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.

Revision as of 21:59, 27 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.

Pearson's Chi-square test

  • We need to filter for locus where blood is significantly different from cancer but not signficantly different from blood.
  • Use the Chi-Square test of independence (Pearson's Chi square) to compare pairwise between blood and each cancer type. Output only the locus count for blood where blood counts is independent from the cancer counts.
  • Perl script : File:ChisqIndep cgLocusAlleles.txt
 chiSqIndep_cgLocusAlleles [cgLocusAlleles file A] [cgLocusAlleles file b]
  • Example ChiSqIndep for Locus: chr6:392409:392417:392428:392435 (between Blood_B1 and PANC1)
 Rows: 2
 Columns: 9
 Degree of Freedom: 8
 Total Count: 136
 Chi-square Statistic: 136
 p-value: 0.00000000000000000000000016093
 Warning: some of the cell counts might be too low. 
       UUUU 	UUUM	UUMM	UMMM	MUUU	MUUM	MUMM	MMUM	MMMM	rtotal
 Blood	30	0	0	0	0	0	0	0	0	30
       (6.62)	(0.88)	(0.66)	(0.22)	(0.44)	(2.65)	(2.87)	(0.66)	(15.00)	
 PANC1	0	4	3	1	2	12	13	3	68	106
       (23.38)	(3.12)	(2.34)	(0.78)	(1.56)	(9.35)	(10.13)	(2.34)	(53.00)	
ctotal	30	4	3	1	2	12	13	3	68	136
  • Example ChiSqIndep for Locus: chr5:1887153:1887159:1887179:1887185 (between Blood_B1 and Blood_H11)
 Rows: 2
 Columns: 10
 Degree of Freedom: 9
 Total Count: 437
 Chi-square Statistic: 24.8046419350707
 p-value: 0.0031951
 Warning: some of the cell counts might be too low. 
       UUUU	UUUM	UUMU	UUMM	UMUU	UMUM	UMMM	MUUU	MUMM	MMUM	rtotal
 B1	350	2	2	1	2	0	0	2	1	2	362
      (345.43)	(4.14)	(2.49)	(1.66)	(1.66)	(0.83)	(1.66)	(1.66)	(0.83)	(1.66)	
 H11	67	3	1	1	0	1	2	0	0	0	75
      (71.57)	(0.86)	(0.51)	(0.34)	(0.34)	(0.17)	(0.34)	(0.34)	(0.17)	(0.34)	
 ctotal	417	5	3	2	2	1	2	2	1	2	437
  • Summary count of independent locus, using Blood_B1 as representative blood profile against Blood_H11 and Blood_D1
  2162 Blood_B1.1.cgLocusAlleles.indepVersus.Blood_D1.1.cgLocusAlleles
  1435 Blood_B1.1.cgLocusAlleles.indepVersus.Blood_H11.1.cgLocusAlleles
  • Summary count of independent locus, Blood_B1 versus cancer haplotypes
  8724 Blood_B1.1.cgLocusAlleles.indepVersus.BE2C_1.1.cgLocusAlleles
 11462 Blood_B1.1.cgLocusAlleles.indepVersus.BXPC3_1.1.cgLocusAlleles
 10516 Blood_B1.1.cgLocusAlleles.indepVersus.PANC1_1.1.cgLocusAlleles
 11428 Blood_B1.1.cgLocusAlleles.indepVersus.T98G_1.1.cgLocusAlleles
 11020 Blood_B1.1.cgLocusAlleles.indepVersus.U87MG_1.1.cgLocusAlleles
  • Generate a final file combining all unique independent locus of blood versus cancer, then remove any locus where Blood_B1 is independent from other blood profiles.
 10738 Blood_Indep_All.cgLocusAlleles

Identify cancer reads by comparing against blood haplotype profile

  • For any read, calculate the probability of seeing a particular allele at a locus in blood sample using the haplotype information for blood:
  • Total probability is obtained by multiplying all of the found probabilities.
  • 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
 Total 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 total probability is high, meaning that this 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 (Blood_B1, Indx6). As a control, I sampled ~200 of another blood sample to compare (Blood_H11, Indx8).
  • Round 1 of <=200 random reads analyzed
MaxScoreFilter BLOOD BE2C U87MG BXPC3 PANC1 T98G
0 1 11 22 21 24 13
0.1 4 58 85 96 97 82
0.15 7 74 97 103 109 90
0.2 24 92 108 112 116 98
Total Reads Analyzed 117 135 132 129 132 118
  • 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