Perl ~/bin/hapinfo2mhl.pl: Difference between revisions
Jump to navigation
Jump to search
>Shicheng (Created page with " #!/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-02-29 us...") |
>Shicheng No edit summary |
||
(One intermediate revision by the same user not shown) | |||
Line 4: | Line 4: | ||
# Contact: Shicheng Guo | # Contact: Shicheng Guo | ||
# Version 1.3 | # Version 1.3 | ||
# Update: 2016- | # Update: 2016-03-31 | ||
use strict; | use strict; | ||
use Cwd; | use Cwd; | ||
my $usage = <<USAGE; | my $usage = <<USAGE; | ||
perl | 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; | ||
my $hapinfList=shift @ARGV; | |||
my @hapInfo_files= | 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 = $ | 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> | 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"; }