tcgaagedx.R

library(SEMMS)
library(QREM)
library(xtable)
QRVS <- function(fn,responseCol,predictors, nnset, qn) {
  dataYXZ <- readInputFile(fn, ycol=responseCol, Zcols=predictors)
  nn <- 10
  n  <- dataYXZ$N
  y0 <- dataYXZ$Y
  K  <- dataYXZ$K
  for (rep in 1:maxReps) {
    # create a subset of the selected columns and run QREM
    #cat(responseCol, K, nnset,"\n")
    preds <- paste(colnames(dataYXZ$Z)[nnset], collapse = "+")
    linmod <- as.formula(paste("Y ~", preds))
    dframetmp <- data.frame(cbind(dataYXZ$Y, dataYXZ$Z[,nnset]))
    colnames(dframetmp) <- c("Y",colnames(dataYXZ$Z)[nnset])
    qremFit <- QREM(lm,linmod, dframetmp, qn, maxInvLambda = 1000)
    # apply the weights found by QREM and rerun SEMMS
    dataYXZtmp <- dataYXZ
    dataYXZtmp$Y <- (dataYXZ$Y-(1-2*qn)/qremFit$weights)
    fittedVSnew <- fitSEMMS(dataYXZtmp,distribution = 'N', mincor=0.8,rnd=F,
                            nnset=nnset, minchange = 1, maxst = 20)
    foundSEMMSnew <- sort(union(which(fittedVSnew$gam.out$lockedOut != 0),
                                fittedVSnew$gam.out$nn))
    if (length(fittedVSnew$gam.out$nn) == length(nnset)) {
      if (all(fittedVSnew$gam.out$nn == nnset)) {
        break
      }
    }
    nnset <- fittedVSnew$gam.out$nn
    if (length(nnset) == 0)
      return(list(selected="none",coef="none"))
  }
  print(summary(qremFit$fitted.mod))
  return(list(selected=dataYXZ$originalZnames[nnset], coef=qremFit$coef$beta))
}


load("TCGAlung.RData")
set.seed(772233)
smokers <- sort(sample(which(clin5$smoking2 =="S"), 101))
nonsmokers <- which(clin5$smoking2 =="N")
sset <- sort(c(smokers,nonsmokers))
dat <- data.frame(clin5$ageDX[sset],clin5$smoking[sset],t(GE5[,sset]))
save(dat,file="agedx.RData")
smk <- which(clin5$smoking2 == "S")
dat <- data.frame(clin5$ageDX[smk],t(GE5[,smk]))
save(dat, file="TCGAsmokers.RData")

fn = "TCGAsmokers.RData"
dataYXZ <- readInputFile(fn, ycol = 1, Zcols = 2:13493)
responseCol = 1
nn =10
maxReps <- 20

# crall <- cor(cbind(dataYXZ$Y, dataYXZ$Z))
# nnset <- setdiff(order(abs(crall[responseCol,]), decreasing = TRUE)[1:nn],responseCol)
nnset <- c(4705,9189,7114,2416,1019,3506,10155,8726,7261)

qn=0.85 
qrvs <- QRVS(fn, responseCol, 2:13493, nnset, qn)
# https://www.creativebiomart.net/symbolsearch_PTGES3.htm?msclkid=4df7d3a778641495f59f8d986030fdc0&utm_source=bing&utm_medium=cpc&utm_campaign=Biomart%201%20Gene%20Family&utm_term=PTGES3&utm_content=PTGES3
# Z9189 PTGES3.10728
#(Intercept)       Z9189 
#  74.111830   -1.260368

qn = 0.5
qrvs <- QRVS(fn, responseCol, 2:13493, nnset, qn)
# none


qn = 0.15
qrvs <- QRVS(fn, responseCol, 2:13493, nnset, qn)
# https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3594794/
# Z10155 SCO1.6341
# (Intercept)      Z10155 
#   58.088002   -2.271105 


preds <- paste(colnames(dataYXZ$Z)[9189], collapse = "+")
linmod <- as.formula(paste("Y ~", preds))
dframetmp <- data.frame(cbind(dataYXZ$Y, dataYXZ$Z[,9189]))
colnames(dframetmp) <- c("Y",colnames(dataYXZ$Z)[9189])
qremFit85 <- QREM(lm,linmod, dframetmp, 0.85, maxInvLambda = 1000)
bcov85 <- bcov(qremFit85, linmod, dframetmp, 0.85)
cat(qremFit85$coef$beta,sqrt(diag(bcov85)),"\n")
QRdiagnostics(dframetmp$Z9189,dataYXZ$originalZnames[9189], qremFit85$ui, 0.85)

plot(dframetmp$Z9189,dframetmp$Y, col=(1+(qremFit85$ui>0)))
abline(qremFit85$coef$beta)

plot(dat[,9190],dat[,1], col=(1+(qremFit85$ui>0)))
abline(qremFit85$coef$beta[1]-qremFit85$coef$beta[2]*mean(dat[,9190])/sd(dat[,9190]), 
       qremFit85$coef$beta[2]/sd(dat[,9190]))

preds <- paste(colnames(dataYXZ$Z)[10155], collapse = "+")
linmod <- as.formula(paste("Y ~", preds))
dframetmp <- data.frame(cbind(dataYXZ$Y, dataYXZ$Z[,10155]))
colnames(dframetmp) <- c("Y",colnames(dataYXZ$Z)[10155])
qremFit15 <- QREM(lm,linmod, dframetmp, 0.15, maxInvLambda = 1000)
bcov15 <- bcov(qremFit15, linmod, dframetmp, 0.15)
cat(qremFit15$coef$beta,sqrt(diag(bcov15)),"\n")
QRdiagnostics(dframetmp$Z10155,dataYXZ$originalZnames[10155], qremFit15$ui, 0.15)

plot(dframetmp$Z10155,dframetmp$Y, col=(1+(qremFit15$ui>0)))
abline(qremFit15$coef$beta)

plot(dat[,10156],dat[,1], col=(1+(qremFit15$ui>0)))
abline(qremFit15$coef$beta[1]-qremFit15$coef$beta[2]*mean(dat[,10156])/sd(dat[,10156]), 
       qremFit15$coef$beta[2]/sd(dat[,10156]))

#####

stage <- as.factor(ifelse(clin5$ajcc_pathologic_tumor_stage[smk] %in% c("Stage III","Stage IIIA","Stage IIIB","Stage IV"), "A","E"))
gender <- droplevels(clin5$gender[smk])
linmod <- as.formula("Y ~ Z9189+stage")
dframetmp <- data.frame(cbind(dataYXZ$Y, dataYXZ$Z[,9189]),stage,gender)
colnames(dframetmp) <- c("Y","Z9189","stage","gender")
qremFit85b <- QREM(lm,linmod, dframetmp, 0.85, maxInvLambda = 1000)
bcov85b <- bcov(qremFit85b, linmod, dframetmp, 0.85)
xtable(cbind(qremFit85b$coef$beta,sqrt(diag(bcov85b))),"\n")
QRdiagnostics(dframetmp$Z9189,dataYXZ$originalZnames[9189], qremFit85$ui, 0.85)