Elementolab/RRBSseeqer Tutorial
From Icbwiki
Elementolab/ Elementolab/RRBSseeqer/
Got to RRBS directory
cd RRBS
Contents |
Detect differentially methylated CGs
./RRBSseeqer_CG -rrbs1 sample.calls.gz -rrbs2 sample_control.calls.gz | sort_column.pl 1 | sort_column_alpnum.pl 0 > cgs.txt
Options:
-rrbs1 FILE (condition 1) -rrbs2 FILE (condition 2) -test [fisher|chi] (Statistical test to be performed) -minfold FLOAT (Def: 1.25) -fdr FLOAT (Def: 0.2)
Detect differentially methylated regions (DMR)
perl RRBSidentifyUpDownDMR.pl --metfile=cgs.txt --dmax=250 --minmetdx=0.1 --minnumcg=5 -outfile=cgs_DMR.txt
Input must come from RRBSeeqer_CG. Parameters indicate maximum distance between diff methylated CGs (--dmax) and minimum number of diff meth CGs for a region to be reported
To create a bed file with all DMRs
perl DMR2bed.pl cgs_DMR.txt > cgs_DMR.bed
To create a bed file that can be uploaded to the Genome Browser, with two different tracks for hyper and hypo DMRs
perl DMR2bed2.pl cgs_DMR.txt
Annotate CGs
With CpG islands
./RRBSannotate -peakfile1 cgs.txt -peakfile2 DATA/hg19/CpGislands_annotation -showpeakdesc 1 -hasid1 0 -outfile cgs_CGI.txt
With RRBSannotate, you can just add annotation columns. Here's how to add shore annotation
./RRBSannotate -peakfile1 cgs_CGI.txt -peakfile2 DATA/hg19/CpGislands_annotation -shores 1 -ext2 1000 -hasid2 0 -showpeakdesc 1 -hasid1 0 -outfile cgs_CGI_shores.txt
With genes
./RRBSannotateGenes -rrbsfile cgs.txt -annotation DATA/hg19/refSeq -chrdata DATA/hg19.chrdata
(you need to download gene annotation files first (e.g., hg18.tar.gz, hg19.tar.gz) from --> http://physiology.med.cornell.edu/faculty/elemento/lab/CS_files/)
Draw methylation across chromosomes
perl DrawRRBSmap2.pl --chrdata=DATA/hg19.chrdata --rrbsdata=cgs.txt
Make BigWig DIFFERENTIAL methylation track
Let's assume ben1_pca1.txt contains a differential methylation analysis generated by RRBSseeqer_CG above. To generate a BigWig track, make sure wigToBigWig is installed, and type something like this
RRBSmakeMethylDiffTrack -methfile ben1_pca1.txt -desc ben1_pca1 -type diff -bigwig 1 -chrdata hg19.chrdata -chromsizes hg19.chroms.sizes -s 25
This command will create a .bw track (binary format), with window size 50bp (2x25) that you can upload to your web server, then make visible to the UCSC browser using:
track type=bigWig name="ben1_pca1" description="ben1_pca1" bigDataUrl=http://physiology.med.cornell.edu/faculty/elemento/lab/files/ben1_pca1.wig.bw \ visibility="full" maxHeightPixels="64:64:11" smoothingWindow="10" viewLimits="0:100" autoScale="off"
A sometimes faster way to visualize a BigWig track is to load it into IGV http://www.broadinstitute.org/igv/
Methylation levels for annotated elements
getMethylationLevelsForAnnotatedElements.pl --callfile=CpG_context_Sample_ln_cap_r.fastq.filtered.noadapt.gz_bismark.txt.calls_NM --headerhasid=0 \ --promlist=refSeq.new.u2000_d2000.PROMOTERS --minnumcg=2 --col=5
