### R代码从vignette源的CNVtools-vignette。Rnw ' ################################################### ### 代码块1号:负载。数据 ################################################### # 源(“. . / CNVtools.r”);dyn.load(“. . / . . / src / CNVtools.so”);加载("../../CNVtools/data/A112. rdata ")库(CNVtools)数据(A112)头(A112) raw。信号<- as。矩阵(A112[， -c(1,2)]) dimnames(raw.signal)[[1]] <- A112$subject mean. ()信号<- apply(raw。信号，MAR=1, FUN=mean) pca。信号<- apply.pca(raw.signal) dir.create("fig") pdf("fig/mean_pca_signal.pdf"， width=10, height=5) par(mfrow=c(1,2)) hist(均值。信号，breaks=50, main='Mean signal'， cex.lab=1.3)信号,减免= 50,主要=第一主成分分析信号,cex.lab dev.off = 1.3) () ################################################### ### 代码块2号:模型。选择 ################################################### 批次< -因子(A112队列美元)示例< -因子(A112主题美元)set.seed(0)的结果< = pca - CNVtest.select.model(信号。信号,批=批次样品=样本,n.H0 = 3,方法=“BIC v.ncomp = 1:5, v.model.component =代表(“高斯”,5),v.model.mean =代表(“~地层(cn)”,5),v.model.var =代表(“~ 1”,5)ncomp < -结果选择美元pdf(“图/ modelselect.pdf”、宽度= 5,身高= 5)情节(结果BIC美元,xlab = n comp, ylab =“BIC类型=“b”,lty = 2, =“红色”,上校pch = ' + ') dev.off () ################################################### ### 代码块3号:集群。主成分分析 ################################################### ncomp < - 3批次< -因子(A112队列美元)示例< -因子(A112主题美元)。pca <- CNVtest。binary ( signal = pca.signal, sample = sample, batch = batches, ncomp = ncomp, n.H0=3, n.H1=0, model.var= '~ strata(cn)') print(fit.pca$status.H0) pdf("fig/pca-fit.pdf", width=10, height=5) par(mfrow=c(1,2)) cnv.plot(fit.pca$posterior.H0, batch = '58C', main = 'Cohort 58C', breaks = 50, col = 'red') cnv.plot(fit.pca$posterior.H0, batch = 'NBS', main = 'Cohort NBS', breaks = 50, col = 'red') dev.off() ################################################### ### code chunk number 4: genotype.assignment ################################################### head(fit.pca$posterior.H0) ################################################### ### code chunk number 5: ldf.improve ################################################### ncomp <- 3 pca.posterior <- as.matrix((fit.pca$posterior.H0)[, paste('P',seq(1:ncomp),sep='')]) dimnames(pca.posterior)[[1]] <- (fit.pca$posterior.H0)$subject ldf.signal <- apply.ldf(raw.signal, pca.posterior) pdf("fig/ldf_pca_signal.pdf", width=10, height=5) par(mfrow=c(1,2)) hist(pca.signal, breaks=50, main='First PCA signal', cex.lab=1.3) hist(ldf.signal, breaks=50, main='LDF signal', cex.lab=1.3) dev.off() ################################################### ### code chunk number 6: test.association ################################################### ncomp <- 3 trait <- ifelse( A112$cohort == '58C', 0, 1) fit.ldf <- CNVtest.binary ( signal = ldf.signal, sample = sample, batch = batches, disease.status = trait, ncomp = ncomp, n.H0=3, n.H1=1, model.var = "~cn") print(fit.ldf$status.H0) print(fit.ldf$status.H1) pdf("fig/ldf-fit.pdf", width=10, height=5) par(mfrow=c(1,2)) cnv.plot(fit.ldf$posterior.H0, batch = '58C', main = 'Cohort 58C', breaks = 50, col = 'red') cnv.plot(fit.ldf$posterior.H0, batch = 'NBS', main = 'Cohort NBS', breaks = 50, col = 'red') dev.off() LR.statistic <- -2*(fit.ldf$model.H0$lnL - fit.ldf$model.H1$lnL) print(LR.statistic) ################################################### ### code chunk number 7: test.association.allelic ################################################### fit.ldf <- CNVtest.binary ( signal = ldf.signal, sample = sample, batch = batches, disease.status = trait, ncomp = 3, n.H0=3, n.H1=1, model.disease = " ~ as.factor(cn)") print(fit.ldf$status.H0) print(fit.ldf$status.H1) LR.statistic <- -2*(fit.ldf$model.H0$lnL - fit.ldf$model.H1$lnL) print(LR.statistic) ################################################### ### code chunk number 8: test.association.qt ################################################### batches <- rep("ALL",length(sample)) qt <- rnorm(length(sample), mean=9.0, sd=1.0) fit.ldf <- CNVtest.qt(signal = ldf.signal, sample = sample, batch = batches, qt = qt, ncomp = ncomp, n.H0=3, n.H1=1, model.var = "~strata(cn)") print(fit.ldf$status.H0) print(fit.ldf$status.H1) LR.statistic <- -2*(fit.ldf$model.H0$lnL - fit.ldf$model.H1$lnL) print(LR.statistic) pdf("fig/qt-fit.pdf", width=15, height=5) qt.plot(fit.ldf) dev.off()