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

From ZhangLabWiki
Jump to navigation Jump to search
Line 198: Line 198:
**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.  
**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)
**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: [[Media: 1407-combined_expanded_step2_batch_command.txt|1407-combined_expanded_step2_batch_command.sh]]

Revision as of 18:44, 14 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