MHL vs Methylation Freqency from PileOMeth
Jump to navigation
Jump to search
# perl #!/usr/bin/perl use Cwd; use strict; my $file="/home/shg047/oasis/DennisLo2015/mhl/Dennis.bed"; open F,$file; while(<F>){ chomp; system("PileOMeth extract -p 5 -q 10 --minDepth 1 -r $_ /home/shg047/oasis/db/hg19/hg19.fa BMT1.read1_val_1.fq.gz_bismark_bt2_pe.sort.bam -o Tmp"); my $amf=qx/grep -v track Tmp_CpG.bedGraph | awk '{cov+=\$5+\$6;nC+=\$5}END{print nC\/cov}'/; print "$_\t$amf"; }
# R library(RColorBrewer) x<-read.table("/home/shg047/oasis/DennisLo2015/sortbam/BMT1.read1_val_1.fq.gz_bismark_bt2_pe.sort.PileOMeth.MF.txt",row.names=1) y<-read.table("/home/shg047/oasis/DennisLo2015/mhl/dennis.hapinfo2mhl.march28.txt",head=T,row.names=1) z<-read.table("/home/shg047/oasis/DennisLo2015/mhl/dennis.hapinfo2mf.march28.txt",head=T,row.names=1) x1<-x[1:nrow(x),1] y1<-y[1:nrow(x),1] z1<-z[1:nrow(x),1] smoothScatter(x1, y1, nrpoints=0.05*length(x1), colramp=colorRampPalette(rev(brewer.pal(11, "RdYlBu"))), pch=19, cex=.3, col = "green1",xlab="PileOMeth-MF",ylab="MHL") dev.off() smoothScatter(x1, z1, nrpoints=0.05*length(x1), colramp=colorRampPalette(rev(brewer.pal(11, "RdYlBu"))), pch=19, cex=.3, col = "green1",xlab="PileOMeth-MF",ylab="Hapinfo-MF") dev.off()