# tissue recognition analysis

# protein data from 12-09-17

library(xlsx)

library(randomForest)

library(bmrm)

library(caret)

library(dplyr)

library(data.table)

library(Boruta)

#=============  

# runs RF for tissue recognition with 10-fold CV.

# use predefined iterations and training-test assignments

# Tissue is in sampleType column of the dat object

# keep cancers with more than 150 as is, for those that are less, repeat them several times and sample the rest

runRF_balancedSample = function(dat, fold, BalanceSampleCount = 150)

{

    require(bmrm)

    require(randomForest)

    require(caret)

    require(dplyr)

#   NIter = unique(iterFold$iteration)

#print(head(dat))

    AllpredictionsNIter_tissue = c()

      MyIndices <- fold

      nTissue = nlevels(as.factor(dat$sampleType))

      Allpredictions = vector(length = length(MyIndices)) 

      Allpredictions_votes = matrix(nrow = length(MyIndices),ncol = nTissue)

      for( i in unique(MyIndices)){

        trainSamples <- which(MyIndices!=i)

        testSamples <- which(MyIndices==i)

#print(table(dat$sampleType[trainSamples]))

    trainSamplesResampled <- tapply(trainSamples,

                    dat$sampleType[trainSamples],

                    FUN = function(x)

                      if (length(x)>BalanceSampleCount)x else {

                        s = rep(x, floor(BalanceSampleCount/length(x)));

                        return(c(s,sample(x,size = BalanceSampleCount-length(s))))})  %>%

                  unlist()

    ThisRF <- randomForest(as.factor(sampleType)~.,

                           data =dat[trainSamplesResampled,],

                           cutoff = rep(1/nTissue,nTissue))

    Allpredictions[testSamples] <- predict(ThisRF,dat[testSamples,]) %>% as.character()

    Allpredictions_votes[testSamples,] <- predict(ThisRF,dat[testSamples,], type = 'vote')

      }

      names(Allpredictions) = rownames(Allpredictions_votes)=rownames(dat)

      colnames(Allpredictions_votes) = levels(as.factor(dat$sampleType))

      conf = confusionMatrix(Allpredictions,reference = dat$sampleType)

      AllpredictionsNIter_tissue <- list(prediction = Allpredictions, 

        confusion = conf, votes = Allpredictions_votes)

    print(conf$overall[1])

    return(AllpredictionsNIter_tissue)

}

# extracts positive samples and folds for an iteration and runs RF for this iteration

runOneIter = function(obj, colNames = colnames(obj), addInfo = NULL, tumSamp, sampType)

{

# positive samples only

    bahmanSamp = rownames(obj)[which(obj$LogitCall)] # 

    s = intersect(tumSamp,rownames(obj)[which(obj$LogitCall)]) # 

    print(length(s))

    if (!is.null(addInfo)) prot = data.frame(sampleType = sampType[s], obj[s,colNames],addInfo[s,]) else 

        prot = cbind(sampleType = obj[s,'Tissue'],obj[s,colNames])

    return(runRF_balancedSample(prot, fold = obj[s,'fold']))

}

#===================================================

# loading results from cancer detection

load('Results/ForLuda8proteins.rda')

#==============================

# read in protein data

protData = read.xlsx2(file = 'MEGA Protein Dataset for Random Forest, December 9, 2017.xlsx', sheetIndex = 1, row.names = 1, check.names = F, stringsAsFactors = F)

# remove Vitronectin, LRG-1, Stage, Sample Type

protOnly = protData %>% select(-`Vitronectin`, -`LRG-1`, -`Stage`, -`Sample Type`)

protOnly = sapply(protOnly, as.numeric)

#protOnly = t(do.call('rbind',protOnly))

rownames(protOnly) = rownames(protData)

colnames(protOnly) = gsub('-','',colnames(protOnly))

colnames(protOnly) = gsub(' ','',colnames(protOnly))

protOnly[is.na(protOnly)] = 0

#==================================

# create tumor type object

sampleType = as.character(protData$`Sample Type`)

# combine esophagus and stomach cancers to one group

sampleType[which(sampleType%in%c('Esophagus','Stomach'))] = 'st_eso'

names(sampleType) = rownames(protData)

# IDs of tumor and normal samples

tumSamp = names(sampleType)[which(sampleType != 'Normal')] # 1005

normSamp = names(sampleType)[which(sampleType == 'Normal')] # 812

#====================================

# patient info

sampInfo <- read.xlsx2("aar3247_Suppl. Excel_seq3_v1.xlsx",sheetIndex = 4,

                                           stringsAsFactors = F, startRow = 3, endRow = 1820)

# gender

gender = sampInfo[c("Plasma.sample.ID..","Sex")]

gender = gender[!duplicated(gender[,1]),]

n = gender[,1]

gender = gender[,2]

names(gender) = n

gender[which(gender %in% c('Female','female'))] = 'F'

gender[which(gender %in% c('Male','male'))] = 'M'

gender[which(gender %in% c('unknown',''))] = NA

#=============================

# add matrix with omegas per gene per patient.

#=============================

# make 10 matrices from Analysis_1

omegaPerPatient = omegaPerPatientK13 = vector(mode = 'list',length = 10)

s = c(tumSamp, normSamp)

for(i in 1:10)

{

    load(paste0('allResults_fromLuMethod_',i,'_20171209.rda'))

    perPat = split(allResults,allResults$Template)

    perPat = lapply(perPat, function(x) x[which(!x[,'blacklist']),c('Gene','omega')])

    # get maximum omega for each gene

    perPatMaxOmega = sapply(perPat, function(x) tapply(x$omega, INDEX = x$Gene, FUN = max, na.rm = T))

    genes = unique(allResults$Gene)

    s1 = union(s,names(perPatMaxOmega))

    omegaPerPatient[[i]] = matrix(0,ncol = length(genes), nrow = length(s1), dimnames = list(s1,genes))

    for(j in names(perPatMaxOmega)) omegaPerPatient[[i]][j,names(perPatMaxOmega[[j]])] = perPatMaxOmega[[j]]

    # separate KRAS codon 13

    resultKRAS13 = allResults %>% filter(Gene == 'KRAS', codonStart == 13)

    resultNoK13 = allResults %>% filter(Gene != 'KRAS' | codonStart != 13)

    resultKRAS13[,'Gene'] = 'KRAS13'

    allResults = rbind(resultNoK13,resultKRAS13)

    # get maximum omega for each gene

    perPatMaxOmega = sapply(perPat, function(x) tapply(x$omega, INDEX = x$Gene, FUN = max, na.rm = T))

    genes = unique(allResults$Gene)

    s1 = union(s,names(perPatMaxOmega))

    omegaPerPatientK13[[i]] = matrix(0,ncol = length(genes), nrow = length(s1), dimnames = list(s1,genes))

    for(j in names(perPatMaxOmega)) omegaPerPatientK13[[i]][j,names(perPatMaxOmega[[j]])] = perPatMaxOmega[[j]]

}

#===========================

# run all proteins, gender and mutations

set.seed(12485)

res8 = vector(length = 10)

for(i in 1:10)

{

    res8[i] = list(runOneIter(cbind(protOnly,ForLuda[[i]][rownames(protOnly),c('LogitCall','fold')]),

            colNames = 1:39, 

            addInfo = data.frame(gender = gender[rownames(protOnly)],omegaPerPatient[[i]][rownames(protOnly),]), 

            tumSamp = tumSamp, sampType = sampleType))

}

# save results for the paper

tab = sapply(res8, FUN = function(x) c(x$confusion$overall[1], x$confusion$byClass[,1])) 

tab = cbind(tab,mean = apply(tab,1,mean))

# mean accuracy 0.6968313

xlsFile = paste('Results/tissueRecog_results_',Sys.Date(),'.xlsx')

write.xlsx(tab, file = xlsFile, sheetName = 'prot +  mut')

# confusion matrix with fractions

conf = res8[[4]]$confusion$table

write.csv(sweep(conf, 2, apply(conf,2,sum),'/'), file = paste('Results/tissueRecog_results_',Sys.Date(),'.csv'), quote = F)

# per samples votes

tab = cbind(res8[[4]]$votes, prediction = res8[[4]]$prediction, tissue = sampleType[names(res8[[4]]$prediction)])

write.xlsx(tab, file = xlsFile, sheetName = 'Iteration 4. Votes', append = T)