#!/usr/bin/perl
=head2 SYNOPSIS

	Description: _focus script is used for process the categorized types by assigning them into a certain subgenome.
             
=head2 AUTHOR

	Kang Zhang  Email:  demut@foxmail.com
	
=cut

#Input
#open SAM,"<$ARGV[0]" or die;
#prefix of the output
$DataType=$ARGV[1];
if(!defined $DataType){
	die "Pls provide the data type! Could be TypeI TypeII...TypeVII\n";
}
$assignGenome=$ARGV[2];
if(!defined $assignGenome){
	die "Pls determine the genome to assign! Could be 1 (first genome), 2 (second genome)\n";
}
#Parameters:
## 归属的标记
$markContig="_JZS"; #as subgenome 2, not match as subgenome 1.


while(<STDIN>){
	chomp;
	if(/^@/){print "$_\n";next;}
	my @arr=split;
	if( ($read ne $arr[0]) and (defined $read)){
		my $length=keys %count; #print STDERR "#$length\n";
		if($length > 100000){
			$countK++;
			print STDERR "Processing $countK"."00k records...";
			&analyzeMultiple;
			print STDERR "Finished\n";
			%markMulti=();%mapInfo=(); %count=(); %whichEnd=(); %flagInfo=(); %mapContigInfo=(); %mapQua=(); %locate=(); %locateCount=(); %multiple=();
			%CoorSub=(); %countCoor=(); %CoorCount=(); %data=(); %appendData=(); undef @orderArray; %occur=(); %processStatus=(); %tmpTargetSub=();	%priMark=();
		}
	}
	$read=shift @arr;
	my $flag=shift @arr;
	$tmpFlag=$flag;
	if($flag >= 512){next;undef $Primary;}
	if($flag >= 256){
		$Primary=0;
		$markMulti{$read}=1;
		$tmpFlag-=256;
	}else{$Primary=1;}
	my $flagCons=&SubFlagTable($tmpFlag);
	#read unmapped? contain 4 is unmapped.
	if($flagCons =~ /\+4/){$mapped=0;}else{$mapped=1;} 
	#first or second read in pair
	if($flagCons =~ /64/){$whichend=1;}elsif($flagCons =~/128/){$whichend=2;}else{print STDERR "Unpaired read found!  $_\n";next;}
	$mapInfo{$read}{$whichend}=$mapped;
	if($mapped==0){next;}
	$count{$read}++;

	#获取first or second 的信息
	$whichEnd{$read}{$count{$read}}=$whichend;
	#获取read 的flag信息
	$flagInfo{$read}{$count{$read}}=$flag;
	#获取比对位置contig信息
	my $tmpContigMap=$arr[0]."\t$arr[1]";
	my $tmpContig=shift @arr;shift @arr; 
	$mapContigInfo{$read}{$count{$read}}=$tmpContigMap;
	#获取质量信息
	$mapQua{$read}{$count{$read}}=shift @arr;
	#获取比对的亚基因组信息
	if($tmpContig =~/$markContig$/){
		$subgenome=2;
	}else{$subgenome=1;} #print STDERR "$read\t$flag\t$tmpContig\t$mapContigInfo{$read}{$count{$read}}\t$mapQua{$read}{$count{$read}}\t---$subgenome\n";
	$locate{$read}{$whichend}{$subgenome}=1;
	$locateCount{$read}{$count{$read}}=$subgenome;
	if($Primary==0){
		$multiple{$read}{$whichend}=1;
	}
	
	#获取coordinate 信息
	if($flagCons =~ /\+2/){$coordinate=1;	$CoorSub{$read}{$subgenome}++;}#统计在哪些subgenome中有coordinate
	else{$coordinate=0;}
	$countCoor{$read}+=$coordinate;
	$CoorCount{$read}{$count{$read}}=$coordinate;
	
	#获取records：
	$data{$read}{$count{$read}}=$_;
	my $tmpStr=join ("\t",@arr);
	$appendData{$read}{$count{$read}}=$tmpStr;
	if(!defined $occur{$read}){
		push (@orderArray,$read);
		$occur{$read}=1;
	}
	
}

&analyzeMultiple;

sub analyzeMultiple {
if($DataType eq "TypeII"){
	foreach my $read (keys %markMulti){
	my $attribute = $assignGenome;
	foreach my $num (keys %{$locateCount{$read}}){
		if($locateCount{$read}{$num} == $attribute){
			$processStatus{$read}{$num}="Stay-TypeII"; #Flag 设置其中一个为primary，其他为multiple 如B9007CACXX:8:1101:5642:2391#ACTTGAAT
			##如果分配后，每一个end只剩下一个且为primary，将其MQ设置为60
		}else{
			$processStatus{$read}{$num}="Dump-TypeII";
		}
	}
	&processStayPrint($read);
	}
}elsif($DataType eq "TypeV" ){
	foreach my $read (keys %markMulti){
	my @tmpCountCoor=keys %{$CoorSub{$read}};
	if($#tmpCountCoor == 0){
		#多个coordinate 在同一个subgenome中
		$tmpTargetSub{$read}=$tmpCountCoor[0];
	}else{
		#在不同的subgenome,要决定归属
		my $attribute = $assignGenome;
		$tmpTargetSub{$read}=$attribute;
	}
	#这样$tmpTargetSub{$read} 就存储了subgenome 的选择信息
	foreach my $num (keys %{$locateCount{$read}}){
		#仅属于选择的subgenome的，且仅保留coordinate的
		if( ($locateCount{$read}{$num} == $tmpTargetSub{$read}) and ($CoorCount{$read}{$num} == 1)){
			$processStatus{$read}{$num}="Stay-TypeV"; #Flag 如果是non-primary，减去-256，MQ 改为60 (本来就有且只有一个保留下来)
		}else{
			$processStatus{$read}{$num}="Dump-TypeV";
		}
	}
	&processStayPrint($read);
	}
}elsif($DataType eq "TypeVI"){
	foreach my $read (keys %markMulti){
	my @tmplocate1=keys %{$locate{$read}{1}}; my @tmplocate2=keys %{$locate{$read}{2}};
	#都不是coordinate,且在不止一个subgenome中有multiple 比对，比如first end 在两个subgenome中有multiple 比对，second在一个subgenome中有multiple
	#或者，两个end 都在两个subgenome中有multiple 比对
	#如果是前者：即其中一个end是在亚基因组内部多比对，那么就不选择，直接归类于这个亚基因组
	if( ($#tmplocate1 + $#tmplocate2) == 1){ #肯定一个是0，一个是1，即一个end在subgenome 内部多比对，一个在两个subgenome 多比对
		#确定哪一个是内部多比对
		if($#tmplocate1 == 0){
			$tmpTargetSub{$read}=$tmplocate1[0];
		}else{
			$tmpTargetSub{$read}=$tmplocate2[0];
		}						
	}else{
		#后者，两个end，两个subgenome
		#如果$#tmplocate1 + $#tmplocate2 >1 ,进行选择，将属于另外一个subgenome 的比对全部去除，选择的这个subgenome则不处理
		$tmpTargetSub{$read} = $assignGenome;
	}
	foreach my $num (keys %{$locateCount{$read}}){
		if( $locateCount{$read}{$num} == $tmpTargetSub{$read}){
			$processStatus{$read}{$num}="Stay-TypeVI"; #Flag：对于first 和 second，各设置其中一个为primary，其他为multiple，如B9007CACXX:8:1101:20240:2433#ACTTGAAT
			#或者 B9007CACXX:8:1101:18177:2513#ACTTGAAT
			#如果分配后，每一个end只剩下一个且为primary，将其MQ设置为60
		}else{
			$processStatus{$read}{$num}="Dump-TypeVI";
		}
	}
	&processStayPrint($read);
	}
}else{
	print STDERR "We don't have this data type: $DataType\n";
	last;
}

}


sub processStayPrint(){
	my $read=$_[0];
undef @printString; undef @StringOrder;
	if( ($processStatus{$read}{1} =~ /TypeII$/) or ($processStatus{$read}{1} =~ /TypeVI$/)){
			%ReadRetain=();
			for(my $i=1;$i<=$count{$read};$i++){
				#print "->$data{$read}{$i}\t$processStatus{$read}{$i}\n";
				#if ($processStatus{$read}{$i} =~/Dump/){print ERR "$processStatus{$read}{$i}\t$data{$read}{$i}\n";}
				next if ($processStatus{$read}{$i} =~/Dump/);
				my $tmpEnd=$whichEnd{$read}{$i};
				if($flagInfo{$read}{$i} < 256){
					#this one is primary
					$priMark{$read}{$tmpEnd}=$i;
				}
				#也有可能保留的records中都没有primary，就要从中选择一个减去256，作为primary
				push (@{$ReadRetain{$tmpEnd}},$i);
			}
			for(my $n=1;$n<=2;$n++){
				if(! defined $priMark{$read}{$n}){
					#保留的没有一个primary,将第一个记录改为primary：
					$flagInfo{$read}{$ReadRetain{$n}[0]}-=256;
					$priMark{$read}{$n}=$ReadRetain{$n}[0];	
				}
				if($#{$ReadRetain{$n}} == 0){
					#如果这个end 只有这个record，就将其MQ设置为60：
					$mapQua{$read}{$ReadRetain{$n}[0]}=60;
				}
				foreach my $num (@{$ReadRetain{$n}}){
					my $tmpString ="$read\t$flagInfo{$read}{$num}\t$mapContigInfo{$read}{$num}\t$mapQua{$read}{$num}\t$appendData{$read}{$num}";
					push (@printString,$tmpString);
					push (@StringOrder,$num);
				}
			}
	}elsif( ($processStatus{$read}{1} =~ /TypeIV$/) or ($processStatus{$read}{1} =~ /TypeV$/)){
		for(my $i=1;$i<=$count{$read};$i++){
				#print "->$data{$read}{$i}\t$processStatus{$read}{$i}\n";
				#if ($processStatus{$read}{$i} =~/Dump/){print ERR "$processStatus{$read}{$i}\t$data{$read}{$i}\n";}
				next if ($processStatus{$read}{$i} =~/Dump/);
				if($flagInfo{$read}{$i} > 256){
					$flagInfo{$read}{$i}-=256;
				}
				#print "$read\t$flagInfo{$read}{$i}\t$mapContigInfo{$read}{$i}\t60\t$appendData{$read}{$i}\n";
				my $tmpString ="$read\t$flagInfo{$read}{$i}\t$mapContigInfo{$read}{$i}\t60\t$appendData{$read}{$i}";
				push (@printString,$tmpString);
				push (@StringOrder,$i);
			}
	}
	&ModifyNHfieldAndPrint($read);
}


##############################################
sub ModifyNHfieldAndPrint {
	my $read=$_[0];
	#其他信息采用全局变量，以避免重复赋值，拖慢速度
	foreach my $num (@StringOrder){
		#先获取处理之后剩余的records数目，每一个end单独分析：
		my $tmpEnd=$whichEnd{$read}{$num};
		$endCount{$read}{$tmpEnd}++;
	}
	for(my $i=0;$i<=$#StringOrder;$i++){
		my $tmpEnd=$whichEnd{$read}{$StringOrder[$i]};
		my $tmpNH= $endCount{$read}{$tmpEnd};
		#这个end 剩余几个records，NH就是几;如果没有NH field，自然不改动
		$printString[$i] =~ s/NH\:i\:(\d+)/NH\:i\:$tmpNH/ if($printString[$i]=~/NH\:i\:/);
		print "$printString[$i]\n";
	}
	
}


sub SubFlagTable
{
	my %TABLE = qw(0   0
      1   1
      2   2
      3   1+2
      4   0+4
      5   1+4
      6   0+2+4
      7   1+2+4
      8   0+8
      9   1+8
      10   0+2+8
      11   1+2+8
      12   0+4+8
      13   1+4+8
      14   0+2+4+8
      15   1+2+4+8
      16   0+16
      17   1+16
      18   0+2+16
      19   1+2+16
      20   0+4+16
      21   1+4+16
      22   0+2+4+16
      23   1+2+4+16
      24   0+8+16
      25   1+8+16
      26   0+2+8+16
      27   1+2+8+16
      28   0+4+8+16
      29   1+4+8+16
      30   0+2+4+8+16
      31   1+2+4+8+16
      32   0+32
      33   1+32
      34   0+2+32
      35   1+2+32
      36   0+4+32
      37   1+4+32
      38   0+2+4+32
      39   1+2+4+32
      40   0+8+32
      41   1+8+32
      42   0+2+8+32
      43   1+2+8+32
      44   0+4+8+32
      45   1+4+8+32
      46   0+2+4+8+32
      47   1+2+4+8+32
      48   0+16+32
      49   1+16+32
      50   0+2+16+32
      51   1+2+16+32
      52   0+4+16+32
      53   1+4+16+32
      54   0+2+4+16+32
      55   1+2+4+16+32
      56   0+8+16+32
      57   1+8+16+32
      58   0+2+8+16+32
      59   1+2+8+16+32
      60   0+4+8+16+32
      61   1+4+8+16+32
      62   0+2+4+8+16+32
      63   1+2+4+8+16+32
      64   0+64
      65   1+64
      66   0+2+64
      67   1+2+64
      68   0+4+64
      69   1+4+64
      70   0+2+4+64
      71   1+2+4+64
      72   0+8+64
      73   1+8+64
      74   0+2+8+64
      75   1+2+8+64
      76   0+4+8+64
      77   1+4+8+64
      78   0+2+4+8+64
      79   1+2+4+8+64
      80   0+16+64
      81   1+16+64
      82   0+2+16+64
      83   1+2+16+64
      84   0+4+16+64
      85   1+4+16+64
      86   0+2+4+16+64
      87   1+2+4+16+64
      88   0+8+16+64
      89   1+8+16+64
      90   0+2+8+16+64
      91   1+2+8+16+64
      92   0+4+8+16+64
      93   1+4+8+16+64
      94   0+2+4+8+16+64
      95   1+2+4+8+16+64
      96   0+32+64
      97   1+32+64
      98   0+2+32+64
      99   1+2+32+64
      100   0+4+32+64
      101   1+4+32+64
      102   0+2+4+32+64
      103   1+2+4+32+64
      104   0+8+32+64
      105   1+8+32+64
      106   0+2+8+32+64
      107   1+2+8+32+64
      108   0+4+8+32+64
      109   1+4+8+32+64
      110   0+2+4+8+32+64
      111   1+2+4+8+32+64
      112   0+16+32+64
      113   1+16+32+64
      114   0+2+16+32+64
      115   1+2+16+32+64
      116   0+4+16+32+64
      117   1+4+16+32+64
      118   0+2+4+16+32+64
      119   1+2+4+16+32+64
      120   0+8+16+32+64
      121   1+8+16+32+64
      122   0+2+8+16+32+64
      123   1+2+8+16+32+64
      124   0+4+8+16+32+64
      125   1+4+8+16+32+64
      126   0+2+4+8+16+32+64
      127   1+2+4+8+16+32+64
      128   0+128
      129   1+128
      130   0+2+128
      131   1+2+128
      132   0+4+128
      133   1+4+128
      134   0+2+4+128
      135   1+2+4+128
      136   0+8+128
      137   1+8+128
      138   0+2+8+128
      139   1+2+8+128
      140   0+4+8+128
      141   1+4+8+128
      142   0+2+4+8+128
      143   1+2+4+8+128
      144   0+16+128
      145   1+16+128
      146   0+2+16+128
      147   1+2+16+128
      148   0+4+16+128
      149   1+4+16+128
      150   0+2+4+16+128
      151   1+2+4+16+128
      152   0+8+16+128
      153   1+8+16+128
      154   0+2+8+16+128
      155   1+2+8+16+128
      156   0+4+8+16+128
      157   1+4+8+16+128
      158   0+2+4+8+16+128
      159   1+2+4+8+16+128
      160   0+32+128
      161   1+32+128
      162   0+2+32+128
      163   1+2+32+128
      164   0+4+32+128
      165   1+4+32+128
      166   0+2+4+32+128
      167   1+2+4+32+128
      168   0+8+32+128
      169   1+8+32+128
      170   0+2+8+32+128
      171   1+2+8+32+128
      172   0+4+8+32+128
      173   1+4+8+32+128
      174   0+2+4+8+32+128
      175   1+2+4+8+32+128
      176   0+16+32+128
      177   1+16+32+128
      178   0+2+16+32+128
      179   1+2+16+32+128
      180   0+4+16+32+128
      181   1+4+16+32+128
      182   0+2+4+16+32+128
      183   1+2+4+16+32+128
      184   0+8+16+32+128
      185   1+8+16+32+128
      186   0+2+8+16+32+128
      187   1+2+8+16+32+128
      188   0+4+8+16+32+128
      189   1+4+8+16+32+128
      190   0+2+4+8+16+32+128
      191   1+2+4+8+16+32+128
      192   0+64+128
      193   1+64+128
      194   0+2+64+128
      195   1+2+64+128
      196   0+4+64+128
      197   1+4+64+128
      198   0+2+4+64+128
      199   1+2+4+64+128
      200   0+8+64+128
      201   1+8+64+128
      202   0+2+8+64+128
      203   1+2+8+64+128
      204   0+4+8+64+128
      205   1+4+8+64+128
      206   0+2+4+8+64+128
      207   1+2+4+8+64+128
      208   0+16+64+128
      209   1+16+64+128
      210   0+2+16+64+128
      211   1+2+16+64+128
      212   0+4+16+64+128
      213   1+4+16+64+128
      214   0+2+4+16+64+128
      215   1+2+4+16+64+128
      216   0+8+16+64+128
      217   1+8+16+64+128
      218   0+2+8+16+64+128
      219   1+2+8+16+64+128
      220   0+4+8+16+64+128
      221   1+4+8+16+64+128
      222   0+2+4+8+16+64+128
      223   1+2+4+8+16+64+128
      224   0+32+64+128
      225   1+32+64+128
      226   0+2+32+64+128
      227   1+2+32+64+128
      228   0+4+32+64+128
      229   1+4+32+64+128
      230   0+2+4+32+64+128
      231   1+2+4+32+64+128
      232   0+8+32+64+128
      233   1+8+32+64+128
      234   0+2+8+32+64+128
      235   1+2+8+32+64+128
      236   0+4+8+32+64+128
      237   1+4+8+32+64+128
      238   0+2+4+8+32+64+128
      239   1+2+4+8+32+64+128
      240   0+16+32+64+128
      241   1+16+32+64+128
      242   0+2+16+32+64+128
      243   1+2+16+32+64+128
      244   0+4+16+32+64+128
      245   1+4+16+32+64+128
      246   0+2+4+16+32+64+128
      247   1+2+4+16+32+64+128
      248   0+8+16+32+64+128
      249   1+8+16+32+64+128
      250   0+2+8+16+32+64+128
      251   1+2+8+16+32+64+128
      252   0+4+8+16+32+64+128
      253   1+4+8+16+32+64+128
      254   0+2+4+8+16+32+64+128
      255   1+2+4+8+16+32+64+128);
	  my $subInput=$_[0];
	  my $subOutput;
	  $subOutput=$TABLE{$subInput};
	  return($subOutput);
}
