Dinh/Dinh 2013/NOTES/2013-12-13

From ZhangLabWiki
Jump to navigation Jump to search

N37+blood probes design v2[edit]

Summary of workflow[edit]

  • We obtained methylation profiles for 10 tissues data from one individual and 1 blood data from three different individuals. Normalization was performed using BSmooth on the methylation values. The coefficient of variance (CV) for each site that was covered by at least 1X in at least 9 tissues was calculated with and without blood tissue. This resulted in measurements for ~88% of CpG sites (24.9 million sites). Another normalization was performed by using 100 bp non-overlapping bins and obtaining the average coefficient of variance across the CpG sites in that bin. We performed segmentation using a three-state Hidden Markov Model, with the genome divided such that consecutive bins were no more than 2 Mbp apart from each other.
  • We were interested in the high CV regions, and to identify those regions with high confident, we performed segmentation on the randomly shuffled CV values. We used an FDR cutoff of 10% and removed bins that has less than 99% state probability or less than 0.212 coefficient of variance. Consecutive bins with the same states and within 2 Kbp of each other were merged into regions. This resulted in totally 1.06 Gbp of low CV, 1.29 Gb of medium CV, and 0.219 high CV regions. Thus, we segmented ~83% of the genome using information from 88% of CpG sites.
  • We performed segmentation with and without the blood tissue and obtained similar distributions of regions and genome coverage (~83%). We merged the high CV regions from the two segmentation results and obtained a set of 147,016 non-overlapping targets to design bisulfite padlock probes. We used a gapsize of 290 to 310 bp and maximum 2 CpG sites in a capture arm. We successfully designed 1,300,119 probes with 97.6% coverage of the targeted regions. To assess the binding affinity of the arms to their targets in silico, we used Bowtie2 (--very-sensitive-local, k=300 valid alignments), and identified uniquely paired arms (paired was defined as being mapped to the same chromosome within 10 kbp). Paired arms must map uniquely to their targets. This resulted in 1,140,796 uniquely paired probes with 95.2% coverage of the targeted regions.
  • We calculated the number of CpG sites and the coefficient of variance across the coverage normalized methylation values within each captured region. We found that probes with less than 3 CpG sites or less than 0.40 CV were depleted in DNaseI hypsensitivity sites, so we rejected those probes. We kept 397,399 probes 60% coverage of the targeted regions. The missing 40% targeted regions might be false positives. While our new probes covers between 42-61% of the previous DMR330K subsets, the DMR330K set as a whole covers only ~19% of the targeted regions. Additionally, the greatest target overlap was for Subset 6 (58%), and Subset 7 (61%) which were the under performing probes from DMR330K. The coverage of CpG islands, shores, and shelves were similarly low for DMR330K and the new probes. The new probes has higher coverage of super-enhancers (92%) although the DMR330K probes also cover super-enhancer regions (85%). However, in terms of base pairs coverage, the new probes has ~1.5 million CpGs in super-enhancers while DMR330K has ~0.5 million CpGs in super-enhancers, thus, new probes has three times as many super-enhancers CpGs as DMR330K. As expected, DMR330K has higher (85%) coverage at Hg19 genes. However, the new probes has higher (50% versus 32%) coverage of lincRNA transcripts. We designed an additional 13,830 padlock probes to cover the remaining genes and lincRNA transcripts.
  • The final probes set has 411,229 probes that is enriched for tissue variability, DHS, TFBS, super-enhancers, and covers 99.45% of genes and 98.64% of lincRNAs. Since 103,931 probes has CpG sites in the capture arm regions, the total number of probes will be 515,160 to account for two versions of the probes with CpG sites.

Perform smoothing on 1x BED files with BSmooth[edit]

  • Use ns=40, h=2000
  • remove sites where less than 9 samples as 1X coverage
  • calculate the CV values at each site.

Perform binning using 100 bp non-overlapping, continuous windows[edit]

  • NOTE: 30,956,951 100bp bins in the human genome, but 16,430,018 (or 53% has no CpGs)
  • Calculate the average CV across the sites within each bin.

Perform segmentation[edit]

  • Segment genome with < 2 Mbp as the gapsize for segments (do not perform HMM prediction over 2 Mbp gap!)
  • Segment N37 only and N37 + blood separately
  • N37 + blood
 fdr = 128680/1352218 = 9.5%
 HMM probability cutoff = 0.99
 CV cutoff = 0.212
  • N37 only
 fdr = 130157/1359506 = 9.6%
 HMM probability cutoff = 0.99
 CV cutoff = 0.212
  • Region merging: Stadler et al defined borders at the middle point between sites, but this might make the target regions for probes design larger & slower
  • Merge regions with less than 10kbp of each other (we can filter non-cg or low cg probes later), its important to not miss anything
  • I calculate the genome coverage of our previous segmentation results, which is ~ 66.5 %, so we had improved the coverage using 1X sites.
  • N37 + blood merged regions:
 Low = 1.06 Gb
 Medium = 1.29 Gb
 High = 0.219 Gb
 ===
 83.1 % genome (3.09 Gb)
  • N37 only merged regions:
 Low = 1.06 Gb
 Medium = 1.29 Gb
 High = 0.221 Gb
 ====
 83.2 % genome (3.09 Gb)
  • The coverage of 100bp bins containing at least 1 CpG was above 90% for all autosomes. The coverage was poor for chrX and chrY, probably because reads were unable to map uniquely.
 File:N37+bloodv2 hmm segmentation chromosomes.png
  • Combine and obtain non-overlapping regions from both sets using bedtools merge.
  • Total target regions: n = 147,016 ( 235.8 Mbp, 7% of genome )
  • 100 Random targeted regions:
 File:N37+bloodv2 DMRs high cv heatmap.png
 File:N37+bloodv2 DMRs high cv density.png

Perform probe design[edit]

  • ppDesigner_v1.1.mod
  • jobFileMaster.pl:
 primerMaxLen=30;
 primerMinLen=24;
 H1_plus_H2_Len = 54;
 maxCG = 2; 
 targetMinLen=290;
 targetMaxLen=310;
  • Number of probes designed: n = 1,300,119 (384.4 Mbp, 97.6% of target regions - 143,523/147,016)

Perform mapping with Bowtie2[edit]

  • --very-sensitive-local, k=300 (maximum valid alignments)
  • Number of unique paired arms within 10kbp of each other: n = 1,140,796 (337.3 Mbp, 95.2% of targets - 139,894/147,016)

Perform filtering on probes with low CG and low coefficient of variation[edit]

  • Filter using min 3 CG and min 0.40 CV to get Keep Probes:

File:N37 blood probes DHS filter-check.png

  • Number of probes left: n = 397,399 (118.0 Mbp, 60 % of targets - 88,212/147,016)
  • Although 40% of the targeted regions were removed they might have very low CG density or were false positive regions.

Statistics of Probes Set[edit]

  • Comparing DMR330K probes, N37 targeted probes, N37 + genes targets probes
  • Overlapping coverage with DMR330K probes
File:Dmr330 n37 overlaps.png
  • Regions coverage
File:Dmr330 n37 cgi superenhancers genes.png
  • Transcripts coverage
File:Dmr330 n37 lincRNA genes.png
  • CpG Coverage
DMR330K New New+genes
Number of probes 244,339 397,399 411,229
Unique CpG sites 1,062,786 3,134,545 3,202,975
Total CpGs 1,449,888 4,741,770 4,811,035
Unique Infinium CpGs 95,443 217,038 223,556
Probe per CpG 1.36 1.51 1.50
CpG per Probe 5.93 11.93 11.70
  • N37+Genes chromosomes
File:N37 probes chromosomes.png
Regions Total DMR330K % new probes % new probes + genes %
N37_targets 147016 27592 18.77% 88210 60.00% 88799 60.40%
Sub1 3350 3350 100.00% 1406 41.97% 1417 42.30%
Sub2 3387 3387 100.00% 1424 42.04% 1434 42.34%
Sub3 28874 28874 100.00% 12263 42.47% 12368 42.83%
Sub4 29259 29259 100.00% 12384 42.33% 12468 42.61%
Sub5 28680 28680 100.00% 12898 44.97% 13006 45.35%
Sub6 12017 12017 100.00% 7005 58.29% 7039 58.58%
Sub7 6427 6427 100.00% 3917 60.95% 3936 61.24%
Exp1-3 132345 132345 100.00% 57772 43.65% 60712 45.87%
Genes 24494 20844 85.10% 19784 80.77% 24360 99.45%
LincRNA 21537 6994 32.47% 10744 49.89% 21244 98.64%
Super Enhancers 58283 49301 84.59% 53636 92.03% 53935 92.54%
CGI 28691 411 1.43% 419 1.46% 420 1.46%
CGI shores 28691 403 1.40% 424 1.48% 425 1.48%
CGI shelves 28691 367 1.28% 402 1.40% 402 1.40%
Genes TSS200 24494 9881 40.34% 16255 66.36% 19616 80.08%
Genes TSS1500 24494 13892 56.72% 16286 66.49% 19090 77.94%
Genes Body 20747 17105 82.45% 17942 86.48% 18502 89.18%
Genes UTR5 21054 11918 56.61% 15863 75.34% 17468 82.97%
Genes UTR3 21164 4638 21.91% 4867 23.00% 5567 26.30%
Genes First Exon 24494 11162 45.57% 16475 67.26% 19033 77.70%
LincRNA TSS200 21537 1660 7.71% 4711 21.87% 11753 54.57%
LincRNA TSS1500 21537 2916 13.54% 5721 26.56% 11784 54.72%
LincRNA body 21527 5432 25.23% 9005 41.83% 13023 60.50%
LincRNA First Exon 21537 1805 8.38% 4698 21.81% 9679 44.94%