Alice:Whole Genome Bisulfite Sequencing Lab Notes: Difference between revisions
Jump to navigation
Jump to search
>Zsakura2 |
>Zsakura2 No edit summary |
||
Line 92: | Line 92: | ||
TO run: | TO run: | ||
/home/dinh/scripts/allBED2Matrix_Jan25.pl [list_bed] [min_sample] [min_depth] [windowSize >= 1] | /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== | ==DMR finding using BSmooth and T-stat== | ||
Line 288: | Line 296: | ||
| 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 | | 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 | ||
|} | |} | ||
Revision as of 01:57, 23 July 2015
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)
- The samples were sequenced in the following runs:
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 also give 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 and iPSC: 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
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:
- 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.
- method #2:
- 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
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 |