Dinh/Dinh 2012/NOTES/2012-2-13: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Dinh
>Dinh
mNo edit summary
Line 19: Line 19:
  Then removed entries/segments which have over-lapping labels -> 3,187,848 entries left
  Then removed entries/segments which have over-lapping labels -> 3,187,848 entries left
  Then removed entries on chr*_random -> '''3,127,206 entries''' -> Non-overlapping repeat regions
  Then removed entries on chr*_random -> '''3,127,206 entries''' -> Non-overlapping repeat regions
* Sort BED and remove duplicates
for f in Indx*.BED
do
        sort -k1.4,1.4n -k2,2n -k3,3n -u -T tmp $f > rmdup.$f
done;


* Strand specific normalization with both IgG and Input controls. Normalization of read depth is taken care of by regression model.
* Strand specific normalization with both IgG and Input controls. Normalization of read depth is taken care of by regression model.
  for chr in chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chrX chrY chrM
  ftumch2 -i rmdup.Indx1.BED -c rmdup.Indx7.BED rmdup.Indx4.BED -o normalized.R.cac.BED -rl 350 -cs 100000 -ws 5 -iv 5
  do
  #ftumch2 -i rmdup.Indx2.BED -c rmdup.Indx7.BED rmdup.Indx4.BED -o normalized.R.fc.BED -rl 350 -cs 100000 -ws 5 -iv 5
        ./getCHR.pl $chr < Indx1.BED | sort -k2,3n > treat.BED
#ftumch2 -i rmdup.Indx3.BED -c rmdup.Indx7.BED rmdup.Indx4.BED -o normalized.R.hmc.BED -rl 350 -cs 100000 -ws 5 -iv 5
        ./getCHR.pl $chr < Indx7.BED | sort -k2,3n > input.BED
#ftumch2 -i rmdup.Indx5.BED -c rmdup.Indx7.BED rmdup.Indx6.BED -o normalized.R.mc.BED -rl 350 -cs 100000 -ws 5 -iv 5
        ./getCHR.pl $chr < Indx4.BED | sort -k2,3n > iggcl.BED
#ftumch2 -i rmdup.Indx8.BED -c rmdup.Indx14.BED rmdup.Indx11.BED -o normalized.S.cac.BED -rl 350 -cs 100000 -ws 5 -iv 5
        ftumch2 -i treat.BED -c input.BED iggcl.BED -o normalized.BED -rl 350 -cs 100000 -ws 5 -iv 5
#ftumch2 -i rmdup.Indx9.BED -c rmdup.Indx14.BED rmdup.Indx11.BED -o normalized.S.fc.BED -rl 350 -cs 100000 -ws 5 -iv 5
        cat normalized.BED >> normalized.R.cac.BED
#ftumch2 -i rmdup.Indx10.BED -c rmdup.Indx14.BED rmdup.Indx11.BED -o normalized.S.hmc.BED -rl 350 -cs 100000 -ws 5 -iv 5
  done >> log.Rcac
#ftumch2 -i rmdup.Indx12.BED -c rmdup.Indx14.BED rmdup.Indx13.BED -o normalized.S.mc.BED -rl 350 -cs 100000 -ws 5 -iv 5
  * repeat this for all the samples.
#ftumch2 -i rmdup.Indx15.BED -c rmdup.Indx19.BED rmdup.Indx18.BED -o normalized.R.dnmt3a.BED -rl 350 -cs 100000 -ws 5 -iv 5
#ftumch2 -i rmdup.Indx16.BED -c rmdup.Indx19.BED rmdup.Indx18.BED -o normalized.R.tdg.BED -rl 350 -cs 100000 -ws 5 -iv 5
#ftumch2 -i rmdup.Indx17.BED -c rmdup.Indx19.BED rmdup.Indx18.BED -o normalized.R.tet2.BED -rl 350 -cs 100000 -ws 5 -iv 5
#ftumch2 -i rmdup.Indx20.BED -c rmdup.Indx24.BED rmdup.Indx23.BED -o normalized.S.dnmt3a.BED -rl 350 -cs 100000 -ws 5 -iv 5
  #ftumch2 -i rmdup.Indx21.BED -c rmdup.Indx24.BED rmdup.Indx23.BED -o normalized.S.tdg.BED -rl 350 -cs 100000 -ws 5 -iv 5
#ftumch2 -i rmdup.Indx22.BED -c rmdup.Indx24.BED rmdup.Indx23.BED -o normalized.S.tet2.BED -rl 350 -cs 100000 -ws 5 -iv 5
 


==BEDtools==
==BEDtools==
* Used bedTools to convert all my BED (mapped reads) to bedGraph format (this is for validation on raw data)
* Used bedTools to convert all my BED (mapped reads) to bedGraph format (this is for validation on raw data)
  for f in sorted*BED
  for f in rmdup*BED
  do
  do
         genomeCoverageBed -bg -i $f -g mm9.fai > $f.bdg
         genomeCoverageBed -bg -i $f -g mm9.fai > $f.bdg

Revision as of 07:56, 15 February 2012

24 ChIP-seq - Stimulated and Resting Neurons

  • Analysis at repeats.
Downloaded RepeatMasker track fro mm9. Then filtered and get Non-overlapping windows:
 2072  awk '{if($12 ~ /DNA/) print $0}' RepeatMasker_mm9.txt | ./getNonOverlapping.pl > DNA_repeats.nonoverlapping &
2073  awk '{if($12 ~ /LINE/) print $0}' RepeatMasker_mm9.txt | ./getNonOverlapping.pl > LINE_repeats.nonoverlapping &
2074  awk '{if($12 ~ /SINE/) print $0}' RepeatMasker_mm9.txt | ./getNonOverlapping.pl > SINE_repeats.nonoverlapping &
2075  awk '{if($12 ~ /LTR/) print $0}' RepeatMasker_mm9.txt | ./getNonOverlapping.pl > LTR_repeats.nonoverlapping &
2076  awk '{if($12 ~ /Satellite/) print $0}' RepeatMasker_mm9.txt | ./getNonOverlapping.pl > Satellite_repeats.nonoverlapping &
File:GetNonOverlapping.txt
  • Note that some repeat segments overlaps between these difference repeat class.
cat *_repeats.nonoverlapping | sort -k1,1 -k2,3n > sorted.Repeats (3,188,316 entries)
./removeOverlapping < sorted.Repeats > tmp
mv tmp sorted.Repeats (3,187,848 entries)
File:RemoveOverlapping.txt
  • Remove chr*_random entries
awk '{if($1 !~ /_random/) print $0}' sorted.Repeats > sorted.Repeats_norandom
  • Final repeat regions file: 3,127,206 entries
Started with 3,188,316 non-overlapping entries for DNA, LINE, SINE, LTR, and Satellite repeats
Then removed entries/segments which have over-lapping labels -> 3,187,848 entries left
Then removed entries on chr*_random -> 3,127,206 entries -> Non-overlapping repeat regions
  • Sort BED and remove duplicates
for f in Indx*.BED 
do
       sort -k1.4,1.4n -k2,2n -k3,3n -u -T tmp $f > rmdup.$f
done;
  • Strand specific normalization with both IgG and Input controls. Normalization of read depth is taken care of by regression model.
ftumch2 -i rmdup.Indx1.BED -c rmdup.Indx7.BED rmdup.Indx4.BED -o normalized.R.cac.BED -rl 350 -cs 100000 -ws 5 -iv 5
#ftumch2 -i rmdup.Indx2.BED -c rmdup.Indx7.BED rmdup.Indx4.BED -o normalized.R.fc.BED -rl 350 -cs 100000 -ws 5 -iv 5
#ftumch2 -i rmdup.Indx3.BED -c rmdup.Indx7.BED rmdup.Indx4.BED -o normalized.R.hmc.BED -rl 350 -cs 100000 -ws 5 -iv 5
#ftumch2 -i rmdup.Indx5.BED -c rmdup.Indx7.BED rmdup.Indx6.BED -o normalized.R.mc.BED -rl 350 -cs 100000 -ws 5 -iv 5
#ftumch2 -i rmdup.Indx8.BED -c rmdup.Indx14.BED rmdup.Indx11.BED -o normalized.S.cac.BED -rl 350 -cs 100000 -ws 5 -iv 5
#ftumch2 -i rmdup.Indx9.BED -c rmdup.Indx14.BED rmdup.Indx11.BED -o normalized.S.fc.BED -rl 350 -cs 100000 -ws 5 -iv 5
#ftumch2 -i rmdup.Indx10.BED -c rmdup.Indx14.BED rmdup.Indx11.BED -o normalized.S.hmc.BED -rl 350 -cs 100000 -ws 5 -iv 5
#ftumch2 -i rmdup.Indx12.BED -c rmdup.Indx14.BED rmdup.Indx13.BED -o normalized.S.mc.BED -rl 350 -cs 100000 -ws 5 -iv 5
#ftumch2 -i rmdup.Indx15.BED -c rmdup.Indx19.BED rmdup.Indx18.BED -o normalized.R.dnmt3a.BED -rl 350 -cs 100000 -ws 5 -iv 5
#ftumch2 -i rmdup.Indx16.BED -c rmdup.Indx19.BED rmdup.Indx18.BED -o normalized.R.tdg.BED -rl 350 -cs 100000 -ws 5 -iv 5
#ftumch2 -i rmdup.Indx17.BED -c rmdup.Indx19.BED rmdup.Indx18.BED -o normalized.R.tet2.BED -rl 350 -cs 100000 -ws 5 -iv 5
#ftumch2 -i rmdup.Indx20.BED -c rmdup.Indx24.BED rmdup.Indx23.BED -o normalized.S.dnmt3a.BED -rl 350 -cs 100000 -ws 5 -iv 5
#ftumch2 -i rmdup.Indx21.BED -c rmdup.Indx24.BED rmdup.Indx23.BED -o normalized.S.tdg.BED -rl 350 -cs 100000 -ws 5 -iv 5
#ftumch2 -i rmdup.Indx22.BED -c rmdup.Indx24.BED rmdup.Indx23.BED -o normalized.S.tet2.BED -rl 350 -cs 100000 -ws 5 -iv 5


BEDtools

  • Used bedTools to convert all my BED (mapped reads) to bedGraph format (this is for validation on raw data)
for f in rmdup*BED
do
       genomeCoverageBed -bg -i $f -g mm9.fai > $f.bdg
done;

SISSRs

  • Peak calling using SISSRs (a simple perl script)..FDR calculated from a Poisson distribution.
for f in normalized*.BED
do
       sissrs.pl -i $f -s 2654911517 -m 0.7 -F 350 -u -D 0.001 -o $f.peaks
done;
-i (input BED file)
-s (genome size in bp)
-m (fraction of genome covered, estimated 0.7)
-F (average fragment length, inferred from sampe.log (BWA) the median is 300, and average could be 350.
-u (report peaks which is not supported by at least two reads in both directions as well)
-D (FDR, 0.001)
-o (output peak file)
  • Peaks for properly paired reads:
 40926 normalized.R.cac.BED.peaks.bed
   819 normalized.R.dnmt3a.BED.peaks.bed
 12266 normalized.R.fc.BED.peaks.bed
 14761 normalized.R.hmc.BED.peaks.bed
 23445 normalized.R.mc.BED.peaks.bed
  1299 normalized.R.tdg.BED.peaks.bed
  2189 normalized.R.tet2.BED.peaks.bed
 61048 normalized.S.cac.BED.peaks.bed
  1431 normalized.S.dnmt3a.BED.peaks.bed
  8990 normalized.S.fc.BED.peaks.bed
 51264 normalized.S.hmc.BED.peaks.bed
 43123 normalized.S.mc.BED.peaks.bed
  1221 normalized.S.tdg.BED.peaks.bed
  1096 normalized.S.tet2.BED.peaks.bed
  • Peaks for improperly paired reads
 15174 normalized.R.cac.BED.peaks.bed
    65 normalized.R.dnmt3a.BED.peaks.bed
  3599 normalized.R.fc.BED.peaks.bed
  1526 normalized.R.hmc.BED.peaks.bed
  1827 normalized.R.mc.BED.peaks.bed
    62 normalized.R.tdg.BED.peaks.bed
    85 normalized.R.tet2.BED.peaks.bed
 15010 normalized.S.cac.BED.peaks.bed
    72 normalized.S.dnmt3a.BED.peaks.bed
  1438 normalized.S.fc.BED.peaks.bed
   981 normalized.S.hmc.BED.peaks.bed
  1091 normalized.S.mc.BED.peaks.bed
    71 normalized.S.tdg.BED.peaks.bed
    69 normalized.S.tet2.BED.peaks.bed

FindPeak 4.0

  • Peak calling using FindPeak 4.0. The reason why MACS1.4 wouldn't work because MACS14 "shifts" read toward the center when reads are pileup on the forward strand and reverse strand. This shifting model wouldn't work for finding enrichment of certain regions in the genome (unlike motif finding, the size of the region is unknown and may be very large).
for f in normalized.*BED
do
       awk '{print $1"\t"$2"\t"$3"\t"$4"\t0\t"$6}' $f > tmp
       java -Xmx2G -jar /usr/local/fp4/FindPeaks.jar -input tmp -aligner bed -output . -wig_step_size 1 - dist_type 0 350 -subpeaks 0.5 -landerwaterman 0.001 -name $f -one_per
done;
  • Number of peaks called for each dataset:
  41109 normalized.R.cac.BED_fixed_350_subpeaks.peaks
    820 normalized.R.dnmt3a.BED_fixed_350_subpeaks.peaks
  12382 normalized.R.fc.BED_fixed_350_subpeaks.peaks
  14825 normalized.R.hmc.BED_fixed_350_subpeaks.peaks
  28826 normalized.R.mc.BED_fixed_350_subpeaks.peaks
   1301 normalized.R.tdg.BED_fixed_350_subpeaks.peaks
   2177 normalized.R.tet2.BED_fixed_350_subpeaks.peaks
  61289 normalized.S.cac.BED_fixed_350_subpeaks.peaks
   1431 normalized.S.dnmt3a.BED_fixed_350_subpeaks.peaks
   9107 normalized.S.fc.BED_fixed_350_subpeaks.peaks
  51385 normalized.S.hmc.BED_fixed_350_subpeaks.peaks
  52640 normalized.S.mc.BED_fixed_350_subpeaks.peaks
   1222 normalized.S.tdg.BED_fixed_350_subpeaks.peaks
   1090 normalized.S.tet2.BED_fixed_350_subpeaks.peaks
  • Things to note: For all properly paired reads there are fewer peaks for ChIP data than for dIP data. This is expected since the ChIP are peaks for DNA binding proteins (they are specific). Whereas for the dIP data, it is possible to have more signals since there are many location for cytosines to be oxydized.
  • I also used the same normalization and peak calling method (except with dist_type 0 110) for the improperly paired reads.
 16759 normalized.R.cac.BED_fixed_110_subpeaks.peaks
    95 normalized.R.dnmt3a.BED_fixed_110_subpeaks.peaks
  4250 normalized.R.fc.BED_fixed_110_subpeaks.peaks
  1720 normalized.R.hmc.BED_fixed_110_subpeaks.peaks
  2047 normalized.R.mc.BED_fixed_110_subpeaks.peaks
    94 normalized.R.tdg.BED_fixed_110_subpeaks.peaks
   130 normalized.R.tet2.BED_fixed_110_subpeaks.peaks
 16148 normalized.S.cac.BED_fixed_110_subpeaks.peaks
   112 normalized.S.dnmt3a.BED_fixed_110_subpeaks.peaks
  1688 normalized.S.fc.BED_fixed_110_subpeaks.peaks
  1097 normalized.S.hmc.BED_fixed_110_subpeaks.peaks
  1212 normalized.S.mc.BED_fixed_110_subpeaks.peaks
   109 normalized.S.tdg.BED_fixed_110_subpeaks.peaks
   104 normalized.S.tet2.BED_fixed_110_subpeaks.peaks
  • Things to note: much fewer peaks found for ChIP data, which is expected, since DNMT3a, TDG, and Tet2 are not known to bind in repeats. Up to 10x greater signals are observed for caC marks in both resting and stimulated neurons.