/home/shg047/oasis/monod/bin/LDR2extent.pl

From ZhangLabWiki
Jump to navigation Jump to search
#!/usr/bin/perl -w
# Transfer Bam to Fastq with samtools command
# Run the script to the Bam directory
# Contact: Shicheng Guo
# Version 1.3
# Update: 2016-02-19
use strict;
use Cwd;
chdir getcwd;
# my $region=shift @ARGV;
& printUsage if scalar(@ARGV)<2;
my $region=shift @ARGV; # chr5:112,039,140-112,047,142
my $r2File = shift @ARGV; # R2 Matrix = "APC.chr5.rsq"
my $cpg_position_file="/home/kunzhang/HsGenome/hg19/HsGenome19.CpG.positions.txt"; # genome-miner
my @region=locationSplit($region);
# my $cpg_position_file="/home/shg047/oasis/db/hg19/meth/bismark/HsGenome19.CpG.positions.txt"; # TSCC
open F,$cpg_position_file;
chomp(my @loci=<F>);
close F;
my $Loci=\@loci;
my @lociByRegion=lociPickUpByRegion($region,$Loci);
my %R2;
foreach my $loc1(@lociByRegion){
foreach my $loc2(@lociByRegion){
$R2{$loc1}{$loc2}="NA";
}
}
open F,$r2File;
my @R2loci;
while(<F>){
chomp;
if (/^\s/){
@R2loci=split /\t/;
next;
}
my @line=split /\s+/;
foreach my $tmp(1..$#line){
$R2{$line[0]}{$R2loci[$tmp]}=$line[$tmp];
}
}
foreach my $loc1(@lociByRegion){
my @tmp;
foreach my $loc2(@lociByRegion){
push(@tmp,$R2{$loc1}{$loc2});
}
my $print=join("\t",@tmp);
print "$print\n";
}
sub locationSplit{
my $location=shift;
$location=~s/,//g;
my ($chr, $start, $end) = split /[:-]/, $location;
return($chr,$start,$end);
}
sub lociPickUpByRegion{
my $region=shift;
my $Loci=shift;
my @loci=@$Loci;
my @lociPickup;
my($chr1,$start1,$end1)=locationSplit($region);
foreach my $loci(@loci){
my ($chr2,$start2)=split/\t/,$loci;
if ($chr2 eq $chr1 and $start2>=$start1 and $start2<= $end1){
push(@lociPickup,$start2);
}
} 
# @lociPickup=sort(@lociPickup);
return(@lociPickup)
}
sub printUsage{
print " Usage:\n";
print " perl $0 <chr:start-end> <R-square matrix file>\n";
print " For example:\n perl $0 chr10:10000873-10001472 APC.chr5.rsq\n\n";
print "-----------------------------------------------------------------------------\n";
print " APC.chr5.rsq format:\n";
print " 10000873 10001472 10001596\n";
print "10000873 0.3 0.8 0.1 \n";
print "10001472 0.5 0.4 0.9 \n";
print "10001596 0.2 0.1 0.8 \n";
print "-----------------------------------------------------------------------------\n";
exit 0;
}