Alice:Whole Genome Bisulfite Sequencing Lab Notes: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Zsakura2
>Zsakura2
No edit summary
Line 1: Line 1:
*[[Alice:Whole_Genome_Bisulfite_Sequencing_Lab_Notes/WGBS|WGBS]]
*[[Alice:Whole_Genome_Bisulfite_Sequencing_Lab_Notes/methylation reporter|methylation reporter knock-in]]
=WGBS of mouse SCNT and iPSC (collaborate with Yang Xu lab)=
=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.
*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.

Revision as of 00:09, 13 January 2016

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: [| Link]
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

  • 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
  6. 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):
  1. 160 are in topological domains (Dixon et al. Nature. 2012)
  2. 47 are in UMR & LMR distal regulatory regions (Stadler et al. Nature. 2009)
  3. 1 is in super enhancer target site (Whyte et al. Cell. 2013)
  4. 8 are in typical enhancer target site, 2 in subset of sites with ChiP-seq density above 1 for Oct4, Sox2, Nanog only (Hnisz et al. Cell. 2013)
  5. 35 are in CpG shore, 13 in CpG shelf

Validation of DMRs identified

  • validate some DMRs identified above in all other mSCNT, miPSC lines
  • samples to be processed for bisulfite conversion:
    • SCNT-ESC (official ID): NB4 (SCNT4), B7 (SCNT5), B8 (SCNT6), B9 (SCNT7), B10 (SCNT8), B11 (SCNT9), B13 (SCNT11)
    • iPSC (official ID): miPS1 (miPSD1), miPS3 (miPSD3), miPS4 (miPSD8), 2eiPS2 (2eiPS2P5)
    • MEF(official ID): B6P1-1(), B6P1-2(), TB_MEF(), ZN_MEF()
  • use 500ng of gDNA as input

Bisulfite conversion

  • use EZ-96 DNA Methylation-Lightning™ MagPrep kit and protocol
    • Add 130 µl of Lightning Conversion Reagent to 20 µl of a DNA sample in a Conversion Plate. Mix the samples by pipetting up and down.
    • Seal the plate with the provided film. Transfer the Conversion Plate to a thermal cycler and perform the following steps:
  1. 98°C for 8 minutes
  2. 54°C for 60 minutes (then 4°C storage for up to 20 hours if needed)
    • Pre-heat a plate heating element to 55°C.
    • Add 600 µl of M-Binding Buffer and 10 µl of MagBinding Beads to each well of a Collection Plate.
    • Transfer the samples from the Conversion Plate into the Collection Plate containing the M-Binding Buffer and MagBinding Beads. Mix by pipetting up and down 3-6 times and/or vortexing at 1,300-1,500 rpm for 30 seconds
    • Let plate stand at room temperature for 5 minutes, then transfer plate to a magnetic stand for an additional 5 minutes or until beads pellet and supernatant is cleared. With the plate on the magnetic stand remove the supernatant and discard.
    • Remove the Collection Plate from the magnetic stand for this and each subsequent buffer addition. Add 400 µl of M-Wash Buffer to the beads. Re-suspend the beads by pipetting up and down or vortexing the plate at 1,300-1,500 rpm for 30 seconds. Replace the plate on the magnetic stand for 3 minutes or until beads pellet. Remove and discard supernatant.
    • Add 200 µl of L-Desulphonation Buffer to the beads. Re-suspend the beads by pipetting up and down or vortexing for 30 seconds. Let plate stand at room temperature (20-30°C) for 15-20 minutes. After the incubation, replace the plate on the magnetic stand for 3 minutes or until beads pellet. Remove and discard supernatant.
    • Add 400 µl of M-Wash Buffer to the beads. Re-suspend the beads by pipetting up and down or vortexing for 30 seconds. Replace the plate on the magnetic stand for 3 minutes or until beads pellet. Remove and discard supernatant. Repeat this wash step.
    • Transfer the plate to a heating element at 55°C for ~15 minutes to dry the beads (be careful don't over-dry) and remove residual M-Wash Buffer.
    • Add 30 µl of M-Elution Buffer directly to the dried beads and pipette or vortex for 30 seconds to re-suspend. Heat the elution at 55°C for 4 minutes then transfer the plate to the magnetic stand for 1 minute or until beads pellet. Remove the supernatant and transfer to a clean Elution Plate.
    • use ~10ng of converted DNA for subsequent PCR

PCR

  • q-PCR to check the bisulfite converted DNA
  • Perform in 50ul reaction with 1ul of bis-cvt DNA
Components volume (ul)
bis-cvt DNA 1.00
10uM Me-Enhancer primer F 1.00
10uM Me-Enhancer primer R 1.00
2x KAPA fast MM 25.00
H2O 21.00
Total 50.00

Program

  • 98C 3min --> (98C 10sec -> 60C 30sec -> 72C 45sec) x 15 cycles --> 72C 3min --> 15C Hold
  • The PCR didn't work, so I quickly stopped and changed annealing temperature
  • 98C 3min --> (98C 10sec -> 58C 30sec -> 72C 45sec) x 12 cycles --> 72C 3min --> 15C Hold
  • PAGE gel results:

Sequencing library preparation

  • In order to be consistent with WGBS data, I will use IQ MM for PCR instead of KAPA
  • The primer design of first round and second round PCR can be found here:
    • briefly, we first design primers that span part (or entire) region of the DMR identified. Then we add part of the Y-adaptor onto the end of both primers.
    • It is the complementary part of the Y-adaptor that will be added onto the end of both primers, then we can use PCR_F and PCR_N2R for 2nd round PCR
    • add in front of all Forward primer: 5' GCGAGAAGGCTAG 3'
    • add in front of all Reverse primer: 5' GCGAGAAGGCTAG 3'
  • sample info:
' sample name N2 barcode
1 NB1 PCR_N2R10
2 SCNTB7 PCR_N2R11
3 SCNTB8 PCR_N2R12
4 SCNTB9 PCR_N2R13
5 SCNTB10 PCR_N2R14
6 SCNTB11 PCR_N2R15
7 SCNTB13 PCR_N2R16
8 B6P1-1 PCR_N2R17
9 B6P1-2 PCR_N2R18
10 miPSD1 PCR_N2R19
11 miPSD3 PCR_N2R20
12 miPSD8 PCR_N2R21
13 ZN_MEF PCR_N2R22
14 TB_MEF PCR_N2R23
15 2eiPS2P5 PCR_N2R24

PCR

Components 1rxn 9 rxn mix (9 primer pairs)
bis-cvt DNA 1.00 0.00
10uM Me-Enhancer primer F 1.00 0.00
10uM Me-Enhancer primer R 1.00 0.00
2x iQ MM (Biorad) 25.00 188
H2O 21.00 160
Total 50.00 375

PCR program

  • after two failed runs with 98C as denaturing temperature, final program used for big batch processing:

95C 2min --> (95C 15sec -> 58C 20sec -> 72C 30sec) x 26-28 cycles --> 72C 3min --> 15C Hold

  • combine the PCR amplicon of different targets into the same tube for the same sample, then perform Qiaquick column purification and elute in 33ul H2O
  • so now there are 15 samples (each with all 9 targets combined into one tube) to run 2nd PCR

2nd PCR with Illumina primers

Components 1rxn 15 rxn mix
bis-cvt DNA targets pooled 10.00 -
10uM PCR F 1.00 16.00
10uM PCR R_N2 barcoded 1.00 -
2x iQ MM (Biorad) 25.00 380
H2O 21.00 315
Total 60.00

PCR program 95C 2min --> (95C 15sec -> 60C 30sec -> 72C 30sec) x 32 cycles --> 72C 3min --> 15C Hold

  • Qiaquick column purification and elute in 42ul H2O each (15 total)
  • loaded 2ul of purified samples for PAGE gel validation
  • PAGE gel results: samples 1,2,4,9,13 had high yield, while the rest of them had very concentrations
  • so I did the pooling based on eye-balling the gel results:
    • 3ul of the samples with high yield (sample 1, 2, 4, 9, 13)
    • 12ul of the samples with relative low yield (sample 10, 11, 14)
    • 15ul of the samples with low yield (sample 3, 5, 6, 7, 8, 12, 15)
  • after pooling, used the total volume ~156ul for gel size selection (~120-220bp bands)
  • ethanol precipitation and elute in total 60ul H2O (blue color due to ~3ul of glycoblue per tube, total of 18ul, was added during first step)
  • gel validation
  • samples were sequenced on Miseq and data is stored under the folder: /home/kunzhang/seqStore/151024_HL017