MHL vs Methylation Freqency from PileOMeth

From ZhangLabWiki
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()