Kun:LabNotes/MONOD/2014-8-12: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
Line 217: Line 217:
**The average read depth is lower in RRBS data. For CpG sites that have low or no coverage in normal plasma, I need to use the methylation level of whole blood WGBS data to estimate the haplotype frequency.
**The average read depth is lower in RRBS data. For CpG sites that have low or no coverage in normal plasma, I need to use the methylation level of whole blood WGBS data to estimate the haplotype frequency.
**I decided to switch to LOD (log of odds ratio) between normal plasma and tumor. Essentially, I created two reference hapInfo file, one for all normal plasma (NC-P) samples combined, another for all primary tumor samples of one cancer type (PC, CC or LC) combined. For each sequence reads in the patient plasma, I asked what's the probability that the haplotype came from blood and the probability it came from tumor. Then I took the log of the odds ratio between these two probabilities, and sum them up for each target region and for the entire genome.
**I decided to switch to LOD (log of odds ratio) between normal plasma and tumor. Essentially, I created two reference hapInfo file, one for all normal plasma (NC-P) samples combined, another for all primary tumor samples of one cancer type (PC, CC or LC) combined. For each sequence reads in the patient plasma, I asked what's the probability that the haplotype came from blood and the probability it came from tumor. Then I took the log of the odds ratio between these two probabilities, and sum them up for each target region and for the entire genome.
*I started with all RRBS targets within the UMRs, but then need to trim the target list, based on which target is informative.
*I started with all RRBS targets within the UMRs, but then need to trim the target list, based on which target is informative. Since the sample size is limited, doing cross validation is not realistic. So I simply pooled all patient plasma data (treating PC, CC, LC separately), and ran mixMethHapAnalysis_19Aug14.pl against pooled normal plasma and pooled tumor to identify the regions that have the highest LOD scores.
*To evaluate the level of false positives due to variability in normal plasma, I took one normal plasma as a test sample to compare with the pool of the other seven normal plasma. I did this for each of the 8 plasma samples, and reported a list of regions that show high LOD (false positives). These regions were excluded from the list of informative markers.

Revision as of 22:54, 24 August 2014

RRBS data analysis (continued)

1. Compile a list of RRBS targets

  • Take all the RRBS data that we generated from primary tumor samples, concatenate all the methylFreq files, and generate a single BED file.
 cat 6-T-1_1.methylFreq 6-T-1_2.methylFreq 6-T-2_1.methylFreq 6-T-2_2.methylFreq 6-T-3_1.methylFreq 6-T-3_2.methylFreq 6-T-4_1.methylFreq 6-T-4_2.methylFreq 6-T-5_1.methylFreq 6-T-5_2.methylFreq 7-T-1_1.methylFreq 7-T-1_2.methylFreq 7-T-2_1.methylFreq 7-T-2_2.methylFreq 7-T-3_1.methylFreq 7-T-3_2.methylFreq 7-T-4_1.methylFreq 7-T-4_2.methylFreq 7-T-5_1.methylFreq 7-T-5_2.methylFreq CTT-FFPE-100ng_1.methylFreq CTT-FFPE-100ng_2.methylFreq CTT-FFPE-5ng_1.methylFreq CTT-FFPE-5ng_2.methylFreq CTT-frozen-100ng_1.methylFreq CTT-frozen-100ng_2.methylFreq CTT-frozen-5ng_1.methylFreq CTT-frozen-5ng_2.methylFreq PC-T-1_1.methylFreq PC-T-1_2.methylFreq PC-T-2_1.methylFreq PC-T-2_2.methylFreq PC-T-4_1.methylFreq PC-T-4_2.methylFreq PC-T-6_1.methylFreq PC-T-6_2.methylFreq PC-T-7_1.methylFreq PC-T-7_2.methylFreq |  /home/dinh/scripts/methylFreq2BED.pl  20 > /home/kunzhang/CpgMIP/Data/MONOD/MONOD_primary_tumor_RRBS_combined.BED.txt
  cat MONOD_primary_tumor_RRBS_combined.BED.txt | sort -k1,1 -k2,2n > MONOD_primary_tumor_RRBS_combined.sorted.BED.txt  
  • Group the sites covered into clusters:
  ./bed2Clusters.pl MONOD_primary_tumor_RRBS_combined.sorted.BED.txt > MONOD_primary_tumor_RRBS_targets.BED.txt
  • A total of 198,439 regions in the autosomes, with a total size of 48.8Mb.
  • I also compiled a more stringent list of target by requiring a minimal read depth of 100 for the CpG sites. This list has 119,966 autosomal regions, with a total size of 25.8Mb.
 awk ' $5>100 { print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9 }' MONOD_primary_tumor_RRBS_combined.sorted.BED.txt > MONOD_primary_tumor_RRBS_combined_RD100.sorted.BED.txt
 ./bed2Clusters.pl MONOD_primary_tumor_RRBS_combined_RD100.sorted.BED.txt > MONOD_primary_tumor_RRBS_targets_RD100.BED.txt

2. Check the total read depth for CpG sites within the RRBS targets

 cat 6-P-10.methylFreq | /home/dinh/scripts/methylFreq2BED.pl 1 |  /home/kunzhang/softwares/bedtools-2.17.0/bin/bedtools intersect -wa -a - -b /home/kunzhang/CpgMIP/Data/MONOD/MONOD_primary_tumor_RRBS_targets.BED.txt | /home/kunzhang/CpgMIP/Data/MONOD/bed_total_RD.pl
 cat 6-P-10.methylFreq | /home/dinh/scripts/methylFreq2BED.pl 1 | /home/kunzhang/CpgMIP/Data/MONOD/bed_total_RD.pl
 
                   On-target read depth     Total CpG read depth     % on-target
 6-P-10             15,965,319/14,617,979      20,663,932                77.3%
 6-P-1              10,145,982/9,347,247       13,357,081                76.0%
 6-T-1_1            42,122,819/40,519,748      51,404,477                81.9%
 6-T-2_1            39,727,871/38,177,035      47,770,689                83.2%
 CTT-frozen-100ng_1 50,274,814/48,129,684      55,387,277                90.8%
 CTT-frozen-5ng_1   14,303,236/13,695,699      15,310,475                93.4%
 CTT-FFPE-100ng_1   41,695,594/43,527,498      47,409,114                91.8%
 PC-P-1              6,363,821/5,972,702        7,727,128                82.3%
 PC-P-10             7,607,573/7,274,238       9,433,521                 80.6%
 PC-T-1_1           44,275,077/42,697,725     53,030,360                 83.4%
  • I believe the on-target rate is an indicator of the degree of DNA fragmentation. If a fraction of DNA are fragmented, they will be included in the sequencing libraries even without MspI digestion. The on-target rate would be lower in such cases. So plasma DNA always have lower on-target rates.
  • There isn't a big difference between the two sets of target regions identified at different level of stringency. Therefore, I checked the on-target rates for all samples using the larger target set.
Samples All sites RD On-target RD On-target Rate
6-P-10 15,965,319 20,663,932 77.3%
6-P-1 10,145,982 13,357,081 76.0%
6-P-2 6,896,186 9,062,234 76.1%
6-P-3 19,268,289 21,802,969 88.4%
6-P-4 5,974,773 7,495,311 79.7%
6-P-5 12,590,365 16,494,290 76.3%
6-P-6 7,512,063 12,272,372 61.2%
6-P-7 11,161,941 15,365,391 72.6%
6-P-8 10,857,495 14,764,214 73.5%
6-P-9 9,301,655 11,828,425 78.6%
6-T-1_1 42,122,819 51,404,477 81.9%
6-T-1_2 55,332,112 68,680,670 80.6%
6-T-2_1 39,727,871 47,770,689 83.2%
6-T-2_2 43,299,367 53,183,818 81.4%
6-T-3_1 36,204,067 45,527,237 79.5%
6-T-3_2 37,657,579 46,334,499 81.3%
6-T-4_1 48,592,836 59,200,663 82.1%
6-T-4_2 44,582,324 52,297,088 85.2%
6-T-5_1 42,943,331 51,084,769 84.1%
6-T-5_2 51,597,080 62,764,249 82.2%
7-P-10 14,120,252 19,830,640 71.2%
7-P-1 18,418,449 22,370,674 82.3%
7-P-2 3,066,869 4,201,304 73.0%
7-P-3 8,092,483 11,658,106 69.4%
7-P-4 16,586,445 20,741,132 80.0%
7-P-5 8,001,688 10,231,049 78.2%
7-P-6 7,444,769 9,792,239 76.0%
7-P-7 10,815,500 13,957,997 77.5%
7-P-8 10,431,532 14,322,889 72.8%
7-P-9 9,590,120 13,362,917 71.8%
7-T-1_1 55,494,618 68,629,065 80.9%
7-T-1_2 42,296,152 51,527,101 82.1%
7-T-2_1 42,767,775 52,191,210 81.9%
7-T-2_2 45,814,173 57,146,402 80.2%
7-T-3_1 43,288,619 53,666,959 80.7%
7-T-3_2 36,712,682 44,754,376 82.0%
7-T-4_1 33,329,767 45,116,118 73.9%
7-T-4_2 33,219,393 41,142,529 80.7%
7-T-5_1 57,347,061 69,371,414 82.7%
7-T-5_2 44,318,977 54,673,022 81.1%
CTT-FFPE-100ng_1 43,527,498 47,409,114 91.8%
CTT-FFPE-100ng_2 42,796,717 46,523,502 92.0%
CTT-FFPE-5ng_1 17,255,895 18,600,287 92.8%
CTT-FFPE-5ng_2 33,184,997 36,355,026 91.3%
CTT-frozen-100ng_1 50,274,814 55,387,277 90.8%
CTT-frozen-100ng_2 55,941,976 60,903,449 91.9%
CTT-frozen-5ng_1 14,303,236 15,310,475 93.4%
CTT-frozen-5ng_2 31,206,451 34,203,612 91.2%
NC-P-1 22,265,516 25,002,054 89.1%
NC-P-2 16,284,869 19,552,139 83.3%
NC-P-3 27,927,900 31,201,357 89.5%
NC-P-5 9,532,973 12,424,591 76.7%
NC-P-6 13,102,631 15,423,838 85.0%
NC-P-7 12,679,131 15,679,757 80.9%
NC-P-8 13,414,256 16,259,364 82.5%
NC-P-9 13,445,816 15,379,652 87.4%
PC-P-10 7,607,573 9,433,521 80.6%
PC-P-1 6,363,821 7,727,128 82.4%
PC-P-2 13,346,408 15,069,125 88.6%
PC-P-3 3,151,205 4,093,467 77.0%
PC-P-4 6,488,210 8,658,320 74.9%
PC-P-5 5,347,960 6,722,768 79.5%
PC-P-6 6,788,449 8,497,332 79.9%
PC-P-7 6,294,024 8,269,872 76.1%
PC-P-8 4,804,772 5,698,389 84.3%
PC-P-9 3,845,849 4,390,032 87.6%
PC-T-1_1 44,275,077 53,030,360 83.5%
PC-T-1_2 51,315,205 59,257,668 86.6%
PC-T-2_1 68,245,319 80,836,193 84.4%
PC-T-2_2 59,886,740 67,965,891 88.1%
PC-T-4_1 56,814,606 65,283,928 87.0%
PC-T-4_2 53,207,679 61,130,443 87.0%
PC-T-6_1 65,535,581 77,482,439 84.6%
PC-T-6_2 53,192,172 63,933,657 83.2%
PC-T-7_1 48,504,470 56,835,764 85.3%
PC-T-7_2 46,509,004 55,296,800 84.1%

3. Haplotype analysis in the UMR regions

  • The previous analysis that I did was focusing on the LMS clusters, which is a small subset of UMRs. Perhaps the target selection was too strict and I might have missed other informative regions. Therefore, I decided to expand the net and search more broadly.
  • I wrote a script to report all haplotypes in a list of target from a RRBS bam file. I then obtained all haplotypes for all samples within the 49Mb RRBS target regions.
  /home/kunzhang/CpgMIP/Data/MONOD/mergedBam2hapInfo.pl /home/kunzhang/CpgMIP/Data/MONOD/MONOD_primary_tumor_RRBS_targets.BED.txt /media/Ext12T/DD_Ext12T/RRBS_MONOD/Bam_Merged/6-P-10.merged.bam >6-P-10.hapInfo.txt &
  Batch processing script: 1407-combined_expanded_step1_batch_command.sh
  • Then I attempted to classify each haplotype as blood unmethylated haplotype or a cancer methylated haplotype, based on the methylation level of the CpG sites within UMR in the healthy plasma.
    • First I merged all the NC plasma methylFreq files, and created a BED file that represent the average methylation for all the NC plasma samples. Then I use bedtools intersect to extract the sites within blood UMRs.
    • Second I wrote a script to report the following three numbers for each target region: (i) the number of haplotypes; (ii) the number of methylated haplotypes; (iii) the sum negative log-p of all haplotypes. All these numbers for one sample are reported in one regionHap.txt file.
    • Finally I wrote another script to combine the regionHap.txt files for multiple samples, and generate two matrix: HMH_load (the fraction of methylated haplotypes, a measurement of cancer load); HMH_NLP (the sum negative log-p, an indicator for significance)
   /home/kunzhang/CpgMIP/Data/MONOD/regionMethHapAnalysis.pl /home/kunzhang/CpgMIP/Data/MONOD/MONOD_nc_plasma_RRBS_combined_UMRs.BED.txt 6-P-10.hapInfo.txt>6-P-10.regionHap.txt
   Batch processing script: 1407-combined_expanded_step2_batch_command.sh
   ../../report_regionHap_matrix.pl UCSD_MCC_PancreasCancer  PC-P-10.regionHap.txt	PC-P-1.regionHap.txt	PC-P-2.regionHap.txt	PC-P-3.regionHap.txt	PC-P-4.regionHap.txt	PC-P-5.regionHap.txt	PC-P-6.regionHap.txt	PC-P-7.regionHap.txt	PC-P-8.regionHap.txt	PC-P-9.regionHap.txt	PC-T-1_1.regionHap.txt	PC-T-1_2.regionHap.txt	PC-T-2_1.regionHap.txt	PC-T-2_2.regionHap.txt	PC-T-4_1.regionHap.txt	PC-T-4_2.regionHap.txt	PC-T-6_1.regionHap.txt	PC-T-6_2.regionHap.txt	PC-T-7_1.regionHap.txt	PC-T-7_2.regionHap.txt
   ../../report_regionHap_matrix.pl UCSD_MCC_LungCancer 7-P-10.regionHap.txt  7-P-1.regionHap.txt  7-P-2.regionHap.txt  7-P-3.regionHap.txt  7-P-4.regionHap.txt  7-P-5.regionHap.txt  7-P-6.regionHap.txt  7-P-7.regionHap.txt  7-P-8.regionHap.txt  7-P-9.regionHap.txt  7-T-1_1.regionHap.txt  7-T-1_2.regionHap.txt  7-T-2_1.regionHap.txt  7-T-2_2.regionHap.txt  7-T-3_1.regionHap.txt  7-T-3_2.regionHap.txt  7-T-4_1.regionHap.txt  7-T-4_2.regionHap.txt  7-T-5_1.regionHap.txt  7-T-5_2.regionHap.txt
   ../../report_regionHap_matrix.pl UCSD_MCC_ColonCancer 6-P-10.regionHap.txt 6-P-1.regionHap.txt 6-P-2.regionHap.txt 6-P-3.regionHap.txt 6-P-4.regionHap.txt 6-P-5.regionHap.txt 6-P-6.regionHap.txt 6-P-7.regionHap.txt 6-P-8.regionHap.txt 6-P-9.regionHap.txt 6-T-1_1.regionHap.txt 6-T-1_2.regionHap.txt 6-T-2_1.regionHap.txt 6-T-2_2.regionHap.txt 6-T-3_1.regionHap.txt 6-T-3_2.regionHap.txt 6-T-4_1.regionHap.txt 6-T-4_2.regionHap.txt 6-T-5_1.regionHap.txt 6-T-5_2.regionHap.txt
   ../../report_regionHap_matrix.pl Plasma  6-P-10.regionHap.txt  6-P-1.regionHap.txt  6-P-2.regionHap.txt  6-P-3.regionHap.txt  6-P-4.regionHap.txt  6-P-5.regionHap.txt  6-P-6.regionHap.txt  6-P-7.regionHap.txt  6-P-8.regionHap.txt  6-P-9.regionHap.txt  7-P-10.regionHap.txt  7-P-1.regionHap.txt  7-P-2.regionHap.txt  7-P-3.regionHap.txt  7-P-4.regionHap.txt  7-P-5.regionHap.txt  7-P-6.regionHap.txt  7-P-7.regionHap.txt  7-P-8.regionHap.txt  7-P-9.regionHap.txt  NC-P-1.regionHap.txt  NC-P-2.regionHap.txt  NC-P-3.regionHap.txt  NC-P-5.regionHap.txt  NC-P-6.regionHap.txt  NC-P-7.regionHap.txt  NC-P-8.regionHap.txt  NC-P-9.regionHap.txt  PC-P-10.regionHap.txt  PC-P-1.regionHap.txt  PC-P-2.regionHap.txt  PC-P-3.regionHap.txt  PC-P-4.regionHap.txt  PC-P-5.regionHap.txt  PC-P-6.regionHap.txt  PC-P-7.regionHap.txt  PC-P-8.regionHap.txt  PC-P-9.regionHap.txt
   ../../report_regionHap_matrix.pl Primary_tumor  6-T-1_1.regionHap.txt 6-T-1_2.regionHap.txt 6-T-2_1.regionHap.txt 6-T-2_2.regionHap.txt 6-T-3_1.regionHap.txt 6-T-3_2.regionHap.txt 6-T-4_1.regionHap.txt 6-T-4_2.regionHap.txt 6-T-5_1.regionHap.txt 6-T-5_2.regionHap.txt 7-T-1_1.regionHap.txt 7-T-1_2.regionHap.txt 7-T-2_1.regionHap.txt 7-T-2_2.regionHap.txt 7-T-3_1.regionHap.txt 7-T-3_2.regionHap.txt 7-T-4_1.regionHap.txt 7-T-4_2.regionHap.txt 7-T-5_1.regionHap.txt 7-T-5_2.regionHap.txt PC-T-1_1.regionHap.txt PC-T-1_2.regionHap.txt PC-T-2_1.regionHap.txt PC-T-2_2.regionHap.txt PC-T-4_1.regionHap.txt PC-T-4_2.regionHap.txt PC-T-6_1.regionHap.txt PC-T-6_2.regionHap.txt PC-T-7_1.regionHap.txt PC-T-7_2.regionHap.txt CTT-FFPE-100ng_1.regionHap.txt CTT-FFPE-100ng_2.regionHap.txt CTT-FFPE-5ng_1.regionHap.txt CTT-FFPE-5ng_2.regionHap.txt CTT-frozen-100ng_1.regionHap.txt CTT-frozen-100ng_2.regionHap.txt CTT-frozen-5ng_1.regionHap.txt CTT-frozen-5ng_2.regionHap.txt
  • To take a quick look at the HMH load in the primary tumors, I took the Primary_tumor.HMH_load.matrix file, manually filtered out the regions that have average load below 0.1 (just to reduce the number of regions from 20k+ to ~5k. Then made a heatmap with R (note that I removed the Biochain CTT samples to avoid variability of sample source). Based on the heatmap, there seems to be a decent level of heterogeneity even for the samples of the same tumor types.
 library(gplots)
 x=read.table("MCC_Primary_tumor_mean_0.1_up.HMH_load.matrix.txt",header=TRUE,row.names=1);
 heatmap.2(as.matrix(x), col=bluered(128), scale="none", cexCol=0.8, cexRow=0.1 ,key=T, symkey=F,density.info="histogram",trace="none",dendrogram="both")
 File:MCC primary tumor mean 0.1 up heatmap.png

4. Mix hap analysis

  • It looks like simple classification of a haplotype into blood or tumor based on the number of methylated CpG is a little bit crude, and it's hard to come up with a robust strategy for detecting tumor haplotypes in patient's plasma.
  • I decided to go back to the mix hap analysis and make a number of adjustments.
    • My previous version of mix hap analysis assumes BSPP data, in which most or all reads overlap perfectly. With RRBS data, many reads in a target region don't overlap. So I needed to be more flexible in comparing haplotypes.
    • The average read depth is lower in RRBS data. For CpG sites that have low or no coverage in normal plasma, I need to use the methylation level of whole blood WGBS data to estimate the haplotype frequency.
    • I decided to switch to LOD (log of odds ratio) between normal plasma and tumor. Essentially, I created two reference hapInfo file, one for all normal plasma (NC-P) samples combined, another for all primary tumor samples of one cancer type (PC, CC or LC) combined. For each sequence reads in the patient plasma, I asked what's the probability that the haplotype came from blood and the probability it came from tumor. Then I took the log of the odds ratio between these two probabilities, and sum them up for each target region and for the entire genome.
  • I started with all RRBS targets within the UMRs, but then need to trim the target list, based on which target is informative. Since the sample size is limited, doing cross validation is not realistic. So I simply pooled all patient plasma data (treating PC, CC, LC separately), and ran mixMethHapAnalysis_19Aug14.pl against pooled normal plasma and pooled tumor to identify the regions that have the highest LOD scores.
  • To evaluate the level of false positives due to variability in normal plasma, I took one normal plasma as a test sample to compare with the pool of the other seven normal plasma. I did this for each of the 8 plasma samples, and reported a list of regions that show high LOD (false positives). These regions were excluded from the list of informative markers.