namespace ddct
{
	class MainClass
	{
		static double estimate_efficiency(Double cp1, Double cp2, Double dil1, Double dil2)
		{
			double ratio=dil1/dil2;
			double shift=cp2-cp1;
			double copyratio=Math.Pow(ratio,1/shift);
			return copyratio;
		}
		static double normalize_cp(Double cp, Double copy_ratio)
		{
			return cp*Math.Log(copy_ratio)/Math.Log(2);
		}
		public static void Main(string[] args)
		{
			/**
			 * Lesson 10: Cleaning up the data by removing everything that is no number
			 * and averaging technical replicates for each group.
			 */
 	  	   Table table = new Table("qpcr.tsv");
           Table cleanedup = new Table();
           List<Table> groups = table.groups("Plate", "CellType", "CellLine", "Gene", "Dilution");
           foreach (Table group in groups)
           {
               Double average = 0;
               int count = 0;
               foreach (Record rec in group)
               {
                   try
                   {
                       String cp_as_text = (String)rec["CP"];
                       Double cp_as_double = Convert.ToDouble(cp_as_text);
                       average += cp_as_double;
                       count++;
                      }
                   catch (FormatException)
                   {
                   }
               }
				if (count>0)
				{               average /= count;
               Record record = group.key.copy();
               record["AverageCP"] = average;
               cleanedup.Add(record);
               }
			}

			/**
			 * Lesson 11: Estimating the efficiencies for each gene 
			 * estimate the efficience per plate, per probe, per celllines, celltype, accross dilutions 
			 * (and of course accross averages). For each combination of dilutions we estimate the 
			 * copyratio and all these estimates are averaged out
			 */
	  		groups=cleanedup.groups("Plate","CellType","CellLine","Gene");
			Table normalized=new Table();
			foreach(Table particulargene in groups)
			{
				Console.Out.Write("Working with group\n" + particulargene);
				// now we want to find all combinations of keys we can imagine
				List<Table> dilutions=particulargene.groups("Dilution");
				double average_copyratio=0;
				int count=0;
				foreach(Table dilution1 in dilutions)
					foreach(Table dilution2 in dilutions)
				{
					Double d1=Convert.ToDouble(dilution1.key["Dilution"]);
					Double d2=Convert.ToDouble(dilution2.key["Dilution"]);
					if (d1==d2) continue;
					Double cp1=(Double)dilution1.single["AverageCP"];
					Double cp2=(Double)dilution2.single["AverageCP"];
					double copyratio=estimate_efficiency(cp1,cp2,d1,d2);
					average_copyratio+=copyratio;
					count++;
				}
				average_copyratio/=count;
				Console.Out.WriteLine(" average efficiency(%) " + 100*(average_copyratio-1.0));
				// We reject copyrates larger than 2.5 and smaller than 1.8
				if (average_copyratio<1.80 || average_copyratio> 2.5) continue;
				
				// at this point we have the average copyration for this gene/celltype/plate/cellline
				// we would like to modify all the AverageCP values to take into account
				// our estimated efficiency and produce a normalized CT value.
				// another possibility could be to calculate the initial concentration directly
				// and continue working with that. We do not do that since it would make it 
				// impossible to present the standard ddct method for analyzing quantitative
				// PCR results.
				foreach(Table dilution in dilutions)
				{
					Record normalizedcp=dilution.single.copy();
					Double newcp=normalize_cp((double)dilution.single["AverageCP"],average_copyratio);
					if (newcp<40)
					{
						normalizedcp["NormalizedCP"]=newcp;
						normalized.Add(normalizedcp);
					}
				}
			}
			Console.Write(normalized);
			
			/**
			 * Lesson 12: Calculating the DCT values from any gene to GADPH.
			 * We first iterate over each element, except for the genes because
			 * we want to have all genes available for each type of experiment,
			 * so that we can find the gadph gene back.
			 */
			groups=normalized.groups("Plate","CellType","CellLine","Dilution");
			Table dct=new Table();
			foreach(Table group in groups)
			{
	  		 Record householdgeneid=new Record();
			 householdgeneid["Gene"]="GADPH";
 	 		 Table householdgene;
			 try
				{
					householdgene=group.group(householdgeneid);
				}
			 catch(KeyNotFoundException)
				{
					// the householdgene was not normalized due to missing data, meaning that we cannot estimate anything else
					// at this point. We can safely skip this dilution and continue with the next
					continue;
				}
			 // go over all the genes we have
			 List<Table> genes=group.groups("Gene");
			 foreach(Table gene in genes)
				{
					// skip if this is the householdgene
					if (gene.key["Gene"].Equals(householdgeneid["Gene"])) continue;
					// otherwise subtract the CT value
					Record dctentry=gene.single.copy();
					dctentry["Dct"]=(Double)householdgene.single["NormalizedCP"]-
						(Double)gene.single["NormalizedCP"];
					dctentry.Remove("AverageCP");
					dctentry.Remove("NormalizedCP");
					dct.Add(dctentry);
				}
			}
			Console.WriteLine("DCT Table\n--------------------");
			Console.WriteLine(dct);
			/**
			 * Lesson 13: Calculating the DDCT values from any celltype to the WT 
			 */
			Table Ddct=new Table();
			groups=dct.groups("Plate","CellLine","Gene","Dilution");
			foreach(Table gene in groups)
			{
				// find the WT cell, in this case it is called Normal
				Record normalkey=new Record();
				normalkey["CellType"]="Normal";
				Record normal;
				try{
					normal=gene.group(normalkey).single;
				}
				catch(KeyNotFoundException)
				{
					continue;
				}
				List<Table> celltypes=gene.groups("CellType");
				foreach(Table celltype in celltypes)
				{
					if (celltype.single["CellType"].Equals(normal["CellType"])) continue;
					Double dct1=(Double)normal["Dct"];
					Double dct2=(Double)celltype.single["Dct"];
					Record ddctentry=celltype.single.copy();
					ddctentry.Remove("Dct");
					ddctentry["Ddct"]=dct1-dct2;
					Ddct.Add(ddctentry);
				}			
			}
			Console.WriteLine("DDCT\n------------"+Ddct);
			/**
			 * Lesson 14: Generating a nice report, averaged over the various
			 * dilutions
			 */
			groups=Ddct.groups("Plate","CellType","CellLine","Gene");
			Table report=new Table();
			foreach(Table gene in groups)
			{
				double ratio=0;
				ArrayList ratios=gene["Ddct"];
				foreach(Double cur_ratio in ratios)
					ratio+=Math.Pow(2,cur_ratio);
				ratio/=ratios.Count;
				
				Record reportentry=gene.key.copy();
				if (ratio>1)
				{
					reportentry["Direction"]="Up";
					reportentry["Ratio"]=ratio;
				}
				else
				{
					reportentry["Direction"]="Down";
					reportentry["Ratio"]=1/ratio;
				}
				report.Add(reportentry);
			}
			Console.WriteLine("REPORT\n------------\n"+report);
			Console.ReadKey();
		}
	}
}