Kun:LabNotes/MONOD/2013-11-22
Jump to navigation
Jump to search
MONOD round 1
Target identification
- Data used:
- Cancer data:
- GBM: U87 (ENCODE RRBS Hudson Alpha; ENCODE 450k)
- Pancreatic cancer: PANC1 (ENCODE RRBS UW & Hudson Alpha; ENCODE 450k)
- Whole blood data: GSE30253 RRBS data; GSE31263 WGBS data
- Cancer data:
- DMSs and DMS clusters identification:
- I wrote a simple script that takes the average of all existing data for each CpG site, and reported the sites at which the methylation difference was greater than 0.8 between cancer and whole blood.
./find_DMS_MONOD_v1.pl > MONOD_v1_DMS.txt
- These sites were then grouped into DMS clusters.
./extract_clusters.pl MONOD_v1_DMS.txt > MONOD_v1_DMS_clusters.txt
Probe design
- Dinh went through one round of design using ppDesigner v1.1. However, upon close inspection, I found that some capturing arms contain long stretches of poly-Ts or poly-Gs.
- I modified get_probes.pl, and filtered out any candidate capturing arm that has over 6 A/T/G/Cs.
- I repeated the probe design.
- PANC_GBM_DMS: Target list
- Set I: gap size = 400-450bp, flanking len = 250bp; ppDesigner Output + probe assembly script V6 = Probe sequences: 3,514 probes
- Set II: gap size = 125-175bp, flanking len = 125bp; ppDesigner Output + probe assembly script V4 = Probe sequences: 6,960 probes
- CRC_DMS: Target list.
- Set I: gap size = 400-450bp, flanking len = 300bp; ppDesigner Output + probe assembly script V6 = Probe sequences: 148 probes
- Set II: gap size = 125-175bp, flanking len = 175bp; ppDesigner Output + probe assembly script V4 = Probe sequences: 149 probes
- PANC_GBM_DMS: Target list
- MONOD v1 oligo pool ordered from Custom Array.
BSPP capture
- Alan performed the first experiment on five cancer cell lines and three whole blood samples. Note that Alan also called this probe set as GP1:
Data analysis
- All 16 sequencing libraries were sequenced in a PE150bp MiSeq run: /home/kunzhang/seqStore/131227_MiSeq_GP1
- Pre-processing of sequencing reads: the 6bp UMI in each of the Read 1 is extracted and placed in the read name, then 27bp of the 5'-ends (corresponding to H1&H2) were trimmed off from Read 1 and Read 2.
./extract_UMI_PE.pl
- The reads were then mapped by bisReadMapperPE19.pl, and the bam files were used for haplotype analysis. The 400bp probe set had very poor capture, so the following analysis focus on the data generated with the 150bp probe set (GP1_V4).
Sample | raw reads | mapped reads | on-target reads | mapping rate | specificity | # CpG called |
BE2C_V4 | 1,386,590 | 915,834 | 862,623 | 62.21% | 94.19% | 31,853 |
BE2C_V6 | 2,933,904 | 489,918 | 38,994 | 1.33% | 7.96% | 28,819 |
BXPC3_V4 | 1,805,276 | 1,256,196 | 1,232,955 | 68.30% | 98.15% | 36,528 |
BXPC3_V6 | 2,753,274 | 376,020 | 24,332 | 0.88% | 6.47% | 18,885 |
PANC1_V4 | 1,264,618 | 807,262 | 763,920 | 60.41% | 94.63% | 32,712 |
PANC1_V6 | 2,522,344 | 386,410 | 29,667 | 1.18% | 7.68% | 22,154 |
T98G_V4 | 1,803,782 | 1,125,066 | 1,107,621 | 61.41% | 98.45% | 37,906 |
T98G_V6 | 3,151,530 | 478,980 | 30,172 | 0.96% | 6.30% | 27,995 |
U87MG_V4 | 1,586,144 | 1,047,246 | 1,011,039 | 63.74% | 96.54% | 34,562 |
U87MG_V6 | 2,856,720 | 419,040 | 40,911 | 1.43% | 9.76% | 26,848 |
UCLA-SZ_B1_V4 | 1,427,496 | 774,492 | 713,916 | 50.01% | 92.18% | 29,093 |
UCLA-SZ_B1_V6 | 2,173,094 | 374,268 | 21,959 | 1.01% | 5.87% | 21,690 |
UCLA-SZ_D1_V4 | 1,409,960 | 523,588 | 472,779 | 33.53% | 90.30% | 26,859 |
UCLA-SZ_D1_V6 | 1,125,166 | 136,234 | 5,584 | 0.50% | 4.10% | 8,815 |
UCLA-SZ_H11_V4 | 1,574,628 | 643,542 | 583,033 | 37.03% | 90.60% | 29,864 |
UCLA-SZ_H11_V6 | 2,165,268 | 268,476 | 15,822 | 0.73% | 5.89% | 16,206 |
- I wrote the bam2hapInfo.pl script to extract all haplotypes from a pair of bam files within all targeted regions.
../bam2hapInfo.pl /home/kunzhang/CpgMIP/MONOD/GBM_PC_CRC_V4_capture_regions_hg19.txt BE2C_V4_R1_UMI_001.fastq.sorted.fwd.bam BE2C_V4_R1_UMI_001.fastq.sorted.rev.bam >BE2C_V4_hapInfo.txt ../bam2hapInfo.pl /home/kunzhang/CpgMIP/MONOD/GBM_PC_CRC_V4_capture_regions_hg19.txt PANC1_V4_R1_UMI_001.fastq.sorted.fwd.bam PANC1_V4_R1_UMI_001.fastq.sorted.rev.bam >PANC1_V4_hapInfo.txt ../bam2hapInfo.pl /home/kunzhang/CpgMIP/MONOD/GBM_PC_CRC_V4_capture_regions_hg19.txt T98G_V4_R1_UMI_001.fastq.sorted.fwd.bam T98G_V4_R1_UMI_001.fastq.sorted.rev.bam >T98G_V4_hapInfo.txt ../bam2hapInfo.pl /home/kunzhang/CpgMIP/MONOD/GBM_PC_CRC_V4_capture_regions_hg19.txt U87MG_V4_R1_UMI_001.fastq.sorted.fwd.bam U87MG_V4_R1_UMI_001.fastq.sorted.rev.bam >U87MG_V4_hapInfo.txt ../bam2hapInfo.pl /home/kunzhang/CpgMIP/MONOD/GBM_PC_CRC_V4_capture_regions_hg19.txt UCLA-SZ_B1_V4_R1_UMI_001.fastq.sorted.fwd.bam UCLA-SZ_B1_V4_R1_UMI_001.fastq.sorted.rev.bam >UCLA-SZ_B1_V4_hapInfo.txt ../bam2hapInfo.pl /home/kunzhang/CpgMIP/MONOD/GBM_PC_CRC_V4_capture_regions_hg19.txt UCLA-SZ_D1_V4_R1_UMI_001.fastq.sorted.fwd.bam UCLA-SZ_D1_V4_R1_UMI_001.fastq.sorted.rev.bam >UCLA-SZ_D1_V4_hapInfo.txt ../bam2hapInfo.pl /home/kunzhang/CpgMIP/MONOD/GBM_PC_CRC_V4_capture_regions_hg19.txt UCLA-SZ_H11_V4_R1_UMI_001.fastq.sorted.fwd.bam UCLA-SZ_H11_V4_R1_UMI_001.fastq.sorted.rev.bam >UCLA-SZ_H11_V4_hapInfo.txt
Sample | # on-target reads | # unique reads | Clonal rate |
BE2C_V4 | 862,623 | 442,154 | 48.7% |
BXPC3_V4 | 1,232,955 | 709,284 | 42.5% |
PANC1_V4 | 763,920 | 268,728 | 64.8% |
T98G_V4 | 1,107,621 | 305,041 | 72.5% |
U87MG_V4 | 1,011,039 | 333,677 | 67.0% |
UCLA-SZ_B1_V4 | 713,916 | 396,377 | 44.5% |
UCLA-SZ_D1_V4 | 472,779 | 86,148 | 81.8% |
UCLA-SZ_H11_V4 | 583,033 | 186,801 | 68.0% |
- Check the capture efficiencies for all probes using bam2probeEfficiency.pl: Batch processing script: get_probeEfficiency_all_files.sh; Merge all files into a matrix:
./get_probeEfficiency_matrix.pl > 131227_MiSeq_V4_probeEfficiency_matrix_UMI_20Mar14.txt
Some probes seemed to have very different capture efficiencies between cancer and whole blood. This could be due to the fact that two probes were synthesized for H1/H2 that having CpG sites and that the "C" probes have higher annealing efficiencies. I plotted the T-stat (between cancer and whole blood) against the number of CpG sites within the capture arms H1/H2, or within the insert.
File:MONOD V1 Capture bias vs CG H1H2.pngFile:MONOD V1 Capture bias vs CG insert.png
- TO BE CONTINUED.