Dinh 2011/NOTES/2011-9-12
Jump to navigation
Jump to search
5fC DIP-seq[edit]
HL102 Dipseq-Input-Control-Sep2 : s1, s2 Dipseq-Input-KD-Sep2 : s3, s4 Dipseq-IP-Control-Sep2 : s5, s6 Dipseq-IP-KD-Sep2 : s7, s8
Mapping[edit]
bowtie -n 2 --best --strata -l 28 -k 1 -m 1 -p 12 --phred64-quals --un $f.unmappable.txt /media/1TB_store1/Rui_Scratch/iGenome/Mus_musculus/UCSC/mm9/Sequence/BowtieIndex/genome $f $f.bowtie.map 2> $f.log &
**Each lane was mapped separately, data was then concatenated.
Mapping rates
' | Unique Hit | No Hit | >1 Hits | TOTAL | %Unique Hit | %No Hit | %>1 Hits |
s1 | 24626804 | 1084945 | 7643512 | 33355261 | 73.83% | 3.25% | 22.92% |
s2 | 25108666 | 1082767 | 7818626 | 34010059 | 73.83% | 3.18% | 22.99% |
s3 | 25759355 | 1188545 | 8379585 | 35327485 | 72.92% | 3.36% | 23.72% |
s4 | 21717915 | 1067567 | 7093600 | 29879082 | 72.69% | 3.57% | 23.74% |
s5 | 19239845 | 1244633 | 7924574 | 28409052 | 67.72% | 4.38% | 27.89% |
s6 | 23397892 | 1180981 | 9638993 | 34217866 | 68.38% | 3.45% | 28.17% |
s7 | 23974239 | 1228065 | 9301572 | 34503876 | 69.48% | 3.56% | 26.96% |
s8 | 22722261 | 1188036 | 8854346 | 32764643 | 69.35% | 3.63% | 27.02% |
- Note that Input has more hits than IP.
cat s_1_sequence.txt.gz.fastq.bowtie.map s_2_sequence.txt.gz.fastq.bowtie.map > 5fC_dIP-seq_Input_Ctrl.bowtie.map & cat s_3_sequence.txt.gz.fastq.bowtie.map s_4_sequence.txt.gz.fastq.bowtie.map > 5fC_dIP-seq_Input_KD.bowtie.map & cat s_5_sequence.txt.gz.fastq.bowtie.map s_6_sequence.txt.gz.fastq.bowtie.map > 5fC_dIP-seq_IP_Ctrl.bowtie.map & cat s_7_sequence.txt.gz.fastq.bowtie.map s_8_sequence.txt.gz.fastq.bowtie.map > 5fC_dIP-seq_IP_KD.bowtie.map &
MACS[edit]
Downloads and installed MACS debian package onto meangenemachine http://liulab.dfci.harvard.edu/MACS/ Followed instructions to install PeakSplitter_v1.0 (copied to /usr/local/bin) (**Results from peak splitter was not used because no p-value given for each peak)
Control[edit]
macs14 -t 5fC_dIP-seq_IP_Ctrl.bowtie.map -c 5fC_dIP-seq_Input_Ctrl.bowtie.map -f BOWTIE -g mm -n Control-Dipseq -w --call-subpeaks > macs14_ctrol &
All peaks called with p-value (based on Poisson distribution with local lambda) and FDR (%) File:Control-Dipseq peaks.xls Negative peaks were called by switching input and ip File:Control-Dipseq negative peaks.xls
Total peaks = 512; Total negative peaks = 25 (5.08% FDR for this method) - I should down-sample hits so that input and ip have equal number of hits.
File for viewing peaks (FDR<5%) on genome browser: File:Control-Dipseq peaks.BED.txt
Knockdown[edit]
macs14 -t 5fC_dIP-seq_IP_KD.bowtie.map -c 5fC_dIP-seq_Input_KD.bowtie.map -f BOWTIE -g mm -n KD-Dipseq -w --call-subpeaks > macs14_kd &
All peaks called with p-value (based on Poisson distribution with local lambda) and FDR (%) File:KD-Dipseq peaks.xls Negative peaks were called by switching input and ip File:KD-Dipseq negative peaks.xls
Total peaks = 317; Total negative peaks = 23 (7.57% FDR for this method)
File for viewing peaks (FDR<5%) on genome browser: File:KD-Dipseq peaks.BED.txt
Enrichment in genomic regions[edit]
(A) Gene regions were taken from iGenome data which Rui downloaded: /media/1TB_store1/Rui_Scratch/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes File:RefFlat.txt.gz File was parsed into five regions files using this script : File:ParseRefFlat.txt mm9_ucsc_genes.3utr.bed - regions from coding start end to transcription end N = 28,108 sort -u | wc -l : 24,052 mm9_ucsc_genes.5utr.bed - regions from transcription start to coding start N = 28,108 sort -u | wc -l : 24,143 mm9_ucsc_genes.tss.bed - regions from -800 to +200 bp of transcription start N = 28,108 sort -u | wc -l : 23,652 mm9_ucsc_genes.exons.bed - coding regions N = 267,543 sort -u | wc -l : 208,522 mm9_ucsc_genes.introns.bed - non-coding regions (between every exons) N = 239,435 sort -u | wc -l : 186,730
**Performed sort -u on all files.
(B) CpGIslands-regions.bed : downloaded from UCSC Table Browser for mm9 (C) mouseESp300binding_LICR_regions.bed : downloaded from UCSC Table Browser for mm9 (D) REST-binding_sites.bed (Repressor Element 1 Silencing Transcription Factor binding sites) : downloaded from UCSC Table Browser for mm9
Criteria for mapping to genomic regions[edit]
- Count frequency of overlaps with a genomic region:
For peaks: Total # bp overlapping with peaks For randomized regions: Total # bp overlapping with random regions
** Note Pastor et al required for larger regions that entire HERGs/MERGs fall within a larger genomic region.
Criteria for randomized regions[edit]
- Randomized regions conserve # peaks, chromosome distribution, and peak widths distribution.
- Performed 1000 randomizations (Imputations).
- Script: File:ImputeRegionEnrichment.txt
Results[edit]
Control sample
Total bp = 400,486bp
Region_File | Test_counts | Imputed_ave | Imputed_std | Empirical_p (greater than) | Empirical_p (less than) |
CpGIslands-regions.bed | 5632 | 1754.636 | 973.5852 | 0 | 1 |
mm9_ucsc_genes.3utr.bed | 6782 | 12116.25 | 3978.4707 | 0.929 | 0.071 |
mm9_ucsc_genes.5utr.bed | 22232 | 16190.186 | 4545.6284 | 0.087 | 0.913 |
mm9_ucsc_genes.exons.bed | 3856 | 9942.256 | 2401.9812 | 0.999 | 0.001 |
mm9_ucsc_genes.introns.bed | 143652 | 149212.281 | 11603.9548 | 0.688 | 0.312 |
mm9_ucsc_genes.tss.bed | 4479 | 3729.626 | 1575.1404 | 0.296 | 0.704 |
mouseESp300binding_LICR_regions.bed | 6007 | 3879.494 | 1621.2735 | 0.112 | 0.888 |
REST-binding_sites.bed | 37 | 28.687 | 57.2239 | 0.238 | 0.762 |
- Control is enriched in CpG Islands and 5'UTR.
- Control is depleted in 3'UTR and exons.
Knock down
Total bp = 220,304bp
Region_File | Test_counts | Imputed_ave | Imputed_std | Empirical_p (greater than) | Empirical_p (less than) |
CpGIslands-regions.bed | 121 | 950.83 | 760.1284 | 0.876 | 0.124 |
mm9_ucsc_genes.3utr.bed | 4465 | 6457.113 | 2854.0275 | 0.742 | 0.258 |
mm9_ucsc_genes.5utr.bed | 5055 | 8971.008 | 3276.3609 | 0.891 | 0.109 |
mm9_ucsc_genes.exons.bed | 436 | 5460.863 | 1840.9045 | 1 | 0 |
mm9_ucsc_genes.introns.bed | 86823 | 81022.857 | 8769.8664 | 0.248 | 0.752 |
mm9_ucsc_genes.tss.bed | 0 | 2064.625 | 1219.295 | 1 | 0 |
mouseESp300binding_LICR_regions.bed | 1278 | 2084.766 | 1255.2735 | 0.707 | 0.293 |
REST-binding_sites.bed | 26 | 17.38 | 46.7585 | 0.169 | 0.831 |
- Knockdown is depleted at exons and tss.
Reads balancing[edit]
head -42637737 5fC_dIP-seq_Input_Ctrl.bowtie.map > sampled_5fC_dIP-seq_Input_Ctrl.bowtie.map head -42637737 5fC_dIP-seq_Input_KD.bowtie.map > sampled_5fC_dIP-seq_Input_KD.bowtie.map head -42637737 5fC_dIP-seq_IP_KD.bowtie.map > sampled_5fC_dIP-seq_IP_KD.bowtie.map ** 5fC_dIP-seq_IP_Ctrl.bowtie.map had the fewest hits **
Call peaks with macs14
nohup macs14 -t 5fC_dIP-seq_IP_Ctrl.bowtie.map -c sampled_5fC_dIP-seq_Input_Ctrl.bowtie.map -f BOWTIE -g mm -n Control-Dipseq -w --call-subpeaks > macs14_ctrol_balanced & nohup macs14 -t sampled_5fC_dIP-seq_IP_KD.bowtie.map -c sampled_5fC_dIP-seq_Input_KD.bowtie.map -f BOWTIE -g mm -n KD-Dipseq -w --call-subpeaks > macs14_kd_balanced &
Peaks called
Positive Negative Result file: Control 354 24 File:Control Dipseq balanced peaks.xls Knockdown 301 21 File:KD Dipseq balanced peaks.xls
FDR < 5% BED file: Control File:Control Dipseq balanced peaks.BED.txt Total 107,302 bp Knockdown File:KD Dipseq balanced peaks.BED.txt Total 124,348 bp
Region_File | Control_Test_counts | Imputed_ave | Imputed_std | Empirical_p (greater than) | Empirical_p (less than) |
CpGIslands-regions.bed | 4991 | 453.284 | 533.1481 | 0 | 1 |
mm9_ucsc_genes.3utr.bed | 162 | 3183.955 | 2565.3198 | 0.935 | 0.065 |
mm9_ucsc_genes.5utr.bed | 5601 | 4194.234 | 2735.7066 | 0.236 | 0.764 |
mm9_ucsc_genes.exons.bed | 748 | 2570.128 | 1355.5966 | 0.954 | 0.0459999999999999 |
mm9_ucsc_genes.introns.bed | 45261 | 38383.062 | 7101.3337 | 0.161 | 0.839 |
mm9_ucsc_genes.tss.bed | 2002 | 949.66 | 869.7907 | 0.123 | 0.877 |
mouseESp300binding_LICR_regions.bed | 1187 | 1065.255 | 915.3372 | 0.387 | 0.613 |
REST-binding_sites.bed | 10 | 7.649 | 40.1048 | 0.085 | 0.915 |
Region_File | KD_Test_counts | Imputed_ave | Imputed_std | Empirical_p (greater than) | Empirical_p (less than) |
CpGIslands-regions.bed | 0 | 465.608 | 533.0195 | 1 | 0 |
mm9_ucsc_genes.3utr.bed | 2560 | 3584.865 | 2209.6552 | 0.64 | 0.36 |
mm9_ucsc_genes.5utr.bed | 1268 | 5157.965 | 2703.5931 | 0.952 | 0.0479999999999999 |
mm9_ucsc_genes.exons.bed | 192 | 3080.453 | 1407.9471 | 0.999 | 0.001 |
mm9_ucsc_genes.introns.bed | 39030 | 45091.197 | 6652.1637 | 0.811 | 0.189 |
mm9_ucsc_genes.tss.bed | 0 | 1169.726 | 903.6226 | 1 | 0 |
mouseESp300binding_LICR_regions.bed | 584 | 1142.344 | 904.8301 | 0.71 | 0.29 |
REST-binding_sites.bed | 11 | 10.553 | 40.473 | 0.116 | 0.884 |
Discussion[edit]
A) Genome browser snap shot of wig files around a called peak.
B) Pie charts showing base pairs genomic distribution of all significant peaks (FDR<5%) with annotated genomic regions.
C) Using wig files from input and IP, the mean occupancy of reads a 100bp sliding window around each of ~19K promoters were calculated. The mean values from input were subtracted from mean values of IP. Genes were separated into 5 groups based on the expression level of replicate #1 in the [GSE13805 | http://www.wip.ncbi.nlm.nih.gov/projects/geo/query/acc.cgi?acc=GSE13805] study on mouse ESCs.
- More peaks were found in knockdown than in wildtype, this supports that the knocked-down protein is responsible for converting 5fC into something else.
- In both control and knockdown, enrichment of peaks were in introns or other genomic regions and depletion of peaks (p<0.05) were in exons.
- In wildtype, occupancy of 5fC around promoters were not correlated with gene activity.
- In knockdown, occupancy of 5fC around promoters appears to be correlated with gene activity. Higher densities of 5fC were found towards the 5'UTR region of the gene than upstream of TSS.