Ｒのhaplo.statsの使い方（３）Haploviewファイルから置換してハプロタイプ推定

• Haploviewの"xxx.ped" ファイルと"xxx.info"ファイルとがある場合に、そのファイルのあるディレクトリを指定し、haplo.statsの"haplo.em"関数の実行用ファイルを作成するperlのコード
• Ｒコンソール上で、上記ファイルを読み込み、haplo.emの結果を表示する関数とその用法
  
  • Haploviewファイルからhaplo.stats用テキストファイル作成ソース

   }elsif($genotype[$i][$numcol-1]==1){

    $line1 .= "1\t1";

   }elsif($genotype[$i][$numcol-1]==2){

    $line1 .= "0\t1";

   }elsif($genotype[$i][$numcol-1]==9){

    $line1 .= "9\t9";

   }
  my $flag=0;

  foreach $haplo (keys(%genohash)){

   if($haplo eq $line1){

    $genohash{$haplo}++;

    $flag=1;

    last;

   }

  }

  if($flag==0){

   $genohash{$line1}=1;

  }

 }

 my $cnt=0;

 foreach $haplo (keys(%genohash)){

  $weightlist[$cnt]=$genohash{$haplo};

  $weightedgenotypelist[$cnt]=$haplo;

  $cnt++;

 }


 return (\@weightlist,\@weightedgenotypelist);
}
sub printOUT($$){

 my ($list,$file)=@_;

 my $string="";

 for(my $i=0;$i<@$list;$i++){

  $string .= "$$list[$i]\n";

 }
 open (OUT, ">>$file");

 print OUT $string;

 close OUT;

}

• Ｒの関数と用法

#関数 dohstats と dohstats_w

#dohstatsは入力パターン1、dohstats_wは入力パターン2

#用法

# dohstats('hogeGeno.txt','hogeInfo.txt')

# dohstats_w('hogeGeno_w.txt','hogeInfo.txt','hogeWeight.txt')

dohstats<-function(file1,file2){

library(haplo.stats)
geno<-read.table(file=file1,FALSE,"\t")

label1<-read.table(file=file2,FALSE,"\t")

save.em.keep<-haplo.em(geno,locus.label=label1$V1,miss.val=9)

print.haplo.em(save.em.keep)

}
dohstats_w<-function(file1,file2,file3){

library(haplo.stats)
geno<-read.table(file=file1,FALSE,"\t")

label1<-read.table(file=file2,FALSE,"\t")

weight<-read.table(file=file3,FALSE,"\t")
save.em.keep<-haplo.em(geno,locus.label=label1$V1,weight=weight$V1,miss.val=9)

print.haplo.em(save.em.keep)

}