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