Alice:Whole Genome Bisulfite Sequencing Lab Notes

From ZhangLabWiki
Revision as of 02:16, 29 July 2015 by >Zsakura2 (→‎DMR finding using DSS)
Jump to navigation Jump to search

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

  • GSE57179: Whole-genome single-base resolution methyl-C sequencing "Abnormalities in human pluripotent cells due to reprogramming mechanisms"
  • 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. add 20 bp to each side of the DMS to convert into 40bp DMR bed file

4. download annotated CpG island-related (CpGi, CpG shore, CpG shelf) and repeat-related (LINE, SINE, LTR, SimpleRepeat) files from UCSC

5. use bedtools to intersect DSS output DMR with annotated functional regions

6. record output # of hits

7. perform random permutation (n=1,000,000) to detect the significance of the enrichment at those regions
     (1) generate a tmp file with # of regions (each have random 40bp 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 1,000,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 mouse SCNT (2 lines cultured with serum) and 4 iPSC lines: 593
  • Here are the # of DMRs we found that overlapped with different functional regions within the mouse genome (annotation downloaded from UCSC table browser):
    • LINE: 52
    • SINE: 64
    • LTR: 52
    • CpG island: 6 (2bp DMS p-value: 0.01381)
    • CpG shore: 30 (2bp DMS p-value: 0.00012, fold of enrichment: 2.082444)
    • CpG shelf: 21 (2bp DMS P-value: 0.0407)
  • 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: )
    1. 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)
      1. of DMS overlapped with topological domains: 457 (2bp DMS, p-value: 0.013291)
    1. of DMS overlapped with DNase I hypersensitive sites (downloaded from UCSC): 88 (2bp DMS, p-value: <10^-6)

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

  • method #1:
  1. download the gene annotation file from UCSC
  2. use "bedtools closest" command to find the closest CpG island shore location of a given gene
  3. trim down the output file using following criteria:
    1. genes on + strand with end of CpG shore- beginning of gene < 0;
    2. genes on - strand with end of gene - beginning of CpG shore > 0
  4. now we have created all the known CpG island shore location for a list of genes
  5. compile a list of genes that are differentially expressed between SCNT-ESC and iPSC or iPSC and ESC (look for literature)
  6. 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.
  • method #2:
  1. find DMS that is located at CpG island shore using bedtools (found 12 out of 111 total DMS)
  2. extend the CpG shore region we have identified from previous step by 2k in both direction (total length of the region is now 6k)
  3. extract the methylation frequency at these locations
  4. find the gene expression level of these genes
  5. do correlation test between the methylation and gene expression level
GENE*Transcript GSM1282330 GSM1378026 GSM1282324 GSM1378018 GSM1282325 GSM1378019 GSM1282326 GSM1378020 GSM1282327 GSM1378021 GSM1282328 GSM1378022 GSM1282329 GSM1378023 GSM1378024 GSM1378025 GSM1282320 GSM1378014 GSM1282321 GSM1378015 GSM1282322 GSM1378016 GSM1282323 GSM1378017
HSPA2*NM_021979 1131.679811 686.4630589 987.4968903 868.651013 820.1849877 751.1453967 492.3872812 535.3525862 203.3407777 372.4418482 144.4719756 402.1298088 279.0146231 472.30402 105.9115375 428.0794859 638.228473 563.2880447 614.8179393 560.9098977 712.4546507 980.9216258 1068.914232 704.9663454
ZNF212*NM_012256 159.2543968 160.4801879 407.4597347 406.3765897 364.3367059 387.8361933 356.002106 271.1569849 384.5609968 267.3194562 403.13993 363.9311311 292.770929 389.4377314 253.5111249 357.5135382 333.7341851 279.2841758 371.8290353 364.5340547 380.132046 443.0885574 423.0458649 396.8179543
SIM1*NM_005068 544.398545 609.9363563 4.728399463 0.649516615 12.91624204 4.531632842 7.788652011 2.744830325 14.34773727 8.141612683 5.884757342 2.586335468 4.782406994 -1.660132734 12.37626581 11.3184976 4.764579097 2.816229225 8.656889871 -1.534563558 6.597318359 3.003275712 7.67742571 0.101714059