1. Preparation
A. Download the R software and Matlab7.0 software from the Internet.（http://www.r-project.org/）
B. Build the formatted excel file with the suffix “csv”, which will be connected to the software in the next step. The formatted excel file should contain three lists as shown in FigS1.A.
FigS1.A
2. R software running (Take ΔCt [β-actin] for example)
Part 1: R code for equation building
setwd("C:/users/Desktop")
——Illustration: We set the software running environment to be the computer path input.
> data1 = read.csv(file="AI.csv",header=T)
> anit=data1[data1$index=="anit"]
> anideltact=data1[data1$index=="anideltact"]
> anipmi=data1[data1$index=="anipmi"]
——Illustration: “data1” connects the file “AI.csv” to the running procedure. The parameters that will be used in the running procedure are meanwhile connected to the corresponding items in the excel file “AI.csv”. The “AI.csv” file contains the animal information; “anit” is short for “animal temperature”, “anideltact” is short for “animal △Ct value” and “anipmi” short for “animal pmi”.
> fit1_A=lm(anideltact~anit+I(anit*anipmi)+anipmi+I(anipmi^2),data=data1);
> summary(fit1_A);
——Illustration: Fitting the animal data in “AI.csv” to the hypothetical equation “anideltact=anit+I(anit*anipmi)+anipmi+I(anipmi^2)” and the R software will summarize the result of the fitting to the surfaced equation. The “I” means the coefficient of each item.
Call:
lm(formula = anideltact ~ anit + I(anit * anipmi) + anipmi + I(anipmi^2), data = data1)
Residuals:
Min 1Q Median 3Q Max
-0.95520 -0.28573 -0.01434 0.28790 1.08491
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.867e+00 1.830e-01 15.668 < 2e-16 ***
anit 4.511e-02 7.600e-03 5.936 4.52e-07 ***
I(anit * anipmi) 5.654e-04 1.134e-04 4.986 1.06e-05 ***
anipmi 3.803e-02 5.240e-03 7.257 5.48e-09 ***
I(anipmi^2) -1.579e-04 3.469e-05 -4.552 4.33e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4338 on 43 degrees of freedom
Multiple R-squared: 0.9417, Adjusted R-squared: 0.9362
F-statistic: 173.5 on 4 and 43 DF, p-value: < 2.2e-16
——Illustration: Whether the items are of statistical significance is shown by the software. Anything with a P value less than 0.01 will get a “*” and the equation will be combined. So, for example, the final equation here will be: “anideltact=0.04511*anit+0.0005654*anit*anipmi+0.03803 *anipmi-0.0001579*anipmi^2+2.867”. The coefficient and the numbers of items to combine in the equation will change if other markers are used.
> par(mfrow=c(2,2))
> plot(fit1_A)
Hit <Return> to see next plot:
——Illustration: This step shows whether the residuals between the actual Ct and fitted Ct exhibit a standard normal distribution.
Part 2: R code for estimation of the difference between real PMI and estimated PMI in the 36 verification samples
> data2 = read.csv(file="36I.csv",header=T);
> 36deltact=data2[data2$index=="36deltact"]
> realpmi=data2[data2$index=="36pmi"]
——Illustration: “data2” connects the file “36I.csv” to the running procedure. The parameters that will be used in the running procedure are meanwhile connected to the corresponding items in the excel file “36I.csv”. The “36I.csv” file contains information about the 36 rat samples. The “realpmi” is the known PMI of the 36 verification samples, while “36deltact” is short for their △Ct value.
> PMI = seq(from=0,to=120,by=0.1)
> PMI
[1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
[12] 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1
……
[1200] 119.9 120.0
> str(PMI)
num [1:1201] 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 ...
——Illustration: This step sets the “PMI” as a series of number from 0 to 120 using an interval of 0.1. In the next step, the software searches these “PMI” until the smallest deviation between its corresponding △Ct value and the known △Ct value is found, outputting the optimal PMI as our estimated PMI.
> estPMI=rep(NA,nrow(data2));
——Illustration: This step sets estPMI as a series of numbers that correspond to each verification samples in data 2.
> for (i in 1:nrow(data2)) {36deltact = data2$36deltact[i];
+ 36t = data2$36t[i];
+ id = which.min(abs(0.04511*36t+0.0005654*36t*PMI+0.03803*PMI-0.0001579*PMI^2+2.867-36deltact));
+ estPMI[i] = PMI[id];}
> estPMI
[1] 78.1 92.9 63.3 85.5 67.0 55.9 62.4 21.8 92.9 74.4 31.9
……
[144] 38.9 54.0 8.6 120.0
——Illustration: This step is a cycle operation that searches “PMI” until it finds the smallest deviation between its corresponding △Ct value and the known △Ct value in data2, and then outputs the optimal PMI as our estimated PMI.
> realPMI=data2$realPMI
> realPMI
[1] 5.50 36.00 25.00 11.00 13.50 5.50 10.00 20.50 14.50 7.00 4.00
……
[144] 20.33 13.00 6.50 10.00
> estPMI-realPMI
[1] 72.60 56.90 38.30 74.50 53.50 50.40 52.40 1.30 78.40 67.40
……
[141] 10.00 112.00 112.33 18.57 41.00 2.10 110.00
> write.csv(estPMI-realPMI, file="pmiESTIMATION.csv")
——Illustration: The final step is to output the deviation between the estimated PMI and the real PMI of each case as an excel file.
3. Matlab code for transforming the equation into visual 3 dimensional mathematical model (Take △Ct [β-actin] for example)
x=0:1:35;
y=0:4:144;
[x,y]=meshgrid(x,y);
z=0.04511*x+0.0005654*x.*y +0.03803*y -0.0001579*y.^2+2.867
surf(x,y,z);
set(gca,'XTick',[5:10:35],'FontWeight','bold');
set(gca,'YTick',[0:24:144],'FontWeight','bold');
xlabel('\bftemperature'),ylabel('\bfPMI'),zlabel('\bf\DeltaCt value');
title('\bfMyocardium(\itActb)');pause;
——Illustration: First, this step should transform the equation from “anideltact=0.04511*anit+0.0005654*anit*anipmi+0.03803*anipmi −0.0001579*anipmi^2+2.867” to “z=0.04511*x+0.0005654*x.*y +0.03803*y−0.0001579*y.^2+2.867”. The horizontal and vertical coordinates are set as the starting point and ending point with the interval value between them. There is no more tautology here regarding the general setting of the font and title. The software will output the model as FigS1.B.
FigS1.B