(*Mathematica:: Version 8.0      *)
(*:Context:: "mfit`"                     *)
(*:Title::  MFit            		*)
(*:Summary::   mfitdata           	*)
(*:References::                         
Yanai, I. Peshkin. L, Jorgenson, P. and Krischner, M.W.  (2011) Mapping gene expression in two Xenopus
     species: evolutionary constraints and developmental flexibility. Developmental Cell, 20(4), 483-496.
      PMID: 21497761. 
 *)
(*Date  2013 Stuart G. Baker   *)

(*Input
  mfitrawdata1.m 
  mfitrawdata2.m
  mfitgenename.m 
  *)

BeginPackage["mfitdata`"]  

GenDataX::usage="GenDataX[set], where  set = 1 or 2 refers to frog data"
GenTestData::usage="GenTestData[fulldataset,set]"

Clear[GenDataX,GenTestData]



Begin["Private`"]	  

	
    
GenDataX[datanum_,meanQ_:True]:=
Module[{data,geneid0,geneid,numtimes,dataset,dataset0,i,timex,
     datax0,datax,datax1,datax2,datax3,geneidx,genenamex,genename,
     xname,yname,plotname,
     datax1a,datax2a,datax3a,
     datax1b,datax2b,datax3b,
     datax1c,datax2c,datax3c,
     datanew,hz,name,name0,name1},
   (*READ IN DATA*)
   xname="time";
   yname="expression";
     If[datanum==1,
        plotname="X.laevis";
        datasetname="frog1";
         data = ReadList["mfitrawdata1.m", Number, RecordLists -> True]];
     If[datanum==2,
       plotname="X.tropicalis";
        datasetname="frog2";
       data = ReadList["mfitrawdata2.m", Number, RecordLists -> True]];
     (*TIMES*)
      timex = {1.5, 5, 7, 9, 13.25, 14.75, 16.25, 18.25, 19.75, 21.75, 24.75, 27.5, 35, 44.5}//N; 
      (*AVERATE OVER TECHNICAL REPLICATE*)   
      datax=Transpose[data];
       datax1=(datax[[Range[14]]] + datax[[Range[14]+14]] + datax[[Range[14]+28]])/3;
       datax2=(datax[[Range[14]+ 14 3]] + datax[[Range[14]+14 4]] + datax[[Range[14]+14 5]])/3;
        datax3=(datax[[Range[14]+ 14 6]] + datax[[Range[14]+14 7]] + datax[[Range[14]+14 8]])/3;
     (*average: already summed*)
    If[meanQ==False,
      datanew=Transpose@Join[datax1,datax2,datax3];  
       numsamples=3;
            times=Flatten@Table[timex,{i,1,numsamples}]];   
      (*MEAN*)
      If[meanQ==True,
       datanew=Transpose@(datax1+datax2+datax3)/3;  
           times=timex];
      
    (*geneid genename*)
       dataname3x = ReadList["mfitgenename.m", {Number, Number, String}, RecordLists -> True];
        dataname3=dataname3x[[1]];  
       {row,geneidx,genenamex}=Transpose[dataname3];   
       (*take top list of genes*); 
          numgenes=Length[datanew];
        genename=RemoveInitialBlankG@Take[genenamex,numgenes];
        geneid=Take[geneidx,numgenes];
        If[meanQ==False,
        datasetnamex=StringJoin["dataset",ToString[datanum]]];
        If[meanQ==True,
	        datasetnamex=StringJoin["datasetm",ToString[datanum]]];
        
       
        Print["DATASET ",datasetname, ": ",datasetnamex];
         dataset={datanew,geneid,genename,times,xname,yname,plotname,datasetname};
    Return[dataset]]



  
  

(*--------------------------test data-------------------------*)

GenTestData[dataset_,datanum_]:=
Module[{geneidvec},
type="B"; (*amount of data*)
geneidvec={20476, 6766, 3737,1474, 2012,18524,25498, 5731, 1320,21154,26635,17239,28040};
geneidvec={20476, 6766, 3737,1474, 2012,18524};
geneidvec={12080,6744,20755,20476};
geneidvec={23255,23254,28040,2012,18524};
geneidvec={2012,18524,24589,22315,11259};
geneidvec= {21061, 9439,20476, 6766, 3737,1474, 2012,18524,25498, 5731, 1320,21154,26635,17239,28040};
geneidvec= {21061, 6766,18052,26427,26223};

data=GenTestData[dataset,datanum,geneidvec];
Return[data]]

GenTestData[dataset_,datanum_,geneidvec_]:=
Module[{datanew,geneidx,genenamex,times,xname,yname,plotname,dataname,
           datanew0,geneidx0,genenamex0,name0,plotname0,posvec0,posvec},
 (*input all data*)          
    {datanew,geneidx,genenamex,times,xname,yname,plotname,dataname}=dataset;
(*name*)
  dataname0=StringJoin["test",type,ToString[datanum]];
 (*find position*)
    posvec0=Flatten@(Position[geneidx,#]&  /@ geneidvec);
   If[type=="A", posvec=posvec0];
   If[type=="B", posvec=posvec0 ~Union~ Range[20]];
 datasetnamex=StringJoin["dataset0",ToString[datanum]];
   Print["Test Data: ",datasetnamex];
   Print["genes ",geneidvec];
 (*take vectors with new position*)
   datanew0=datanew[[posvec]];
  geneidx0=geneidx[[posvec]];
  genenamex0=genenamex[[posvec]];
  plotname0=StringJoin["test",type,ToString[datanum]];
  dataset0={datanew0,geneidx0,genenamex0,times,xname,yname, plotname0,dataname0};
Return[dataset0]]


RemoveInitialBlankG[genename_]:=RemoveInitialBlankX[#]& /@ genename
 
RemoveInitialBlankX[genename_]:=
Module[{res,test,g2,g3},
g2=StringJoin[genename,"QQQ"];
g3=StringSplit[g2][[1]];
If[g3=="QQQ",test=True,test=False];
If[test==True, res="no name",
res=StringSplit[genename][[1]]];
Return[res]]


End[] 
EndPackage[]

datasetm1=GenDataX[1,True];
datasetm2=GenDataX[2,True];
dataset01=GenTestData[datasetm1,1];
dataset02=GenTestData[datasetm2,2];

