Chr3Dir=$HomeDirectory<>"\\Desktop\\AfterDisert.Work\\ChrThreeFullHaplosFiles\\";
PopulFiles=FileNames["*",Chr3Dir];
NewSNPsDir=$HomeDirectory<>"\\Desktop\\AfterDisert.Work\\Chr3SNPsPotentials\\";

PopulSnps[Populno_]:=(
HaplocationFiles=FileNames["Haplocation*",PopulFiles[[Populno]]];
OnePopulSNPs[Populno]={};
For[Snp=1,Snp≤Length[HaplocationFiles],Snp++,
ImSNPs=Import[HaplocationFiles[[Snp]],"Table"];
DroImSNPs=Drop[ImSNPs,1];
OnePopulSNPsLocations=Table[DroImSNPs[[loc,2]],{loc,1,Length[DroImSNPs]}];
AppendTo[OnePopulSNPs[Populno],OnePopulSNPsLocations];
];
Flatten[OnePopulSNPs[Populno]]
);
AllSNPS=Table[PopulSnps[Populno],{Populno,1,Length[PopulFiles]}];
RepetedSNPs=Intersection[PopulSnps[1],PopulSnps[2],PopulSnps[3],PopulSnps[4],PopulSnps[5],PopulSnps[6],PopulSnps[7],PopulSnps[8]];

RsBlockfilename[popno_]:=FileNames["Haploblock*",FileNames["*",Chr3Dir][[popno]]];
RsLocationfilename[popno_]:=FileNames["Haplocation*",FileNames["*",Chr3Dir][[popno]]];
popslashposition[popno_]:=StringPosition[FileNames["*",Chr3Dir][[popno]],"\\"][[-1,-1]];
popname[popno_]:=StringDrop[FileNames["*",Chr3Dir][[popno]],popslashposition[popno]];
SNPname[snpno_]:=RepetedSNPs[[snpno]];
AllpopTEs={{"CLM",1.1255766729746808`},{"FIN",1.1255515144314452`},{"IBS",1.1107920563525304`},{"JPT",1.1088685528507198`},{"KHV",1.1065608191416791`},{"MSL",1.2377165220900395`},{"PEL",1.1087444589495752`},{"YRI",1.1959850154565101`},
{"LWK",1.2251228197944872`,1.7`},{"TSI",1.122627188367992`,1.35`},{"CHB",1.1088685528507198`,1.31`}};

(*SNPLocator has output imported={location file no, location data, block data}*)
SNPLocator[snpno_,popno_]:=(
pos={};
jj=0;
While[pos=={},
jj=jj+1;
locdata=Import[RsLocationfilename[popno][[jj]],"Table"];
pos=Position[locdata,ToString[RepetedSNPs[[snpno]]]];
];
{jj,locdata,Import[StringReplace[RsLocationfilename[popno][[jj]],{"Haplocation"->"Haploblock"}],"Table"]}
);
(*SNPelem takes imported SNPLocator[snpno,popno] *)
SNPelem[snpno_,imported_]:= (
SNPLocationdata=imported;
SNPpos=Position[SNPLocationdata[[2]],ToString[RepetedSNPs[[snpno]]]][[1,1]];
SNPLocationdata[[2]][[SNPpos]]
);
(*SNPinfo extracts information from SNPelem[snpno,SNPLocator[snpno,popno]]*)
SNPinfo[snpelem_]:= (
{{snpelem[[1]],snpelem[[2]],snpelem[[3]]},{snpelem[[-1]],snpelem[[-2]]}}
);
SNPInfoList[snpno_,imported_]:=(
Infolist=SNPinfo[SNPelem[snpno,imported]];
blockdata=imported[[3]];
Markerpos=Position[blockdata,"MARKERS:"];
Markerblock=(
jj=1;
While[(jj<=Length[Markerpos])&&(Position[blockdata[[Markerpos[[jj,1]]]],Infolist[[1,1]]]== {}), 
jj++
];
jj
);
(* WE NEED TO BYPASS WHAT FOLLOWS IF SNP NOT IN LD Markerblock>Length[Markerpos] *)
LDlist:=(
MarkerblockReal=1.0*Markerblock;
Actualblkline=(
kk=1;
blkno=Position[blockdata,MarkerblockReal];
While[(kk≤Length[blkno])&&(blockdata[[blkno[[kk,1]],1]]≠"BLOCK"),
kk++
];
blkno[[kk,1]]
);
Finaldatalineno=Which[
Length[Markerpos]>Position[Markerpos,Actualblkline][[1,1]],Markerpos[[Position[Markerpos,Actualblkline][[1,1]]+1,1]]-2,
True,Length[blockdata]
];

Markersno=Drop[blockdata[[Actualblkline]],3];
Haplist={};
For[hn=Actualblkline+1,hn<=Finaldatalineno,hn++,
AppendTo[Haplist, {ToString[blockdata[[hn,1]]],ToExpression[blockdata[[hn,2]]]}];
];
{Infolist,Markersno,Haplist}
);
Which[Markerblock>Length[Markerpos],{Infolist},Markerblock<=Length[Markerpos],LDlist]
);
Print[" Defining Entropy"]
NonLinkEntropy[infoln_]:=(
maf=infoln[[1,2,2]];
Which[maf==0,0,
True,- maf*Log[2,maf]-(1-maf)*Log[2,(1-maf)]
]
);
LDEntropy[infoln_]:=(
(* infoln[[3,n,2]] is the nth probability *)
Sum[- infoln[[3,n,2]]*Log[2,infoln[[3,n,2]]],{n,1,Length[infoln[[3]]]}]
);
BlockEntropy[infoln_]:= Which[
Length[infoln]==1,NonLinkEntropy[infoln],
True, LDEntropy[infoln]
];
Print[" Defining Genodynamic Potentials "];
(* SNP Potential*)
NonLinkSNPpotential[infoln_,TE_]:=Which[
infoln[[1,2,2]]==0,1-TE,
True,1-TE+TE*NonLinkEntropy[infoln]
];
(*Block Potential*)
LDblockpotential[infoln_,TE_]:=(
numSNPs=StringLength[infoln[[3,1,1]]];(* same for all Haplotypes *)
(1-TE)*numSNPs+TE*LDEntropy[infoln]
);
Blockpotential[infoln_,TE_]:= Which[
Length[infoln]==1,NonLinkSNPpotential[infoln,TE],
True,LDblockpotential[infoln,TE]
];
(* Allelic Potential*)
NonLinkMu[infoln_,TE_,alno_]:=(
maf=infoln[[1,2,2]];
Which[
alno==2,(1-TE)-TE *Log[2,maf],
alno==1,(1-TE)-TE *Log[2,(1-maf)]
]
);
(* Haplotype Potential*)
LDMu[infoln_,TE_,Hnum_]:=(
numSNPs=StringLength[infoln[[3,1,1]]];
(1-TE)*numSNPs-TE*Log[2,infoln[[3,Hnum,2]]]
);
(*MAFs[snpno_,popno_] gives outcome {maf of snpno,sum mafs in the block} *)
MAFs[snpno_,imported_]:=(
locdata=imported[[2]];
snpinflst=SNPInfoList[snpno,imported];
markernums=snpinflst[[2]];
snpmaf=snpinflst[[1,2,2]];
locmarkers=Drop[locdata,1][[markernums,1]];
Which[
StringLength[snpinflst[[3,1,1]]]!=Length[markernums],Print["Haploblock Length is not same as number of markers "],
locmarkers==markernums,{snpmaf,Tr[Drop[locdata,1][[markernums,-2]]],Length[markernums]},
True,Print["Error in Marker Locations"]
]
);
(*DistributedSNPPotenial gives the distributed potential of snpno in a Haploblock*)
DistributedSNPPotenial[snpno_,imported_,TE_]:=(
mafvalue=MAFs[snpno,imported];
snplist=SNPInfoList[snpno,imported];
(1-TE)+(Blockpotential[snplist,TE]-mafvalue[[3]](1-TE))*(mafvalue[[1]]/mafvalue[[2]])
);
(*BindingEnergy = DistributedSNPPotenial-NonLinkSNPpotential*)
BindingEnergy[snpno_,imported_,TE_]:= (
snplist=SNPInfoList[snpno,imported];
DistributedSNPPotenial[snpno,imported,TE]- NonLinkSNPpotential[snplist,TE]
);
(*DistributedAllelicpotential[snpno,popno,TE,allele string] gives NonLinkMu+BindingEnergy*)
DistributedAllelicpotential[snpno_,imported_,TE_,al_]:=(
snplist=SNPInfoList[snpno,imported];
alleleno=Which[
StringPosition[snplist[[1,2,1]],al][[1,1]]==3,2,
StringPosition[snplist[[1,2,1]],al][[1,1]]==1,1,
True, Print[" Not an allele in this SNP "]
];
nldMu=NonLinkMu[snplist,TE,alleleno];
benergy=BindingEnergy[snpno,imported,TE];
nldMu+benergy
);
(*SinglePopSNPInformation[snpno,imported] has output {{popno,chrno,rs},{{TE,1000G},{TE,1000G}},{{a1:a2,{SNP pot, a1 pot, 2 pot, binding E}},{{a1:a2,{SNP pot, a1 pot, 2 pot, binding E}}}} *)
SinglePopSNPInformation[snpno_,popno_,imported_]:=(
Inflst=SNPInfoList[snpno,imported];
TE=AllpopTEs[[Position[AllpopTEs,popname[popno]][[1,1]]]][[2]];
InfoPots=Which[Length[Inflst]≤1,{NonLinkSNPpotential[Inflst,TE],NonLinkMu[Inflst,TE,1],NonLinkMu[Inflst,TE,2]},
True,
{DistributedSNPPotenial[snpno,imported,TE],DistributedAllelicpotential[snpno,imported,TE,StringDrop[Inflst[[1,2,1]],-2]],
(DistributedAllelicpotential[snpno,imported,TE,StringDrop[Inflst[[1,2,1]],2]]),BindingEnergy[snpno,imported,TE]}
];
{{popname[popno],"Chrno",Inflst[[1,1,2]]},{{TE,"1000Genome"},{TE,"1000Genome"}},{{Inflst[[1,2,1]],InfoPots},{Inflst[[1,2,1]],InfoPots}}}
);
(*SNPInformation[snpno] gives SinglePopSNPInformation for all populations*)
SNPInformation[snpno_]:=(
InformTable=Table[SNPLocator[snpno,popno],{popno,1,Length[PopulFiles]}];
Table[SinglePopSNPInformation[snpno,popno,InformTable[[popno]]],{popno,1,Length[PopulFiles]}]
);
(*NewInfoSNPslist[snpno] saves SNPInformation into NewSNPsDir<>SNPname[snpno]*)
(*It takes 12 min for 10 SNPs*)
NewInfoSNPslist[snpno_]:=(
SNPname[snpno];
NewInfolist=SNPInformation[snpno];
Put[NewInfolist,NewSNPsDir<>SNPname[snpno]]
);

(*Every time add one to the previous RunNo*)
(*As there is 375769 SNPs, the last run is 375000 to 375769 - Do it by hand-*)
RunNo=106;
NoSNPs=1000;
Print["Running SNP no ", 1+RunNo*NoSNPs," to ", (RunNo+1)*NoSNPs ];
Print["start date =",strt=Date[]];
For[snpno=(NoSNPs*RunNo)+1,snpno≤NoSNPs*(RunNo+1),snpno++,NewInfoSNPslist[snpno]];
Print["Finish date =",fin=Date[],". Time Elapsed=",fin-strt ];