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

From ZhangLabWiki
Jump to navigation Jump to search
 
(19 intermediate revisions by the same user not shown)
Line 219: Line 219:
*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.
*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.
  ../../mixMethHapAnalysis_19Aug14.pl  NC-P-ALL.hapInfo.txt 6-T-ALL.hapInfo.txt 6-P-ALL.hapInfo.txt /home/kunzhang/CpgMIP/MONOD/Public_data/WB_middle_age_UMR.BED.txt | sort -k 2,2nr > RRBS_targets_LOD_CC.txt &
  ../../mixMethHapAnalysis_19Aug14.pl  NC-P-ALL.hapInfo.txt 6-T-ALL.hapInfo.txt 6-P-ALL.hapInfo.txt /home/kunzhang/CpgMIP/MONOD/Public_data/WB_middle_age_UMR.BED.txt | sort -k 2,2nr > RRBS_targets_LOD_CC.txt &
  ../../mixMethHapAnalysis_19Aug14.pl  NC-P-ALL.hapInfo.txt 6-T-ALL.hapInfo.txt 6-P-ALL.hapInfo.txt /home/kunzhang/CpgMIP/MONOD/Public_data/WB_middle_age_UMR.BED.txt  | sort -k 2,2nr  > RRBS_targets_LOD_CC.txt &
  ../../mixMethHapAnalysis_19Aug14.pl  NC-P-ALL.hapInfo.txt 7-T-ALL.hapInfo.txt 7-P-ALL.hapInfo.txt /home/kunzhang/CpgMIP/MONOD/Public_data/WB_middle_age_UMR.BED.txt  | sort -k 2,2nr  > RRBS_targets_LOD_LC.txt &
  ../../mixMethHapAnalysis_19Aug14.pl  NC-P-ALL.hapInfo.txt 6-T-ALL.hapInfo.txt 6-P-ALL.hapInfo.txt /home/kunzhang/CpgMIP/MONOD/Public_data/WB_middle_age_UMR.BED.txt | sort -k 2,2nr > RRBS_targets_LOD_CC.txt &
  ../../mixMethHapAnalysis_19Aug14.pl  NC-P-ALL.hapInfo.txt PC-T-ALL.hapInfo.txt PC-P-ALL.hapInfo.txt /home/kunzhang/CpgMIP/MONOD/Public_data/WB_middle_age_UMR.BED.txt | sort -k 2,2nr > RRBS_targets_LOD_PC.txt &
 


*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. It appeared that whether there are false positive calls in a region also depends on which tumor reference file is used. For now, I just identified these regions for each of the three tumor types separately. In the future, we might need to process a lot of normal plasma and primary tumor samples to fully cover the biological variability.
*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. It appeared that whether there are false positive calls in a region also depends on which tumor reference file is used. For now, I just identified these regions for each of the three tumor types separately. In the future, we might need to process a lot of normal plasma and primary tumor samples to fully cover the biological variability.
Line 228: Line 227:
     ./filter_targets.pl  RRBS_targets_LOD_PC.txt NC-P_high-PC-LOD_targets.txt  | sort -k 2,2nr  > RRBS_targets_LOD_PC_top_347.txt
     ./filter_targets.pl  RRBS_targets_LOD_PC.txt NC-P_high-PC-LOD_targets.txt  | sort -k 2,2nr  > RRBS_targets_LOD_PC_top_347.txt
     ./filter_targets.pl  RRBS_targets_LOD_LC.txt NC-P_high-LC-LOD_targets.txt | sort -k 2,2nr > RRBS_targets_LOD_LC_top_346.txt
     ./filter_targets.pl  RRBS_targets_LOD_LC.txt NC-P_high-LC-LOD_targets.txt | sort -k 2,2nr > RRBS_targets_LOD_LC_top_346.txt
     ./filter_targets.pl  RRBS_targets_LOD_LC.txt NC-P_high-CC-LOD_targets.txt  | sort -k 2,2nr >  RRBS_targets_LOD_CC_top_302.txt
     ./filter_targets.pl  RRBS_targets_LOD_CC.txt NC-P_high-CC-LOD_targets.txt  | sort -k 2,2nr >  RRBS_targets_LOD_CC_top_302.txt
*Then I searched for the tumor haplotypes in each plasma samples within the marker regions.
*Then I searched for the tumor haplotypes in each plasma samples within the marker regions.
     ./get_all_CC-P_LOD.pl > CC-P_LOD_out.txt &
     ./get_all_CC-P_LOD.pl > CC-P_LOD_out.txt &
Line 237: Line 236:
**Tumor loads are detectible in patient's plasma, and the load varies from one sample to another, more so than from one tumor type to another.
**Tumor loads are detectible in patient's plasma, and the load varies from one sample to another, more so than from one tumor type to another.
**I was only able to achieve this after removing all regions that contain false-positive haplotypes in normal plasma. Cross validation is not possible with so few normal plasma samples. So we need more normal plasma data to have a real assessment of the false positive rate.
**I was only able to achieve this after removing all regions that contain false-positive haplotypes in normal plasma. Cross validation is not possible with so few normal plasma samples. So we need more normal plasma data to have a real assessment of the false positive rate.
**Are there the regions that contain tumor haplotypes across multiple samples?? A total 241 regions were found in 2 or more samples. The regions that appeared in 5 or more samples are listed below.
**Are there the regions that contain tumor haplotypes across multiple samples?? A total 241 regions were found in 2 or more samples. These regions are enriched for DNA binding proteins, especially HOX genes.
  [[Image:Common_high_LOD_regions_InterPro_enrichment.png|800px]]
  [[Image:Common_high_LOD_regions_GO_MolecularFunction_enrichment.png|800px]]
**The regions that appeared in 5 or more samples are listed below.


{| {{table}}
{| {{table}}
Line 270: Line 272:
| chr2:177036288-177036973||7||HOXD3
| chr2:177036288-177036973||7||HOXD3
|-
|-
| chr22:19710616-19712149||7||SETP5
| chr22:19710616-19712149||7||GP1BB
|-
|-
| chr11:555351-556680||6||LMNTD2
| chr11:555351-556680||6||LMNTD2
Line 330: Line 332:
|  
|  
|}
|}
*A region related to SPN gene
                Normal Plasma                                            Primary Tumor(LC-4)                          Matched Patient Plasma(LC-4)
  [[Image:chr16.29675696-29676014_SPN_NCP.png|400px]][[Image:chr16.29675696-29676014_SPN_LC-T4.png|400px]][[Image:chr16.29675696-29676014_SPN_LC-P4.png|200px]]
====Repeat the analysis using a refined set of targets====
*The analyses above were performed based on a set of RRBS targets compiled based on the read depth of CpG sites in tumor data. Some of the target regions are fairly large, and contain many non-overlapping reads that cannot be directly compared in the haplotype analysis.
*A different way to define the RRBS targets is based on the read depth distribution (not just CpG sites), which might better separate non-overlapping reads.
*To test this idea, I compiled a new set of RRBS targets with the following steps.
  samtools merge /home/kunzhang/CpgMIP/Data/MONOD/1407-combined_RRBS/Primary_tumor_ALL.merged.bam 6-T-1_1.merged.bam 6-T-1_2.merged.bam 6-T-2_1.merged.bam 6-T-2_2.merged.bam 6-T-3_1.merged.bam 6-T-3_2.merged.bam 6-T-4_1.merged.bam 6-T-4_2.merged.bam 6-T-5_1.merged.bam 6-T-5_2.merged.bam 7-T-1_1.merged.bam 7-T-1_2.merged.bam 7-T-2_1.merged.bam 7-T-2_2.merged.bam 7-T-3_1.merged.bam 7-T-3_2.merged.bam 7-T-4_1.merged.bam 7-T-4_2.merged.bam 7-T-5_1.merged.bam 7-T-5_2.merged.bam PC-T-1_1.merged.bam PC-T-1_2.merged.bam PC-T-2_1.merged.bam PC-T-2_2.merged.bam PC-T-4_1.merged.bam PC-T-4_2.merged.bam PC-T-6_1.merged.bam PC-T-6_2.merged.bam PC-T-7_1.merged.bam PC-T-7_2.merged.bam
  samtools sort Primary_tumor_ALL.merged.bam Primary_tumor_ALL.sorted &
  nohup /home/kunzhang/softwares/bedtools-2.17.0/bin/bedtools genomecov -bg -split -ibam Primary_tumor_ALL.sorted.bam >    Pimary_tumor_ALL.genomecov.bed &
  awk '$4>50 { print $1"\t"$2"\t"$3}'  Primary_tumor_ALL.genomecov.bed | /home/kunzhang/softwares/bedtools-2.17.0/bin/bedtools merge -d 10 -i - > Primary_tumor_ALL.genomecov.RD50.merged.bed
  awk '$3-$2>80 {print $1"\t"$2"\t"$3"\t"$3-$2+1}' Primary_tumor_ALL.genomecov.RD50.merged.bed> Primary_tumor_ALL.genomecov.RD50_80UP.merged.bed
  sh [[media:1407-combined_expanded_31Aug14_step1_batch_command.txt|1407-combined_expanded_31Aug14_step1_batch_command.sh]] &
  ../../mixMethHapAnalysis_19Aug14.pl  NC-P-ALL.hapInfo.txt 6-T-ALL.hapInfo.txt 6-P-ALL.hapInfo.txt /home/kunzhang/CpgMIP/MONOD/Public_data/WB_middle_age_UMR.BED.txt | sort -k 2,2nr > RRBS_targets_LOD_CC.txt &
  ../../mixMethHapAnalysis_19Aug14.pl  NC-P-ALL.hapInfo.txt 7-T-ALL.hapInfo.txt 7-P-ALL.hapInfo.txt /home/kunzhang/CpgMIP/MONOD/Public_data/WB_middle_age_UMR.BED.txt  | sort -k 2,2nr  > RRBS_targets_LOD_LC.txt &
  ../../mixMethHapAnalysis_19Aug14.pl  NC-P-ALL.hapInfo.txt PC-T-ALL.hapInfo.txt PC-P-ALL.hapInfo.txt /home/kunzhang/CpgMIP/MONOD/Public_data/WB_middle_age_UMR.BED.txt | sort -k 2,2nr > RRBS_targets_LOD_PC.txt &
  ./find_NC-P_high-LOD_targets.pl > NC-P_high-LOD_targets.txt
    ./filter_targets.pl  RRBS_targets_LOD_PC.txt NC-P_high-PC-LOD_targets.txt  | sort -k 2,2nr  > RRBS_targets_LOD_PC_top_1914.txt
    ./filter_targets.pl  RRBS_targets_LOD_LC.txt NC-P_high-LC-LOD_targets.txt | sort -k 2,2nr > RRBS_targets_LOD_LC_top_1629.txt
    ./filter_targets.pl  RRBS_targets_LOD_LC.txt NC-P_high-CC-LOD_targets.txt  | sort -k 2,2nr >  RRBS_targets_LOD_CC_top_1487.txt
    ./get_all_CC-P_LOD.pl > CC-P_LOD_out.txt &
    ./get_all_LC-P_LOD.pl > LC-P_LOD_out.txt &
    ./get_all_PC-P_LOD.pl > PC-P_LOD_out.txt &
*A total of 822 regions that contain cancer haplotypes in two or more samples were identified.
    [[Image:Colon_cancer_UMR_tumor_haplotype_load_31Aug14.png|350px]][[Image:Lung_cancer_UMR_tumor_haplotype_load_31Aug14.png|350px]][[Image:Pancreatic_cancer_UMR_tumor_haplotype_load_31Aug14.png|350px]]
    [[Image:CC_LC_PC_common_Go_biological_process_31Aug14.png|1200px]]
    [[Image:CC_LC_PC_common_Go_molecular_function_31Aug14.png|1200px]]
    [[Image:CC_LC_PC_common_InterPro_31Aug14.png|800px]]
    [[Image:CC_LC_PC_common_MSigDB_cancer_neighborhood_31Aug14.png|800px]]


==5. Searching for primary tumor specific haplotypes==
==5. Searching for primary tumor specific haplotypes==
Line 339: Line 374:
   [[Image:Primary_tumor_specific_haplotypes.png|400px]]
   [[Image:Primary_tumor_specific_haplotypes.png|400px]]
*The tumor loads seems to be more consistently across different samples with this analysis. Will need to zoom into specific regions, and confirm the haplotype visually.
*The tumor loads seems to be more consistently across different samples with this analysis. Will need to zoom into specific regions, and confirm the haplotype visually.
              Normal Plasma
  [[Image:chr18.19756471-19757245_GATA6_NCP.png|350px]]
            LC Patient #4 plasma                              LC Patient #4 primary tumor
  [[Image:chr18.19756471-19757245_GATA6_LC-P4.png|350px]][[Image:chr18.19756471-19757245_GATA6_LC-T4.png|350px]]
            PC Patient #7 plasma                              PC Patient #7 primary tumor
  [[Image:chr18.19756471-19757245_GATA6_PC-P7.png|350px]][[Image:chr18.19756471-19757245_GATA6_PC-T7.png|350px]]

Latest revision as of 19:37, 26 October 2014

RRBS data analysis (continued)[edit]

1. Compile a list of RRBS targets[edit]

  • 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[edit]

 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[edit]

  • 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[edit]

  • 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.
../../mixMethHapAnalysis_19Aug14.pl  NC-P-ALL.hapInfo.txt 6-T-ALL.hapInfo.txt 6-P-ALL.hapInfo.txt /home/kunzhang/CpgMIP/MONOD/Public_data/WB_middle_age_UMR.BED.txt | sort -k 2,2nr > RRBS_targets_LOD_CC.txt &
../../mixMethHapAnalysis_19Aug14.pl  NC-P-ALL.hapInfo.txt 7-T-ALL.hapInfo.txt 7-P-ALL.hapInfo.txt /home/kunzhang/CpgMIP/MONOD/Public_data/WB_middle_age_UMR.BED.txt  | sort -k 2,2nr  > RRBS_targets_LOD_LC.txt &
../../mixMethHapAnalysis_19Aug14.pl  NC-P-ALL.hapInfo.txt PC-T-ALL.hapInfo.txt PC-P-ALL.hapInfo.txt /home/kunzhang/CpgMIP/MONOD/Public_data/WB_middle_age_UMR.BED.txt | sort -k 2,2nr > RRBS_targets_LOD_PC.txt &
  • 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. It appeared that whether there are false positive calls in a region also depends on which tumor reference file is used. For now, I just identified these regions for each of the three tumor types separately. In the future, we might need to process a lot of normal plasma and primary tumor samples to fully cover the biological variability.
  ./find_NC-P_high-LOD_targets.pl > NC-P_high-LOD_targets.txt
  • I manually split the NC-P_high-LOD_targets.txt file into three based on the tumor types, then extracted the informative marker lists based on the differences of LOD scores between patient plasma and normal plasma.
   ./filter_targets.pl  RRBS_targets_LOD_PC.txt NC-P_high-PC-LOD_targets.txt  | sort -k 2,2nr  > RRBS_targets_LOD_PC_top_347.txt
   ./filter_targets.pl  RRBS_targets_LOD_LC.txt NC-P_high-LC-LOD_targets.txt | sort -k 2,2nr > RRBS_targets_LOD_LC_top_346.txt
   ./filter_targets.pl  RRBS_targets_LOD_CC.txt NC-P_high-CC-LOD_targets.txt  | sort -k 2,2nr >  RRBS_targets_LOD_CC_top_302.txt
  • Then I searched for the tumor haplotypes in each plasma samples within the marker regions.
   ./get_all_CC-P_LOD.pl > CC-P_LOD_out.txt &
   ./get_all_LC-P_LOD.pl > LC-P_LOD_out.txt &
   ./get_all_PC-P_LOD.pl > PC-P_LOD_out.txt &
 File:Colon cancer UMR tumor haplotype load.pngFile:Lung cancer UMR tumor haplotype load.pngFile:Pancreatic cancer UMR tumor haplotype load.png
  • Some obervations:
    • Tumor loads are detectible in patient's plasma, and the load varies from one sample to another, more so than from one tumor type to another.
    • I was only able to achieve this after removing all regions that contain false-positive haplotypes in normal plasma. Cross validation is not possible with so few normal plasma samples. So we need more normal plasma data to have a real assessment of the false positive rate.
    • Are there the regions that contain tumor haplotypes across multiple samples?? A total 241 regions were found in 2 or more samples. These regions are enriched for DNA binding proteins, especially HOX genes.
 File:Common high LOD regions InterPro enrichment.png
 File:Common high LOD regions GO MolecularFunction enrichment.png
    • The regions that appeared in 5 or more samples are listed below.
Target_region # samples Gene
chr16:29675696-29676014 13 SPN
chr11:72295386-72295938 12 PDE2A
chr16:67428179-67429159 12 ZDHHC1
chr17:76127129-76128704 12 TMC6
chr2:177016416-177016749 12 HOXD4
chr19:1074331-1075173 10 HMHA1
chr6:100897213-100897619 10 SIM1
chr2:127839497-127839766 8 BIN1
chr22:20784703-20786333 8 SCARF2
chr6:42071960-42072652 8 C6orf132
chr12:124941780-124941976 7 NCOR2
chr19:18980616-18981821 7 GDF1/CERS1
chr2:177036288-177036973 7 HOXD3
chr22:19710616-19712149 7 GP1BB
chr11:555351-556680 6 LMNTD2
chr12:49942621-49943178 6 KCNH3
chr14:37127773-37128202 6 PAX9
chr16:54970220-54970503 6 CRNDE/IRX5
chr17:5404000-5404600 6 LOC728392
chr2:111876191-111877448 6 ACOXL
chr4:169798910-169799623 6 PALLD
chr7:127671170-127671464 6 SND1
chr7:155246283-155248093 6 EN2
chr8:80804055-80804317 6 GDF6
chr8:97156647-97158028 6 GDF6
chr10:118030474-118031344 5 GFRA1
chr1:151810190-151811682 5 C2CD4D
chr12:115105189-115105423 5 TBX3 (downstream)
chr1:225865165-225865313 5 ENAH
chr13:112547537-112548375 5 LINC00354
chr15:96895306-96895504 5 NR2F2 (downstream)
chr16:86612197-86613526 5 FOXL1
chr17:70111952-70112499 5 SOX9-AS1
chr1:91300884-91301262 5 #not near a gene#
chr19:4328686-4329420 5 STAP2
chr21:36900927-36901032 5 LOC100506403
chr2:171678137-171678305 5 GAD1
chr2:177023686-177024305 5 HOXD3/4
chr2:233925159-233925345 5 INPP5D
chr5:88184868-88185224 5 MEF2C
chr6:10385059-10385426 5 TFAP2A
chr8:142528135-142529145 5 MROH5
  • A region related to SPN gene
               Normal Plasma                                            Primary Tumor(LC-4)                           Matched Patient Plasma(LC-4)
 File:Chr16.29675696-29676014 SPN NCP.pngFile:Chr16.29675696-29676014 SPN LC-T4.pngFile:Chr16.29675696-29676014 SPN LC-P4.png

Repeat the analysis using a refined set of targets[edit]

  • The analyses above were performed based on a set of RRBS targets compiled based on the read depth of CpG sites in tumor data. Some of the target regions are fairly large, and contain many non-overlapping reads that cannot be directly compared in the haplotype analysis.
  • A different way to define the RRBS targets is based on the read depth distribution (not just CpG sites), which might better separate non-overlapping reads.
  • To test this idea, I compiled a new set of RRBS targets with the following steps.
  samtools merge /home/kunzhang/CpgMIP/Data/MONOD/1407-combined_RRBS/Primary_tumor_ALL.merged.bam 6-T-1_1.merged.bam 6-T-1_2.merged.bam 6-T-2_1.merged.bam 6-T-2_2.merged.bam 6-T-3_1.merged.bam 6-T-3_2.merged.bam 6-T-4_1.merged.bam 6-T-4_2.merged.bam 6-T-5_1.merged.bam 6-T-5_2.merged.bam 7-T-1_1.merged.bam 7-T-1_2.merged.bam 7-T-2_1.merged.bam 7-T-2_2.merged.bam 7-T-3_1.merged.bam 7-T-3_2.merged.bam 7-T-4_1.merged.bam 7-T-4_2.merged.bam 7-T-5_1.merged.bam 7-T-5_2.merged.bam PC-T-1_1.merged.bam PC-T-1_2.merged.bam PC-T-2_1.merged.bam PC-T-2_2.merged.bam PC-T-4_1.merged.bam PC-T-4_2.merged.bam PC-T-6_1.merged.bam PC-T-6_2.merged.bam PC-T-7_1.merged.bam PC-T-7_2.merged.bam
  samtools sort Primary_tumor_ALL.merged.bam Primary_tumor_ALL.sorted &
  nohup /home/kunzhang/softwares/bedtools-2.17.0/bin/bedtools genomecov -bg -split -ibam Primary_tumor_ALL.sorted.bam >    Pimary_tumor_ALL.genomecov.bed &
  awk '$4>50 { print $1"\t"$2"\t"$3}'  Primary_tumor_ALL.genomecov.bed | /home/kunzhang/softwares/bedtools-2.17.0/bin/bedtools merge -d 10 -i - > Primary_tumor_ALL.genomecov.RD50.merged.bed
  awk '$3-$2>80 {print $1"\t"$2"\t"$3"\t"$3-$2+1}' Primary_tumor_ALL.genomecov.RD50.merged.bed> Primary_tumor_ALL.genomecov.RD50_80UP.merged.bed
  sh 1407-combined_expanded_31Aug14_step1_batch_command.sh &
  ../../mixMethHapAnalysis_19Aug14.pl  NC-P-ALL.hapInfo.txt 6-T-ALL.hapInfo.txt 6-P-ALL.hapInfo.txt /home/kunzhang/CpgMIP/MONOD/Public_data/WB_middle_age_UMR.BED.txt | sort -k 2,2nr > RRBS_targets_LOD_CC.txt &
  ../../mixMethHapAnalysis_19Aug14.pl  NC-P-ALL.hapInfo.txt 7-T-ALL.hapInfo.txt 7-P-ALL.hapInfo.txt /home/kunzhang/CpgMIP/MONOD/Public_data/WB_middle_age_UMR.BED.txt  | sort -k 2,2nr  > RRBS_targets_LOD_LC.txt &
  ../../mixMethHapAnalysis_19Aug14.pl  NC-P-ALL.hapInfo.txt PC-T-ALL.hapInfo.txt PC-P-ALL.hapInfo.txt /home/kunzhang/CpgMIP/MONOD/Public_data/WB_middle_age_UMR.BED.txt | sort -k 2,2nr > RRBS_targets_LOD_PC.txt &
  ./find_NC-P_high-LOD_targets.pl > NC-P_high-LOD_targets.txt
   ./filter_targets.pl  RRBS_targets_LOD_PC.txt NC-P_high-PC-LOD_targets.txt  | sort -k 2,2nr  > RRBS_targets_LOD_PC_top_1914.txt
   ./filter_targets.pl  RRBS_targets_LOD_LC.txt NC-P_high-LC-LOD_targets.txt | sort -k 2,2nr > RRBS_targets_LOD_LC_top_1629.txt
   ./filter_targets.pl  RRBS_targets_LOD_LC.txt NC-P_high-CC-LOD_targets.txt  | sort -k 2,2nr >  RRBS_targets_LOD_CC_top_1487.txt
   ./get_all_CC-P_LOD.pl > CC-P_LOD_out.txt &
   ./get_all_LC-P_LOD.pl > LC-P_LOD_out.txt &
   ./get_all_PC-P_LOD.pl > PC-P_LOD_out.txt &
  • A total of 822 regions that contain cancer haplotypes in two or more samples were identified.
    File:Colon cancer UMR tumor haplotype load 31Aug14.pngFile:Lung cancer UMR tumor haplotype load 31Aug14.pngFile:Pancreatic cancer UMR tumor haplotype load 31Aug14.png
    File:CC LC PC common Go biological process 31Aug14.png
    File:CC LC PC common Go molecular function 31Aug14.png
    File:CC LC PC common InterPro 31Aug14.png
    File:CC LC PC common MSigDB cancer neighborhood 31Aug14.png

5. Searching for primary tumor specific haplotypes[edit]

  • In the analysis above, I used the pooled normal plasma and pooled primary tumor data as the negative and positive controls. For each haplotype from a patient's plasma, I determined the likelihood of it being present in normal plasma versus primary tumor, then calculated the tumor haplotype load.
  • The other approach, assuming the matched primary tumor is available, is to look for the haplotypes in the patient's plasma that exactly matched the ones seen in primary tumor and not present in normal plasma. This could be more specific as we make fewer assumption on what haplotype should look like tumor haplotypes.
  • For this purpose, I used a modified version of mixMethHapAnalysis script to identify haplotypes that have near exact match to primary tumor haplotypes. mixMethHapAnalysisStringent_27Aug14.pl
  • Note that in this case in expanded the search across all RRBS targets, beyond the UMR regions used in the previous analysis.
  ./get_tumor_haps.sh &
  File:Primary tumor specific haplotypes.png
  • The tumor loads seems to be more consistently across different samples with this analysis. Will need to zoom into specific regions, and confirm the haplotype visually.
              Normal Plasma
 File:Chr18.19756471-19757245 GATA6 NCP.png
           LC Patient #4 plasma                              LC Patient #4 primary tumor
 File:Chr18.19756471-19757245 GATA6 LC-P4.pngFile:Chr18.19756471-19757245 GATA6 LC-T4.png
           PC Patient #7 plasma                              PC Patient #7 primary tumor
 File:Chr18.19756471-19757245 GATA6 PC-P7.pngFile:Chr18.19756471-19757245 GATA6 PC-T7.png