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 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
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
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
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(); } } }