# Simulation functions for telomere measurement error project # Daniel Nettle, April 2018 ### TS ratio calculation function ts()### # Calculates TS ratio given cqs ts=function(tel, scg, ref.tel, ref.scg) { delta.cq.sample=tel-scg delta.cq.reference=ref.tel - ref.scg delta.delta.cq=delta.cq.sample-delta.cq.reference ts=2^(-delta.delta.cq) return(ts) } ###Data generation function generate.one.dataset() ##### # Generates a data under given parameters # Allows comparison of true and measured values generate.one.dataset=function( n=10000, b=10, a=1000, f=28, var.sample.size=1, telomere.var=0.1, error.telo=0.05, error.scg=0.05, error.cor=0){ # Calculate the true DNA and Cq values true.dna.scg=rnorm(n, b, sd=var.sample.size) true.dna.scg[true.dna.scg<=0]=0.0000001 true.cq.scg=f-log2(true.dna.scg) true.telo.var=rnorm(n, 1, sd=telomere.var) true.dna.telo=a*true.dna.scg*true.telo.var true.telo.cq=f - log2(true.dna.telo) # Now add measurement errors scaled.error.for.scg=rnorm(n) scaled.error.for.telo=error.cor*scaled.error.for.scg + sqrt(1-error.cor^2)*rnorm(n) error.for.scg=error.scg*scaled.error.for.scg error.for.telo=error.telo*scaled.error.for.telo measured.cq.scg=true.cq.scg + error.for.scg measured.cq.telo=true.telo.cq + error.for.telo # Now calculate T/S ratios true.ts=ts(tel=true.telo.cq, scg=true.cq.scg, ref.tel=mean(true.telo.cq), ref.scg = mean(true.cq.scg)) measured.ts=ts(tel=measured.cq.telo, scg=measured.cq.scg, ref.tel=mean(true.telo.cq), ref.scg=mean(measured.cq.scg)) # And the errors error.computed=measured.ts - true.ts error.analytic = (2^(error.for.scg - error.for.telo) - 1)*true.telo.var data=data.frame(true.telo.var, true.dna.scg, true.telo.cq, true.cq.scg, true.ts, measured.cq.telo, measured.cq.scg, error.for.telo, error.for.scg, measured.ts, error.analytic, error.computed) return(data) # End function } ####Generated repeated measures function generate.repeated.measure()#### # Generates a dataset where the same individuals are measured twice with no true change generate.repeated.measure=function( n=10000, b=10, a=1000, f=28, var.sample.size=1, telomere.var=0.1, error.telo=0.05, error.scg=0.05, error.cor=0){ # Calculate the true DNA and cq values at each time point (true rtl is the same) true.dna.scg=rnorm(n, b, sd=var.sample.size) true.dna.scg[true.dna.scg<=0]=0.0000001 true.dna.scg.2=rnorm(n, b, sd=var.sample.size) true.dna.scg.2[true.dna.scg.2<=0]=0.0000001 true.cq.scg=f-log2(true.dna.scg) true.cq.scg.2=f-log2(true.dna.scg.2) true.telo.var=rnorm(n, 1, sd=telomere.var) true.dna.telo=a*true.dna.scg*true.telo.var true.telo.cq=f - log2(true.dna.telo) true.dna.telo.2=a*true.dna.scg.2*true.telo.var true.telo.cq.2=f - log2(true.dna.telo.2) # Now add measurement error for two assays scaled.error.for.scg=rnorm(n) scaled.error.for.telo=error.cor*scaled.error.for.scg + sqrt(1-error.cor^2)*rnorm(n) error.for.scg=error.scg*scaled.error.for.scg error.for.telo=error.telo*scaled.error.for.telo measured.cq.scg=true.cq.scg + error.for.scg measured.cq.telo=true.telo.cq + error.for.telo scaled.error.for.scg.2=rnorm(n) scaled.error.for.telo.2=error.cor*scaled.error.for.scg.2 + sqrt(1-error.cor^2)*rnorm(n) error.for.scg.2=error.scg*scaled.error.for.scg.2 error.for.telo.2=error.telo*scaled.error.for.telo.2 measured.cq.scg.2=true.cq.scg.2 + error.for.scg.2 measured.cq.telo.2=true.telo.cq.2 + error.for.telo.2 # Now calculate T/S ratios true.ts.1=ts(tel=true.telo.cq, scg=true.cq.scg, ref.tel = mean(true.telo.cq), ref.scg = mean(true.cq.scg)) measured.ts.1=ts(tel=measured.cq.telo, scg=measured.cq.scg, ref.tel = mean(measured.cq.telo), ref.scg = mean(measured.cq.scg)) true.ts.2=ts(tel=true.telo.cq.2, scg=true.cq.scg.2, ref.tel = mean(true.telo.cq.2), ref.scg = mean(true.cq.scg.2)) measured.ts.2=ts(tel=measured.cq.telo.2, scg=measured.cq.scg.2, ref.tel = mean(measured.cq.telo.2), ref.scg = mean(measured.cq.scg.2)) # And the errors error.computed=measured.ts.1 - true.ts.1 error.analytic = (2^(error.for.scg - error.for.telo) - 1)*a*true.telo.var error.computed.2=measured.ts.2 - true.ts.2 error.analytic.2 = (2^(error.for.scg.2 - error.for.telo.2) - 1)*a*true.telo.var data=data.frame(true.telo.var, true.dna.scg, true.dna.scg.2, true.dna.telo, true.telo.cq, true.cq.scg, error.for.telo, error.for.scg, error.for.telo.2, error.for.scg.2, measured.cq.telo, measured.cq.scg, measured.cq.telo.2, measured.cq.scg.2, true.ts.1, true.ts.2, measured.ts.1, measured.ts.2, error.computed, error.analytic, error.computed.2, error.analytic.2) return(data) # End function } ###calculate.repeatability() function#### # Calculates the ICC of the TS ratio when the same individuals are run twice calculate.repeatability=function(n = 10000, error.telo, error.scg, error.cor=0, a=1000, var.sample.size=1){ require(irr) d=generate.repeated.measure(n = n, error.telo=error.telo, error.scg=error.scg, error.cor=error.cor, a=a, var.sample.size=var.sample.size) m=icc(data.frame(d$measured.ts.1, d$measured.ts.2)) return(m$value) } ###compare.repeatability() function#### # For comparing repeatability of T/S and Cq telo values compare.repeatibility=function(n=10000, error.telo, error.scg, error.cor=0, a=1000, var.sample.size=1) { require(irr) d=generate.repeated.measure(n=n, error.telo=error.telo, error.scg=error.scg, error.cor=error.cor, a=a, var.sample.size = var.sample.size) m1=icc(data.frame(d$measured.cq.telo, d$measured.cq.telo.2)) m2=icc(data.frame(d$measured.cq.scg, d$measured.cq.scg.2)) m3=icc(data.frame(d$measured.ts.1, d$measured.ts.2)) return(c(m1$value, m2$value, m3$value)) # End function }