BEGIN{FS=OFS="\t";
    getline;
    getline;
    if (ctl1==""){ctl1="PHR24r-IgG";}
    if (ctl2==""){ctl2="PHR19r-Input";}
    ns=NF;
    for (i=1;i<=ns;i++){
	sampID[i]=$i;
	if ($i==ctl1){ctl1i=i;print "Found ctr1"}
	if ($i==ctl2){ctl2i=i;print "Found ctr2"}	
    }
}{
    for (i=7;i<=ns;i++){
	x=$i;
	nr[i]+=x;fl[NR]=$6;
	r[NR,i]=x;l[NR]=$0;
    }
}END{
    nl=NR;
    print "Sample","Total_mapped_reads"
    ofn1[ctl1i]="RPKM_"ctl1"_vs_"ctl2;
    for (i=7;i<=ns;i++){
	print sampID[i],nr[i]
	rpm[i]=nr[i]/1000000.0;
	if ((i!=ctl1i)&&(i!=ctl2i)){
	    ofn1[i]="RPKM_"sampID[i]"_vs_"ctl1"_"ctl2;
#	    ofn2[i]="RPKM_"sampID[i]"_vs_"ctl2;
	}
    }
    print "control sample A is",ctl1
    print "control sample B is",ctl2

    for (i=7;i<=ns;i++){#headers
	if ((i!=ctl1i)&&(i!=ctl2i)){
	    printf("%s",sampID[1])>ofn1[i];
	    for (j=2;j<=ns;j++){
		printf("\t%s",sampID[j])>>ofn1[i];		
	    }
	    printf("\tRPKM_%s",sampID[i])>>ofn1[i];		
	    printf("\tRPKM_%s",sampID[ctl1i])>>ofn1[i];		
	    printf("\tfc_vs_%s",sampID[ctl1i])>>ofn1[i];		
	    printf("\tRPKM_%s",sampID[ctl2i])>>ofn1[i];		
	    printf("\tfc_vs_%s\n",sampID[ctl2i])>>ofn1[i];		
	}
    }
    i=ctl1i;
    printf("%s",sampID[1])>ofn1[i];
    for (j=2;j<=ns;j++){
	printf("\t%s",sampID[j])>>ofn1[i];		
    }
    printf("\tRPKM_%s",sampID[ctl1i])>>ofn1[i];		
    printf("\tRPKM_%s",sampID[ctl2i])>>ofn1[i];		
    printf("\tfc_vs_%s\n",sampID[ctl2i])>>ofn1[i];		

    for (p=3;p<=nl;p++){
	kb=fl[p]*0.001;
	for (i=7;i<=ns;i++){
	    rpkm[i]=(r[p,i]/rpm[i])/kb;
	}
	i=ctl1i;
	fc2=(1+rpkm[i])/(1+rpkm[ctl2i]);
	print l[p],rpkm[ctl1i],rpkm[ctl2i],fc2>>ofn1[i];
	
	for (i=7;i<=ns;i++){
	    if ((i!=ctl1i)&&(i!=ctl2i)){
		if (1){##		if (rpkm[i]>=1){
		    if (rpkm[ctl1i]==0){fc1="NA";}else{fc1=rpkm[i]/rpkm[ctl1i];}
		    if (rpkm[ctl2i]==0){fc2="NA";}else{fc2=rpkm[i]/rpkm[ctl2i];}
		    fc1=(1+rpkm[i])/(1+rpkm[ctl1i]);
		    fc2=(1+rpkm[i])/(1+rpkm[ctl2i]);
		    print l[p],rpkm[i],rpkm[ctl1i],fc1,rpkm[ctl2i],fc2>>ofn1[i];
		}
	    }
	}
    }    
 }