#!	/usr/bin/perl	-w

=head1	header

	###########################################################
	#  
	#  This script cut SNPs from reference haplotype data based on a flag file
	#  
	#  perl cut1000GHapwithflag.pl Chr22_Ref_Ex Chr22_Ref_Ex_flag.txt
	#  
	#  Input : Chr22_Ref_Ex.hap, Chr22_Ref_Ex.map, Chr22_Ref_Ex.snps
	#  Output : Chr22_Ref_Ex_cut.hap, Chr22_Ref_Ex_cut.map, Chr22_Ref_Ex_cut.snps
	#  
	#  1KG haplotype data is available from http://www.sph.umich.edu/csg/abecasis/MACH/download/
	#  
	#  Any questions to Yukinori Okada (http://plaza.umin.ac.jp/~yokada/datasource/software.htm   yokada@broadinstitute.org)

=cut



&main(@ARGV);
exit;

=head1	main

=cut

sub	main(@){
	my ($input_head,$flag_file)=@_;
	my $ped = $input_head.".hap";
	my $info = $input_head.".snps";
	my $map = $input_head.".map";
	my $cutped = $input_head."_cut.hap";
	my $cutinfo = $input_head."_cut.snps";
	my $cutmap = $input_head."_cut.map";
	my $keepflag = 1;


	open (INPUTflag, "$flag_file");
	my @flag;
	my $row=0;
	while(<INPUTflag>){
		chomp;
		$flag[$row]=$_;
		$row++;
	}

	&linecut($info, $cutinfo, @flag);
	&linecut($map, $cutmap, @flag);
#	for (my $i=0;$i<6;$i++) {
#		unshift @flag,$keepflag;
#	}
	&rowcut($ped, $cutped, @flag);

	sub linecut{
		my ($input_file, $outputfile, @flag)=@_;
		my $counter=0;
		open (INPUT, "$input_file");
		open (OUT, "> $outputfile");
		while(<INPUT>){
#			chomp;
			if ($flag[$counter]==$keepflag) {
				print OUT $_;
			}
			$counter++;
		}
		close INPUT;
		close OUT;
	}

	sub rowcut{
		my ($input_file, $outputfile, @flag)=@_;
		my $counter=0;
		open (INPUT, "$input_file");
		open (OUT, "> $outputfile");
		while(<INPUT>){
			chomp;
#			my @inline = split(/\t/);
#			my @inline2 = split(//, $inline[1]);
			my @inline = split(/ /);
			my @inline2 = split(//, $inline[2]);
			my $out = $inline[0]." ".$inline[1]." ";
			for(my $i=0;$i<@inline2;$i++){
				if ($flag[$i]==$keepflag) {
					$out .= "$inline2[$i]";
				}
			}
			print OUT "$out\n";
		}
		close INPUT;
		close OUT;
	}

}
