WGBS of mouse SCNT and iPSC (collaborate with Yang Xu lab)
- The goal of this project is to examine the methylation pattern of mouse SCNT and iPSC against mESC using whole genome bisulfite sequencing, and to characterize epigenetic signatures that resulted from different stem cell reprogramming methods.
- detailed experimental procedures can be found here: [[1]]
sample ID
|
N2 Index used
|
cell type
|
progenitor line
|
progenitor/donor mouse strain
|
2A4F1 |
N2 index 1 |
iPSC |
N/A (used TB_B6MEF in exome study) |
C57BL/6J
|
2A4F33 |
N2 index 2 |
iPSC |
N/A (used TB_B6MEF in exome study) |
C57BL/6J
|
miPS B3 |
N2 index 4 |
iPSC |
ZN_B6MEF |
C57BL/6J
|
1E12P20 |
N2 index 5 |
iPSC |
TB_B6MEF |
C57BL/6J
|
nt ES P7 C |
N2 index 6 |
SCNT |
B6MEF P0-3 |
C57BL/6J
|
nt ES P8 B |
N2 index 7 |
SCNT |
B6MEF P0-3 |
C57BL/6J
|
NB 3 |
N2 index 9 |
SCNT |
B6-P1-1 |
C57BL/6J
|
SCNT B12 |
N2 index 10 |
SCNT |
B6-P1-2 |
C57BL/6J
|
B6 ES |
N2 index 11 |
ES |
- |
C57BL/6J
|
129 ES |
N2 index 12 |
ES |
- |
-
|
Sequencing statistics
- Data Source
- The samples were sequenced in the following runs:
- Hiseq 130104_SN1001
- GAIIx 121226_HL140 (lane 1 only)
Sample ID
|
N2 Index
|
total reads
|
mapped reads
|
mapping rate
|
total # of CpG sites covered (coverage >=5)
|
mean coverage (sites with coverage >=5)
|
mean methylation level (sites with coverage >=5)
|
2A4F1 |
1 |
375578848 |
283699999 |
76% |
17510624 |
12.5 |
0.76
|
2A4F33 |
2 |
360720634 |
253128940 |
70% |
16863956 |
11.3 |
0.74
|
mipsB3 |
4 |
342349949 |
211416203 |
62% |
15215962 |
9.8 |
0.71
|
1E12P20 |
5 |
328678606 |
227406985 |
69% |
15865114 |
10.6 |
0.73
|
nt ES P7C |
6 |
315952780 |
205976115 |
65% |
15025471 |
9.6 |
0.57
|
nt ES P8B |
7 |
461470678 |
333879642 |
72% |
17852991 |
12.1 |
0.61
|
NB3 |
9 |
297479356 |
205398733 |
69% |
14210572 |
9.3 |
0.79
|
SCNTB12 |
10 |
322193316 |
228779506 |
71% |
16016176 |
11.1 |
0.73
|
B6ES |
11 |
310311719 |
214278139 |
69% |
15464407 |
10.5 |
0.74
|
129ES |
12 |
271922039 |
163956125 |
60% |
7433945 |
7.7 |
0.64
|
Pre-processing scripts
- data and processed files can be found on meangenemachine:/home/zhl002/ZL_1TB
- Dinh also has processed files on genome-miner:/media/LTS_15T/Dinh_LTS/WGBS/SCNT_mm9/Soap
Mapping and calculate MethylFreq
perform mapping using SOAP
filtering and processing of alignments using SAMtools and Picard
then calculate methylation frequency
make BED format files
/home/dinh/scripts/methylFreq2BED.pl
TO run:
/home/dinh/scripts/methylFreq2BED.pl [sample_name] [min_depth] < [methylFreq file] > Sample_name.bed.txt
Quality control: calculate forward and reverse correlation
/home/dinh/scripts/frMethylCorr.pl
TO run:
/home/dinh/scripts/frMethylCorr.pl [min_depth] < [methylFreq file] > Sample_name.corr
generate methylation matrix
/home/dinh/scripts/allBED2Matrix_Jan25.pl
TO run:
/home/dinh/scripts/allBED2Matrix_Jan25.pl [list_bed] [min_sample] [min_depth] [windowSize >= 1]
additional public datasets
- Previous studies will be used in my analysis to cross check my results:
- GSE57179: Whole-genome single-base resolution methyl-C sequencing "Abnormalities in human pluripotent cells due to reprogramming mechanisms"
- GSE41923: Whole-genome bisulfite sequencing of two distinct interconvertible DNA methylomes of mouse embryonic stem cells
- GSM1027571 DNA_Methylation_serum_LIF_E14 link
- GSM1027572 DNA_Methylation_E14_adapted_2i link
- GSM1127953 DNA_Methylation_E14_serum_LIF_replica link
- GSM1127954 DNA_Methylation_E14_adapted_2i_replica link
- Related dataset not used in my analysis:
- GSE53060: Illumina HumanMethylation450 BeadChip array "Abnormalities in human pluripotent cells due to reprogramming mechanisms"
- GSE61461: Illumina HumanMethylation450 BeadChip array "Comparable frequencies of coding mutations and loss of imprinting in human pluripotent cells derived by nuclear transfer and defined factors"
- GSE53094: RNA-seq "Abnormalities in human pluripotent cells due to reprogramming mechanisms"
- GSE61656: RNA-seq "Comparable frequencies of coding mutations and loss of imprinting in human pluripotent cells derived by nuclear transfer and defined factors"
DMR finding using BSmooth and T-stat
perform smoothing on the data using BSmooth
BSmooth R script
bsseq package usage and description: Link
UserGuide
TO run:
script_dir="/media/1TB_storeB/ZL_meangenemachine/Habibi_smoothing"
bed="/media/1TB_storeB/ZL_meangenemachine/BEDfiles"
for c in chr#
do -> for f in list_of_file_names
do -> grep $c [space]$bed/$f.bed.txt > $f.$c
for ns in number_CpG_sites
do -> for h in min_window_size
do -> sed "s/TESTNS/$ns/g" runSmoothing.R | sed "s/TESTH/$h/g" | sed "s/CHROM/$c/g" > run.R
output files:
file_name.chr#
chr#.ns#_h#.rda
perform T-statistics on smoothed data
Tstat R script
TO run:
for f in *rda
do -> sed "s/RDA/$f/g" runTstat.R > Run.R
output files:
chr#.ns#_h#.rda_comparison_coveredsites
chr#.ns#_h#.rda_comparison_DRMs.txt
chr#.ns#_h#.rda_comparison_plotDMRs.pdf
highlight:
group1 and group2 contain the sample names of the two groups being compared
estimate.var describes which sample are being used to estimate variability, and choices are:
gr# (e.g normal sample when compared again cancer data)
same(assume same variability in each group)
paired(do a paired T-test)
In my analysis, I used gr2 to estimate variability, which are the ES lines
plotting the DMRs
T-stat will generate txt file that identified potential DMRs between two group of samples
Dinh modified the methyl matrix perl script which allows it to use chr_interval as input and generate matrix only within the interval specified
on meangenemachine: /home/zhl002/ZL_1TB/allBED2Matrix_Jan222014.pl
TO run: perl allBED2Matrix_Jan222014.pl [list_bed] [min_sample] [min_depth] [windowSize >= 1] chr#:#-#
output: MethylMatrix_bed_list.min#.md#.win#
Then use plotting DMR R script to generate the plot
use vim to change the first line of the script to use the output matrix that we just generated from last step
results
used previously published WGBS on ES cell cultured under 2i and serum conditions Link
do BSmooth on all the data: ns=20 and h=2000 (20 CpG sites in the window that range from minimum =2000 and maximum =10,000)
then run Tstat to do the comparison:
SCNT-serum (indx 9 and 10) against ES-serum (indx11 and Habibi_ES_serum) (ES serum is group2 and also used to estimate variability)
iPSC-serum (indx 1, 2, 4, 5) against ES-serum (indx11 and Habibi_ES_serum) (ES serum is group2 and also used to estimate variability)
SCNT-2i (indx 6 and 7) against ES-serum (Habibi_ES_2i_adapted and Habibi_ES_2i) (ES 2i is group2 and also used to estimate variability)
Among the candidate DMRs generated from the above steps, we were able to detect true culture-specific DMRs, where there are differentially methylated regions between 2i and serum conditions. However, because this method relies on biological replicates to do smoothing and calculate T-statistics, so we have decided to use other published method to find DMRs.
DMR finding using MOAB
- detailed method and pre-liminary analysis of the pipeline can also be found here: [Dinh's wiki page]
DMR plots and method comparison
All 4 ESCs vs IPS comparisons overlap at chr19:53899387-53899586 (Rbm20 gene)
File:Methylgraph chr19 53898966-53899834.pdf
I plotted one of the locus where all 4 ESCs vs SCNT comparisons overlap at 3 locus on chr2
chr2:152511745-152513298
File:Methylgraph chr2 152511745-152513298.pdf
SCNT 2i versus SCNT serum:
large # of DMRs were observed
previous study has reported methylome of male ESCs in serum parallels postimplantation blastocyst cells, and 2i stalls ESCs in a hypomethylated, ICM-like state
increase the minCredibleDiff to 0.5 and reduce the amount of DMRs for validation purpose
Discussion
- we focused on the DMRs generated by method 2 and method 3. Although method 3 has the lowest FDR, it might be too stringent because the candidate DMRs seem to be biased towards samples with higher coverage, therefore method 2 would be more appropriate
- by down-sampling the published ESC data, we are expected to observe less of strain-specific DMRs
- we will repeat the comparison between individual ESC and iPSC/SCNT
DMR finding using DSS
- Yupeng (Joe Ecker group student) suggested to try DSS for DMR finding, as it was used by the Ecker group on the human SCNT paper. (DSS is an R library performing differntial analysis for count-based sequencing data. It detectes differentially expressed genes (DEGs) from RNA-seq, and differentially methylated loci or regions (DML/DMRs) from bisulfite sequencing (BS-seq). The core of DSS is a new dispersion shrinkage method for estimating the dispersion parameter from Gamma-Poisson or Beta-Binomial distributions. LINK)
- Shicheng helped with modifying and writing script
1. run DSS (defined fdr_cutoff=0.05) between 2 groups: group1=miPS_B3, miPS_1E12P20; group2=SCNT_NB3, SCNT_B12
2. DMS (differentially methylated sites) output from DSS:
chr;
pos;
mu1 (mean methylation level of group1);
mu2 (mean methylation level of group2);
diff (difference between the means);
diff.se (diff.se - standard error of the difference);
stat (stat - test statistic);
phi1 (phi of group1);
phi2 (phi of group2);
pval (p-value);
fdr (false discovery rate);
postprob.overThreshold (reprsenting the posterior probability of the difference in methylation greater than delta)
3. download annotated CpG island-related (CpGi, CpG shore, CpG shelf) and repeat-related (LINE, SINE, LTR, SimpleRepeat) files from UCSC
4. use bedtools to intersect DSS output DMS with annotated functional regions
5. record output # of hits
6. perform random permutation (n=100,000) to detect the significance of the enrichment at those regions
(1) generate a tmp file with # of regions (each have random 2bp overlap within the mouse genome) matching with # of DMR calculated by DSS
(2) intersect the tmp file with each functional region of interest (e.g. CpGi, CpG shore, etc...)
(3) do a word count of the # lines that would show up in the functional regions
(4) repeat for 100,000 times
(5) assuming data from random permutation is normally distributed, then calculate the p-value for the level of enrichment in our DMR set
Results
- total # of DMS between human SCNT and iPSC (hg19, Ma et al): 111
- Here are the # of DMRs we found that overlapped with different functional regions within the human genome (annotation downloaded from UCSC table browser):
- CpG island: 7 (2bp DMS p-value: <10^-6)
- CpG shore: 12 (2bp DMS p-value: 0.00063, fold of enrichment: 3.02 )
- CpG shelf: 10 (2bp DMS p-value: 0.00524
- total # of DMS between mouse SCNT (2 lines cultured with serum) and 4 iPSC lines: 489 (0.05 FDR)
- Here are the # of DMRs we found that overlapped with different functional regions within the mouse genome (annotation downloaded from UCSC table browser):
- CpG island: (2bp DMS p-value: )
- CpG shore: 28 (2bp DMS p-value: 5.214671e-05; fold of enrichment:)
- CpG shelf: (2bp DMS P-value: )
- DNase I hypersensitive sites (UCSC ENCODE data from UW): ~88 (2bp DMS, p-value: <10^-6)
- Enhancer for Oct4, Sox2, Nanog, Med1, H3K27Ac, Klf4, and Esrrb in mESC: 10 (2bp DMS, p-value: 0.00053) (reference: Master Transcription Factors and Mediator Establish Super-Enhancers at Key Cell Identity Genes. Whyte et al. Cell. 2013)
- superenhancer and typical enhancer ChiP-Seq data: (reference: Super-Enhancers in the Control of Cell Identity and Disease, doi:10.1016/j.cell.2013.09.053)
- total # of DMR we identified overlapped with topological domains that were previously published by Ren lab (source: Dixon JR et al. Topological domains in mammalian genomes identified by analysis of chromatin interactions. Nature. 2012;485:376-380)
- J1 mESC data (HindIII_combined,mm9)
- total # of DMS overlapped with topological domains: 457 (2bp DMS, p-value: 0.013291)
Discussion
- We have found significantly level of differences between iPSC and SCNT at CpG shore and shelf
- Although we have confirmed some DMS by eye balling, however, we can't rule out the possibility of epigenetic memory from progenitor, sample-specific phenomenon, etc.
- we will perform the same analysis on human SCNT, iPSC sample sets to determine whether we can find the consistent pattern.
- the discrepancy does not mean what we found were not biologically meaningful. hESCs and mESCs are very different in terms of DNA methylation. For example, abundant non-CG methylation does not show up (or appear at extremely low level) in mESCs. Thus, discrepancy in methylation patterns between hESCs and mESCs is expected. We will keep looking at overlap of these DMRs with DHSs (from Dnase-seq), putative enhancers (based on chromatin data) and transcription factor binding sites (TF ChIP-seq data).
correlation between CpG island shore methylation and gene expression levels
- download the gene annotation file from UCSC
- use "bedtools closest" command to find the closest CpG island shore location of a given gene
- trim down the output file using following criteria:
- genes on + strand with end of CpG shore- beginning of gene < 0;
- genes on - strand with end of gene - beginning of CpG shore > 0
- now we have created all the known CpG island shore location for a list of genes
- compile a list of genes that are differentially expressed between SCNT-ESC and iPSC or iPSC and ESC (look for literature)
- find the methylation level of the CpG shore for these genes
- how to calculate the p-value?
- split CpG island shore methylation into 2 groups: differentially methylated and non-differential methylation
- find genes that are differentially expressed in the above two groups
- chi-square test for relationship between CpG shore methylation and gene expression
- Discussion: this method isn't going to work because we have only observed enrichment of DMR at CpG shore, it is not a very strong signal such that high # of DMRs are all located at the shore, so the observations we can fit into the above groups would be very small. So we should try to find another way to look for correlation between the groups.
- find DMS that is located at CpG island shore using bedtools (found 12 out of 111 total DMS)
- extend the CpG shore region we have identified from previous step by 2k in both direction (total length of the region is now 6k)
- extract the methylation frequency at these locations
- find the gene expression level of these genes
- do correlation test between the methylation and gene expression level
- the genes below are the three candidates identified by above method, but the results weren't informative to make any conclusion
GENE*Transcript
|
HSPA2*NM_021979
|
ZNF212*NM_012256
|
SIM1*NM_005068
|
Results
- Some thoughts on the results we have so far:
- comparable global DNA methylation level between mES, miPSC, and mSCNT ESC (cluster dendrogram, histogram/density plot)
- DNA methylation aberration in both iPS & SCNT or one cell type only compare to ESC (heat map, DMR plot)
- differentially methylated sites/regions is enriched in regulatory elements such as CpG shore, enhancers, and DNaseI hypersensitive sites (Random sampling simulation to determine the distribution of the expected number of genes that are targets of the indicated regulatory element and have aberrant DNA methylation CpGs overlapping the site; histogram with box plot)
- these regulatory elements are also crucial for dynamics genome-wide de novo methylation/demethylation of serum and 2i state of SCNT ESCs
- DMRs identified in (2) have the following features (total # is 173):
- 160 are in topological domains (Dixon et al. Nature. 2012)
- 47 are in UMR & LMR distal regulatory regions (Stadler et al. Nature. 2009)
- 1 is in super enhancer target site (Whyte et al. Cell. 2013)
- 2 are in typical enhancer target site, subset of sites with ChiP-seq density above 1 for Oct4, Sox2, Nanog only (Hnisz et al. Cell. 2013)
- 35 are in CpG shore, 13 in CpG shelf