#mix
#details @ http://nbviewer.ipython.org/github/sr320/ipython_nb/blob/master/fish546/BSseq_workflow.ipynb
!head /Volumes/web/cnidarian/YE_mix_22sm_methratio_out.txt

#control
#details @ http://nbviewer.ipython.org/github/sr320/ipython_nb/blob/master/fish546/BSseq_YE_control.ipynb
!head /Volumes/web/cnidarian/YE_control_22sm_methratio_out.txt

!grep "[A-Z][A-Z]CG[A-Z]" </Volumes/web/cnidarian/YE_mix_22sm_methratio_out.txt> /Volumes/web/cnidarian/YE_mix_22sm_methratio_CG.txt 

!head /Volumes/web/cnidarian/YE_mix_22sm_methratio_CG.txt

!grep "[A-Z][A-Z]CG[A-Z]" </Volumes/web/cnidarian/YE_control_22sm_methratio_out.txt> /Volumes/web/cnidarian/YE_control_22sm_methratio_CG.txt
!head /Volumes/web/cnidarian/YE_control_22sm_methratio_CG.txt

!awk '{print $1,$2,$2+1,$3,$8,($7/$8),(1-($7/$8))}' </Volumes/web/cnidarian/YE_mix_22sm_methratio_CG.txt> /Volumes/web/cnidarian/YE_mix_22sm_methkit_outCG.txt 

!head /Volumes/web/cnidarian/YE_mix_22sm_methkit_outCG.txt 

#convert space to tab
!tr ' ' "\t" </Volumes/web/cnidarian/YE_mix_22sm_methkit_outCG.txt> /Volumes/web/cnidarian/YE_mix_22sm_methkit.txt 

!head /Volumes/web/cnidarian/YE_mix_22sm_methkit.txt 

#processing control sample
!awk '{print $1,$2,$2+1,$3,$8,($7/$8),(1-($7/$8))}' </Volumes/web/cnidarian/YE_control_22sm_methratio_CG.txt> /Volumes/web/cnidarian/YE_control_22sm_methkit_outCG.txt
!tr ' ' "\t" </Volumes/web/cnidarian/YE_control_22sm_methkit_outCG.txt> /Volumes/web/cnidarian/YE_control_22sm_methkit.txt 
!head /Volumes/web/cnidarian/YE_control_22sm_methkit.txt 

#filter mk files for 3x coverage
!awk '{if ($5 >= 3) print $1,$2,$3,$4,$5,$6,$7}' </Volumes/web/cnidarian/YE_mix_22sm_methkit.txt> /Volumes/web/cnidarian/YE_mix_22sm_methkit3x.txt 
!head /Volumes/web/cnidarian/YE_mix_22sm_methkit3x.txt 

#filter mk files for 3x coverage
!awk '{if ($5 >= 3) print $1,$2,$3,$4,$5,$6,$7}' </Volumes/web/cnidarian/YE_control_22sm_methkit.txt> /Volumes/web/cnidarian/YE_control_22sm_methkit3x.txt 
!head /Volumes/web/cnidarian/YE_control_22sm_methkit3x.txt 

%pylab inline

%load_ext rpy2.ipython

%R library(methylKit)


%R library(data.table)


%R library(GenomicRanges)


!head /Volumes/web/cnidarian/YE_mix_22sm_methkit3x.txt 

%%R file.list <- list 
('/Volumes/web/cnidarian/YE_mix_22sm_methkit3x.txt',
 '/Volumes/web/cnidarian/YE_control_22sm_methkit3x.txt')

%R myobj=read(file.list,sample.id=list("mix","control"),assembly="v9",treatment=c(0,1))

%%R
meth<-unite(myobj)

%%R
meth<-unite(myobj)
head(meth)
nrow(meth)
getCorrelation(meth,plot=T)
hc<- clusterSamples(meth, dist="correlation", method="ward", plot=T)
PCA<-PCASamples(meth)