Perl ~/bin/hapinfo2mhl.pl: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Shicheng
No edit summary
>Shicheng
No edit summary
 
Line 4: Line 4:
  # Contact: Shicheng Guo
  # Contact: Shicheng Guo
  # Version 1.3
  # Version 1.3
  # Update: 2016-03-16
  # Update: 2016-03-31
  use strict;
  use strict;
  use Cwd;
  use Cwd;
  my $usage = <<USAGE;
  my $usage = <<USAGE;
  perl haplo2hml.pl Directory_of_haploinfo
  perl $0 Hapinfo_File_List.txt > MHL.output.txt
  USAGE
  USAGE
  die $usage if @ARGV <1;
  die $usage if @ARGV <1;
Line 14: Line 14:
  my %probe_HMH_samples;
  my %probe_HMH_samples;
  my %hap_count_matrix;
  my %hap_count_matrix;
  chdir shift @ARGV;
  my $hapinfList=shift @ARGV;
  my @hapInfo_files= glob("*.hapInfo.txt");
  open FF,$hapinfList;
chomp(my @hapInfo_files=<FF>);
close FF;
  my @sample_list;
  my @sample_list;
  foreach my $hapInfo_file(@hapInfo_files){
  foreach my $hapInfo_file(@hapInfo_files){
  my $sample_name = $hapInfo_file;
my @line=split /\//,$hapInfo_file;
  my $sample_name = $line[$#line];
  $sample_name =~ s/.hapInfo.txt//;
  $sample_name =~ s/.hapInfo.txt//;
  push(@sample_list, $sample_name);
  push(@sample_list, $sample_name);
Line 41: Line 44:
  foreach my $hapString (keys(%{$hap_count_matrix{$probeID}->{$sample_name}})){
  foreach my $hapString (keys(%{$hap_count_matrix{$probeID}->{$sample_name}})){
  for(my $word_size = 1; $word_size<=length($hapString); $word_size++){
  for(my $word_size = 1; $word_size<=length($hapString); $word_size++){
  next if($word_size>5);
  next if($word_size>9);
  for(my $i=0; $i<=length($hapString)-$word_size; $i++){
  for(my $i=0; $i<=length($hapString)-$word_size; $i++){
  my $sub_hapString = substr($hapString,$i,$word_size);
  my $sub_hapString = substr($hapString,$i,$word_size);
Line 69: Line 72:
  }
  }
  print "Probe_id\t", join("\t", @sample_list), "\n";
  print "Probe_id\t", join("\t", @sample_list), "\n";
  foreach my $probeID (keys(%mch_load_matrix)){
  foreach my $probeID (sort keys(%mch_load_matrix)){
  print "$probeID";
  print "$probeID";
  foreach my $sample_name(@sample_list){
  foreach my $sample_name(sort @sample_list){
  $mch_load_matrix{$probeID}->{$sample_name}="NA" if(!$mch_load_matrix{$probeID}->{$sample_name});
  $mch_load_matrix{$probeID}->{$sample_name}="NA" if(! defined($mch_load_matrix{$probeID}->{$sample_name}));
  print "\t", $mch_load_matrix{$probeID}->{$sample_name};
  print "\t", $mch_load_matrix{$probeID}->{$sample_name};
  }
  }
  print "\n";
  print "\n";
  }
  }

Latest revision as of 21:55, 31 March 2016

#!/usr/bin/perl -w
# Hapinfo to methylation haplotype load (MHL)
# Run the script to the Hapinfo directory
# Contact: Shicheng Guo
# Version 1.3
# Update: 2016-03-31
use strict;
use Cwd;
my $usage = <<USAGE;
perl $0 Hapinfo_File_List.txt > MHL.output.txt
USAGE
die $usage if @ARGV <1;
my %mch_load_matrix;
my %probe_HMH_samples;
my %hap_count_matrix;
my $hapinfList=shift @ARGV;
open FF,$hapinfList;
chomp(my @hapInfo_files=<FF>);
close FF;
my @sample_list;
foreach my $hapInfo_file(@hapInfo_files){
my @line=split /\//,$hapInfo_file;
my $sample_name = $line[$#line];
$sample_name =~ s/.hapInfo.txt//;
push(@sample_list, $sample_name);
open(INFILE, "$hapInfo_file") || die("Error in opening $hapInfo_file!");
while(my $line = <INFILE>){
chop($line);
my @fields = split(/\t/, $line);
next if(scalar(@fields)<4);
my $probeID = $fields[0];
my $hapString = $fields[1];
next if(length($hapString)<1);
$hap_count_matrix{$probeID}->{$sample_name}->{$hapString}=$fields[2];
}
close(INFILE);
}
my @unmethylated_haps= ("T", "TT", "TTT", "TTTT", "TTTTT","TTTTTT","TTTTTTT","TTTTTTTT","TTTTTTTTT");
my @methylated_haps = ("C", "CC", "CCC", "CCCC", "CCCCC","CCCCCC","CCCCCCC","CCCCCCCC","CCCCCCCCC");
foreach my $probeID (keys(%hap_count_matrix)){
foreach my $sample_name (keys(%{$hap_count_matrix{$probeID}})){
my %k_mer_counts;
my $mc_hap_load=0;
foreach my $hapString (keys(%{$hap_count_matrix{$probeID}->{$sample_name}})){
for(my $word_size = 1; $word_size<=length($hapString); $word_size++){
next if($word_size>9);
for(my $i=0; $i<=length($hapString)-$word_size; $i++){
my $sub_hapString = substr($hapString,$i,$word_size);
next if($sub_hapString =~ /[NAG]/i);
$k_mer_counts{$word_size}->{$sub_hapString}+=$hap_count_matrix{$probeID}->{$sample_name}->{$hapString};
}
}
}
my $norm_factor=0;
foreach my $word_size (keys(%k_mer_counts)){
$k_mer_counts{$word_size}->{$unmethylated_haps[$word_size-1]}=0 if(!$k_mer_counts{$word_size}->{$unmethylated_haps[$word_size-1]});
$k_mer_counts{$word_size}->{$methylated_haps[$word_size-1]}=0 if(!$k_mer_counts{$word_size}->{$methylated_haps[$word_size-1]});
my $total_count=0;
foreach my $allele (keys(%{$k_mer_counts{$word_size}})){
$total_count+=$k_mer_counts{$word_size}->{$allele};
}
next if($total_count<1);
my $mh_fraction = $k_mer_counts{$word_size}->{$methylated_haps[$word_size-1]}/$total_count;
my $weight = $word_size;
$mc_hap_load += $weight*$mh_fraction;
$norm_factor+=$weight;
}
next if(!$norm_factor);
$mc_hap_load/=$norm_factor;
$mch_load_matrix{$probeID}->{$sample_name}=$mc_hap_load;
}
}
print "Probe_id\t", join("\t", @sample_list), "\n";
foreach my $probeID (sort keys(%mch_load_matrix)){
print "$probeID";
foreach my $sample_name(sort @sample_list){
$mch_load_matrix{$probeID}->{$sample_name}="NA" if(! defined($mch_load_matrix{$probeID}->{$sample_name}));
print "\t", $mch_load_matrix{$probeID}->{$sample_name};
}
print "\n";
}