~/oasis/monod/bin/hapinfo2R2ByBed.pl

From ZhangLabWiki
Jump to navigation Jump to search
#!/usr/bin/perl -w
# this script reads in a haploInfo file and output the rsq value for each pair of cg sites
# Transfer Bam to Fastq with samtools command
# Run the script to the Bam directory
# Contact: Shicheng Guo
# Version 1.3
# Update: 2016-02-25
use strict;
use strict;
use warnings;
use Sort::Array qw/Sort_Table/;
my %cgTable;
my %RsqTable;
my %RsqPvalTable;
my %posList;
my %methylTable;
my $min_coverage = 2;
my $min_rsq = 0.1;
my $pvalue_cutoff = 0.1;
my $distance_threshold = 200;
my %rcTable;
$rcTable{'A'}='T';
$rcTable{'T'}='A';
$rcTable{'G'}='C';
$rcTable{'C'}='G';
$rcTable{'N'}='N';
$rcTable{'D'}='D';
printUsage() if(!$ARGV[0]);
my ($target_chr, $target_start, $target_end) = split /[:-]/, $ARGV[1];
sub main{
while(my $line = <STDIN>){
next if $line=~/^\s+$/;
chomp($line);
my @tmp = split /\s+/, $line;
next if scalar(@tmp)<4;
my ($chr, $start, $end) = split /[:-]/, $tmp[0];
next if($chr ne $target_chr || $start > $target_end || $end < $target_start);
my ($hapString, $hapCount) = ($tmp[1], $tmp[2]);
my @cgPos = split ",", $tmp[3];
next if(length($tmp[1]) ne scalar(@cgPos));
my %indices;
for(my $i = 0; $i < scalar(@cgPos); $i++){
my $pos = $cgPos[$i];
my $call = substr($hapString,$i,1);
if($indices{$pos}){
# pos has already been assigned?
die("Getting ambiguous call for position $pos for:\n$line\n");
delete $indices{$pos};
}else{
next if($call !~ /[CT]/);
$indices{$pos} = $call;
}
}
my @ordered_pos = sort {$a<=>$b} keys %indices;
next if(scalar(@ordered_pos) < 2);
#print join(",", @ordered_pos), "\n";
for(my $i = 0; $i < scalar(@ordered_pos); $i++){
for(my $j = $i+1; $j < scalar(@ordered_pos); $j++){
my $value_key = $chr.":".$ordered_pos[$i]."-".$ordered_pos[$j];
my $value_call = $indices{$ordered_pos[$i]}.$indices{$ordered_pos[$j]};
#$value_call =~ tr/T/U/;
#$value_call =~ tr/C/M/;
my ($mm, $mu, $um, $uu) = (0,0,0,0);
my ($rsq, $pval, $counts, $totalHaps) = ("NA", 0, "NA", 0);
if($RsqTable{$value_key}){
($rsq, $pval, $counts,$totalHaps) = split "\t", $RsqTable{$value_key};
($mm, $mu, $um, $uu) = split ":", $counts;
}
$mm+=$hapCount if($value_call eq "CC");
$mu+=$hapCount if($value_call eq "CT");
$um+=$hapCount if($value_call eq "TC");
$uu+=$hapCount if($value_call eq "TT");
$totalHaps+=$hapCount;
$counts = "$mm:$mu:$um:$uu";
if($totalHaps >= $min_coverage){
# uniform site is where one cpg have maf of 0
my $mafA = ($um + $uu) < ($mu + $mm) ? $um + $uu : $mu + $mm;
my $mafB = ($mu + $uu) < ($um + $mm) ? $mu + $uu : $um + $mm;
if($mafA ne 0 and $mafB ne 0){
($rsq, $pval) = updateRsq($counts);
$RsqTable{$value_key} = "$rsq\t$pval\t$counts\t$totalHaps";
}else{
$rsq = 1 if($mafA + $mafB == 0);
$rsq = ($mu+$um)/$mm if($uu == 0 && $mafA + $mafB > 0);
$rsq = ($mu+$um)/$uu if($mm == 0 && $mafA + $mafB > 0);
if($rsq ne "NA"){
$rsq = 1/$rsq if($rsq > 1);
$rsq = $rsq > (1-$rsq) ? $rsq : 1 - $rsq;
my $out_rsq = sprintf("%4.3f", $rsq);
$RsqTable{$value_key} = "$out_rsq\tMAF_0\t$counts\t$totalHaps";
}else{
$RsqTable{$value_key} = "NA\tMAF_0\t$counts\t$totalHaps";
}
}
}else{
$RsqTable{$value_key} = "$rsq\t$pval\t$counts\t$totalHaps";
}
}
}
undef(@ordered_pos);
undef(%indices);
}
outputRsqTable();
}
sub readCGTable{
open(IN, "$ARGV[1]") || die("Error opening $ARGV[1]\n");
while(my $line = <IN>){
chomp($line);
my @tmp = split /\t/, $line;
$tmp[0] =~ s/:W//;
$cgTable{$tmp[0].":".$tmp[1]} = 1;
}
close(IN);
}
sub reportBlocks{
my $filename = shift;
my $h_blocks = shift;
my %blocktable = %{$h_blocks};
open(OUTFILE, ">$filename") || die("Error in opening file $filename\n");
print OUTFILE "block id\tstart locus\tend locus\n";
foreach my $id (sort keys(%blocktable)) {
next if($blocktable{$id}->{'end'} - $blocktable{$id}->{'start'} < 100);
print OUTFILE $id, "\t", $blocktable{$id}->{'start'}, "\t", $blocktable{$id}->{'end'}, "\n";
}
close(OUTFILE);
}
sub outputRsqTable{
# open(RSQ_OUT, ">$ARGV[0]") || die("error writing to file $ARGV[0]\n");
my %table;
#print RSQ_OUT "index\tRsq\tpval\tmm:mu:um:uu\ttotalHaps\n";
foreach my $pairs(keys %RsqTable){
my ($chr, $start, $end) = split /[:-]/, $pairs;
my ($rsq, $pval, $counts, $totalHaps) = split "\t", $RsqTable{$pairs};
next if($totalHaps < $min_coverage);
$posList{$chr}->{$start} = 1;
$posList{$chr}->{$end} = 1;
#print RSQ_OUT $pairs, "\t", $RsqTable{$pairs},"\n";
$table{$start}->{$end} = $rsq;
$table{$end}->{$start} = $rsq;
}
foreach my $chr (keys %posList){
open(RSQ_OUT, ">$ARGV[0].$chr.rsq") || die("error writing to file $ARGV[0].chr.rsq\n");
my @sorted_list = sort {$a<=>$b} keys %{$posList{$chr}};
print RSQ_OUT "\t", join("\t", @sorted_list), "\n";
for(my $i = 0; $i < scalar(@sorted_list); $i++){
print RSQ_OUT $sorted_list[$i];
for(my $j = 0; $j < scalar(@sorted_list); $j++){
if(!$table{$sorted_list[$i]}->{$sorted_list[$j]}){
print RSQ_OUT "\tNA";
next;
}
print RSQ_OUT "\t", $table{$sorted_list[$i]}->{$sorted_list[$j]};
}
print RSQ_OUT "\n";
}
close(RSQ_OUT);
}
}
sub updateRsq{
my $counts = shift;
my ($mm, $mu, $um, $uu) = split ":", $counts;
my $total = $mm+$mu+$um+$uu;
#calculate p1 = prob(1M), p2=prob(1U), q1=prob(2M), q2=prob(2U)
#calculate D = freq(1M/2M) - prob(1M)*prob(2M)
#calculate r^2 = D^2/(p1*p2*q1*q2)
my ($p1, $q1, $p2, $q2) = ($mm+$mu, $mm+$um, $um+$uu, $mu+$uu);
$p1 /=$total;
$q1 /=$total;
$p2 /=$total;
$q2 /=$total;
$mm /=$total;
$mu /=$total;
$um /=$total;
$uu /=$total;
my $D = $mm*$uu - $um*$mu;
my $abs_Dprime;
my $Dmax = $D > 0 ? ($p1*$q1 < $p2*$q2 ? $p1*$q1:$p2*$q2) : ($p1*$q2 < $p2*$q1 ? $p1*$q2:$p2*$q1);
my $rsq = 0;
if($D == 0){
$abs_Dprime = ($Dmax == 0.0) ? 1.0 : 0.0;
$rsq = ($p1*$p2*$q1*$q2) == 0 ? 1.0 : 0.0;
#$abs_d = ($q1*$q2 == 0) ? 1.0 : 0.0;
#$abs_Q = ($mm*$uu + $um*$mu == 0.0) ? 1.0 : 0.0;
}else{
$abs_Dprime = abs($D/$Dmax);
$rsq = ($D*$D)/($p1*$p2*$q1*$q2);
#$abs_d = abs($D/($q1*$q2));
#$abs_Q = abs($D/($mm*$uu + $um*$mu));
}
my $out_rsq = sprintf("%4.3f", $rsq);
my $out_pval;
my ($mm_s, $mu_s, $um_s, $uu_s) = (int($min_coverage*$mm), int($min_coverage*$mu), int($min_coverage*$um), int($min_coverage*$uu));
my $counts_s = "$mm_s:$mu_s:$um_s:$uu_s";
if(defined($RsqPvalTable{$counts_s})){
return $out_rsq, $RsqPvalTable{$counts_s};
}
# begins permutation:
my $num_perm = 500;
$p1 = sprintf("%.3f", $p1);
$q1 = sprintf("%.3f", $q1);
my $ptotal = $min_coverage;
my $greater_eq_than = 0;
for(my $i = 0; $i < $num_perm; $i++){
my ($r_mm, $r_mu, $r_um, $r_uu) = (0,0,0,0);
for(my $j = 0; $j < $ptotal; $j++){
my $a = int(rand(1001))/1000;
my $b = int(rand(1001))/1000;
$r_mm++ if($a <= $p1 and $b <= $q1);
$r_mu++ if($a <= $p1 and $b > $q1);
$r_um++ if($a > $p1 and $b <= $q1);
$r_uu++ if($a > $p1 and $b > $q1);
}
$r_mm /= $ptotal;
$r_mu /= $ptotal;
$r_um /= $ptotal;
$r_uu /= $ptotal;
my $r_D = $r_mm*$r_uu - $r_um*$r_mu;
my $r_rsq = ($r_D*$r_D)/($p1*$p2*$q1*$q2);
$greater_eq_than++ if($r_rsq >= $rsq);
}
$out_pval = sprintf("%4.3f", $greater_eq_than/$num_perm);
$RsqPvalTable{$counts_s} = $out_pval;
return $out_rsq, $out_pval;
}
sub lookupLD{
my $chr = shift;
my $pos1 = shift; # first cpg position
my $pos2 = shift; # last cpg position
my $rsq;
if(!$RsqTable{$chr.":".$pos1."-".$pos2}){
$rsq = "NA"; # no evidence for LD so split block
}else{
my ($rsq_out, $pval, $allele_count, $totalHaps) = split "\t", $RsqTable{$chr.":".$pos1."-".$pos2};
$rsq = $rsq_out;
$rsq = 0 if($pval ne "MAF_0" and $pval > $pvalue_cutoff);
}
return $rsq;
}
sub getBlocksGreedy{
my $sorted_list = shift;
my $chr = shift;
my $start = shift;
my $end = shift;
my $threshold = shift;
my $i = $start;
my @blocks;
my %blockSetTable;
print "Total nummber of CpGs is $end for $chr\n";
while($i <= $end){
# continue adding CpGs to set if there is another CpG within distance threshold.
while(($i < $end) and @{$sorted_list}[$i+1] - @{$sorted_list}[$i] < $distance_threshold){
$i++;
}
print $i, "\n";
push(@blocks, getAllBlockInRegion($sorted_list, $chr, $start, $i, $threshold));
$i++;
$start = $i;
}
my @sortedBlocks = Sort_Table(
cols => '2',
field => '1',
sorting => 'ascending',
structure => 'csv',
separator => '\:',
data => \@blocks,
);
for(my $i=0; $i<scalar(@blocks); $i++){
my $id = sprintf("B%05d", $i+1);
my @words = split(/:/, $sortedBlocks[$i]);
$blockSetTable{$id}->{'start'} = $words[0];
$blockSetTable{$id}->{'end'} = $words[1];
}
return %blockSetTable;
}
sub getAllBlockInRegion{
my $sorted_list = shift;
my $chr = shift;
my $start = shift;
my $end = shift;
my $threshold = shift;
my @blocks;
if($start == $end){
push(@blocks, @{$sorted_list}[$start] . ":" . @{$sorted_list}[$end]);
return @blocks;
}
my ($block_start, $block_end) = findMaxBlockInRegion($sorted_list, $chr, $start, $end, $threshold);
push(@blocks, @{$sorted_list}[$block_start] . ":" . @{$sorted_list}[$block_end]) if($block_end - $block_start > 2);
# if the largest block is not all inclusive:
if($block_start > $start){
my @sub_blocks = getAllBlockInRegion($sorted_list, $chr, $start, $block_start-1, $threshold);
push(@blocks, @sub_blocks);
}
if($block_end < $end){
my @sub_blocks = getAllBlockInRegion($sorted_list, $chr, $block_end+1, $end, $threshold);
push(@blocks, @sub_blocks);
}
return @blocks;
}
sub findMaxBlockInRegion(){
my $sorted_list = shift;
my $chr = shift;
my $start = shift;
my $end = shift;
my $threshold = shift;
my $max_block_start=$start;
my $max_block_end=$start;
return ($start, $end) if($start == $end);
for(my $size = $end-$start+1; $size >1; $size--){ # try to get the largest size first
my $good_block=0;
for(my $i= $start; $i<= $end-$size+1; $i++){ # allow every pos to be starting point
$good_block = 1;
for(my $j= $i; $j<$i+$size; $j++){ # every pos to every other pos within block must be in LD!
for(my $k = $j+1; $k<$i+$size; $k++){
if(lookupLD($chr, @{$sorted_list}[$j], @{$sorted_list}[$k]) ne "NA") {
$good_block = 0 if(lookupLD($chr, @{$sorted_list}[$j], @{$sorted_list}[$k]) < $threshold);
last;
}
}
}
if($good_block){
$max_block_start = $i;
$max_block_end = $i + $size -1;
last;
}
}
last if($good_block);
}
return ($max_block_start, $max_block_end);
}
sub printUsage{
print " Usage:\n";
print " perl ~/monod/bin/$0 <rlt_prefix> <chr:start-end> < haploinfo.input.txt\n";
print " For example:\n perl $0 rlt chr10:10000873-10001472 < haploinfo.input.txt\n\n";
print "-----------------------------------------------------------------------------\n";
print " haploinfo.input.txt format:\n";
print " chr10:10000873-10001472 CCC 1 10001056,10001082,10001168\n\n";
print "-----------------------------------------------------------------------------\n";
exit 0;
}
main();