Kun:LabNotes/MONOD/2014-8-12
Jump to navigation
Jump to search
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.