Dinh 2011/NOTES/2011-9-12

From ZhangLabWiki
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]

File:5fC Figure1.pdf

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.