Comprehensive analysis of interactions between brain aging and late-onset psychoses using heuristic mapping models.

1. #' This function performs variable selection on gene expression matrices.
2. #' It can, for instance, remove genes with low variation.
3. #' @param exprMat A matrix of gene expression levels. rownames() are genes, and colnames() are samples.
4. #' @param removeLowVaryingGenes The proportion of low varying genes to be removed.The default is .2
5. #' @return A vector of row/genes to keep.
6. #' @export
7. doVariableSelection <- function(exprMat, removeLowVaryingGenes=.2)
8. {
9. vars <- apply(exprMat, 1, var)
10. return(order(vars, decreasing=TRUE)[seq(1:as.integer(nrow(exprMat)*(1-removeLowVaryingGenes)))])
11. }
12. #'This function takes two gene expression matrices (like trainExprMat and testExprMat) and returns homogenized versions of the matrices by employing the homogenization method specified.
13. #'By default, the Combat method from the sva library is used.
14. #'In both matrices, genes are row names and samples are column names.
15. #'It will deal with duplicated gene names, as it subsets and orders the matrices correctly.
16. #'@param testExprMat A gene expression matrix for samples on which we wish to predict a phenotype.Genes are rows, samples are columns.
17. #'@param trainExprMat A gene expression matrix for samples for which the phenotype is already known.Genes are rows, samples are columns.
18. #'@param batchCorrect The type of batch correction to be used. Options are 'eb' for Combat, 'none', or 'qn' for quantile normalization.
19. #'#The default is 'eb'.
20. #'@param selection This parameter can be used to specify how duplicates are handled. The default value of -1 means to ask the user.
21. #'#Other options include '1' to summarize duplicates by their mean, and '2'to discard all duplicated genes.
22. #'@param printOutput To suppress output, set to false. Default is TRUE.
23. #'@import sva
24. #'@import preprocessCore
25. #'@keywords Homogenize gene expression data.
26. #'@return A list containing two entries $train and $test, which are the homogenized input matrices.
27. #'@export
28. homogenizeData<-function (testExprMat, trainExprMat, batchCorrect = "eb", selection = -1, printOutput = TRUE)
29. {
30. #Check the batchCorrect parameter
31. if (!(batchCorrect %in% c("eb", "qn", "none",
32. "rank", "rank_then_eb", "standardize")))
33. stop("\"batchCorrect\" must be one of \"eb\", \"qn\", \"rank\", \"rank_then_eb\", \"standardize\" or \"none\"")
35. #Check if both row and column names have been specified.
36. if (is.null(rownames(trainExprMat)) || is.null(rownames(testExprMat))) {
37. stop("ERROR: Gene identifiers must be specified as \"rownames()\" on both training and test expression matrices. Both matices must have the same type of gene identifiers.")
38. }
40. #Check that some of the row names overlap between both datasets (print an error if none overlap)
41. if (sum(rownames(trainExprMat) %in% rownames(testExprMat)) == 0) {
42. stop("ERROR: The rownames() of the supplied expression matrices do not match. Note that these are case-sensitive.")
43. }
44. else {
45. if (printOutput)
46. message(paste("\n", sum(rownames(trainExprMat) %in% rownames(testExprMat)), " gene identifiers overlap between the supplied expression matrices... \n", paste = ""))
47. }
49. #If there are duplicate gene names, give the option of removing them or summarizing them by their mean.
50. if ((sum(duplicated(rownames(trainExprMat))) > 0) || sum(sum(duplicated(rownames(testExprMat))) > 0)) {
51. if (selection == -1) {
52. message("\nExpression matrix contain duplicated gene identifiers (i.e. duplicate rownames()), how would you like to proceed:")
53. message("\n1. Summarize duplicated gene ids by their mean value (acceptable in most cases)")
54. message("\n2. Disguard all duplicated genes (recommended if unsure)")
55. message("\n3. Abort (if you want to deal with duplicate genes ids manually)\n")
56. }
57. while (is.na(selection) | selection <= 0 | selection > 3) {
58. selection <- readline("Selection: ")
59. selection <- ifelse(grepl("[^1-3.]", selection), -1, as.numeric(selection))
60. }
62. message("\n")
64. if (selection == 1) #Summarize duplicates by their mean.
65. {
66. if ((sum(duplicated(rownames(trainExprMat))) > 0)) {
67. trainExprMat <- summarizeGenesByMean(trainExprMat)
68. }
69. if ((sum(duplicated(rownames(testExprMat))) > 0)) {
70. testExprMat <- summarizeGenesByMean(testExprMat)
71. }
72. }
73. else if (selection == 2) #Disguard all duplicated genes.
74. {
75. if ((sum(duplicated(rownames(trainExprMat))) > 0)) {
76. keepGenes <- names(which(table(rownames(trainExprMat)) == 1))
77. trainExprMat <- trainExprMat[keepGenes, ]
78. }
79. if ((sum(duplicated(rownames(testExprMat))) > 0)) {
80. keepGenes <- names(which(table(rownames(testExprMat)) == 1))
81. testExprMat <- testExprMat[keepGenes, ]
82. }
83. }
84. else {
85. stop("Exectution Aborted!")
86. }
87. }
89. #Subset and order gene ids on the expression matrices.
90. commonGenesIds <- rownames(trainExprMat)[rownames(trainExprMat) %in%
91. rownames(testExprMat)]
92. trainExprMat <- trainExprMat[commonGenesIds, ]
93. testExprMat <- testExprMat[commonGenesIds, ]
95. #Subset and order the 2 expression matrices.
96. if (batchCorrect == "eb") {
97. #Subset to common genes and batch correct using ComBat.
98. dataMat <- cbind(trainExprMat, testExprMat)
99. mod <- data.frame(`(Intercept)` = rep(1, ncol(dataMat)))
100. rownames(mod) <- colnames(dataMat)
101. whichbatch <- as.factor(c(rep("train", ncol(trainExprMat)),
102. rep("test", ncol(testExprMat))))
104. # Added
105. # Filter out genes with low variances to make sure comBat run correctly
106. dataMat <- cbind(trainExprMat, testExprMat)
107. gene_vars = apply(dataMat, 1, var)
108. genes<-as.vector(gene_vars)
110. if (length(which(genes <= 1e-3) != 0)){ #If some genes have low variances (if the variance is not 0), remove them.
111. dataMat = dataMat[-(which(genes <= 1e-3)),]
112. }
114. # End added
116. combatout <- ComBat(dataMat, whichbatch, mod = mod)
117. return(list(train = combatout[, whichbatch == "train"],
118. test = combatout[, whichbatch == "test"], selection = selection))
119. }
120. else if (batchCorrect == "standardize") #Standardize to mean 0 and variance 1 in each dataset using a non EB based approach.
121. {
122. for (i in 1:nrow(trainExprMat)) {
123. row <- trainExprMat[i, ]
124. trainExprMat[i, ] <- ((row - mean(row))/sd(row))
125. }
126. for (i in 1:nrow(testExprMat)) {
127. row <- testExprMat[i, ]
128. testExprMat[i, ] <- ((row - mean(row))/sd(row))
129. }
130. return(list(train = trainExprMat, test = testExprMat,
131. selection = selection))
132. }
133. else if (batchCorrect == "rank") #The random-rank transform approach, that may be better when applying models to RNA-seq data.
134. {
135. for (i in 1:nrow(trainExprMat)) {
136. trainExprMat[i, ] <- rank(trainExprMat[i, ], ties.method = "random")
137. }
138. for (i in 1:nrow(testExprMat)) {
139. testExprMat[i, ] <- rank(testExprMat[i, ], ties.method = "random")
140. }
141. return(list(train = trainExprMat, test = testExprMat,
142. selection = selection))
143. }
144. else if (batchCorrect == "rank_then_eb") #Rank-transform the RNAseq data, then apply ComBat
145. {
146. #First, rank transform the RNA-seq data.
147. for (i in 1:nrow(testExprMat)) {
148. testExprMat[i, ] <- rank(testExprMat[i, ], ties.method = "random")
149. }
150. #Subset to common genes and batch correct using ComBat.
151. dataMat <- cbind(trainExprMat, testExprMat)
152. mod <- data.frame(`(Intercept)` = rep(1, ncol(dataMat)))
153. rownames(mod) <- colnames(dataMat)
154. whichbatch <- as.factor(c(rep("train", ncol(trainExprMat)),
155. rep("test", ncol(testExprMat))))
156. combatout <- ComBat(dataMat, whichbatch, mod = mod)
157. return(list(train = combatout[, whichbatch == "train"],
158. test = combatout[, whichbatch == "test"], selection = selection))
159. }
160. else if (batchCorrect == "qn")
161. {
162. dataMat <- cbind(trainExprMat, testExprMat)
163. dataMatNorm <- normalize.quantiles(dataMat)
164. whichbatch <- as.factor(c(rep("train", ncol(trainExprMat)),
165. rep("test", ncol(testExprMat))))
166. return(list(train = dataMatNorm[, whichbatch == "train"],
167. test = dataMatNorm[, whichbatch == "test"],
168. selection = selection))
169. }
170. else {
171. return(list(train = trainExprMat, test = testExprMat,
172. selection = selection))
173. }
174. }
175. #'This function takes a gene expression matrix and if duplicate genes are measured, summarizes them by their means.
176. #'@param exprMat A gene expression matrix with genes as rownames() and samples as colnames().
177. #'@return A gene expression matrix that does not contain duplicate genes.
178. #'@keywords Summarize duplicate genes by their mean.
179. #'@export
180. summarizeGenesByMean <- function(exprMat)
181. {
182. geneIds <- rownames(exprMat)
183. t <- table(geneIds) #How many times is each gene name duplicated.
184. allNumDups <- unique(t)
185. allNumDups <- allNumDups[-which(allNumDups == 1)]
187. #Create a *new* gene expression matrix with everything in the correct order....
188. #Starrt by just adding stuff that isn't duplicated
189. exprMatUnique <- exprMat[which(geneIds %in% names(t[t == 1])), ]
190. gnamesUnique <- geneIds[which(geneIds %in% names(t[t == 1]))]
192. #Add all the duplicated genes to the bottom of "exprMatUniqueHuman", summarizing as you go
193. for(numDups in allNumDups)
194. {
195. geneList <- names(which(t == numDups))
197. for(i in 1:length(geneList))
198. {
199. exprMatUnique <- rbind(exprMatUnique, colMeans(exprMat[which(geneIds == geneList[i]), ]))
200. gnamesUnique <- c(gnamesUnique, geneList[i])
201. # print(i)
202. }
203. }
205. if(class(exprMatUnique) == "numeric")
206. {
207. exprMatUnique <- matrix(exprMatUnique, ncol=1)
208. }
210. rownames(exprMatUnique) <- gnamesUnique
211. return(exprMatUnique)
212. }
213. #'This function predicts a phenotype (drug sensitivity score) when provided with microarray or bulk RNAseq gene expression data of different platforms.
214. #'The imputations are performed using ridge regression, training on a gene expression matrix where phenotype is already known.
215. #'This function integrates training and testing datasets via a user-defined procedure, and power transforming the known phenotype.
216. #'@param trainingExprData The training data. A matrix of expression levels. rownames() are genes, colnames() are samples (cell line names or cosmic ides, etc.). rownames() must be specified and must contain the same type of gene ids as "testExprData"
217. #'@param trainingPtype The known phenotype for "trainingExprData". This data must be a matrix of drugs/rows x cell lines/columns or cosmic ids/columns. This matrix can contain NA values, that is ok (they are removed in the calcPhenotype() function).
218. #'@param testExprData The test data where the phenotype will be estimated. It is a matrix of expression levels, rows contain genes and columns contain samples, "rownames()" must be specified and must contain the same type of gene ids as "trainingExprData".
219. #'@param batchCorrect How should training and test data matrices be homogenized. Choices are "eb" (default) for ComBat, "qn" for quantiles normalization or "none" for no homogenization.
220. #'@param powerTransformPhenotype Should the phenotype be power transformed before we fit the regression model? Default to TRUE, set to FALSE if the phenotype is already known to be highly normal.
221. #'@param removeLowVaryingGenes What proportion of low varying genes should be removed? 20 percent be default
222. #'@param minNumSamples How many training and test samples are required. Print an error if below this threshold
223. #'@param selection How should duplicate gene ids be handled. Default is -1 which asks the user. 1 to summarize by their or 2 to disguard all duplicates.
224. #'@param printOutput Set to FALSE to supress output.
225. #'@param pcr Indicates whether or not you'd like to use pcr for feature (gene) reduction. Options are 'TRUE' and 'FALSE'. If you indicate 'report_pc=TRUE' you need to also indicate 'pcr=TRUE'
226. #'@param removeLowVaringGenesFrom Determine method to remove low varying genes. Options are 'homogenizeData' and 'rawData'.
227. #'@param report_pc Indicates whether you want to output the training principal components. Options are 'TRUE' and 'FALSE'.
228. #'@param cc Indicate if you want correlation coefficients for biomarker discovery.
229. #'@param percent Indicate percent variability (of the training data) you'd like principal components to reflect if pcr=TRUE. Default is 80 for 80%
230. #'These are the correlations between a given gene of interest across all samples vs. a given drug response across samples.
231. #'These correlations can be ranked to obtain a ranked correlation to determine highly correlated drug-gene associations.
232. #'@param rsq Indicate whether or not you want to output the R^2 values for the data you train on from true and predicted values.
233. #'These values represent the percentage in which the optimal model accounts for the variance in the training data.
234. #'Options are 'TRUE' and 'FALSE'.
235. #'@return .txt files will be saved into your working directory. Depending on the parameter specified, the .txt file outputs of this function can include the estimated phenotype/sensitivity predictions,
236. #'the R^2 data, and the correlation coefficients. Principal components are stored as .RData files for each drug in your drug dataset.
237. #'@import sva
238. #'@import ridge
239. #'@import car
240. #'@import ridge
241. #'@import glmnet
242. #'@import tidyverse
243. #'@import utils
244. #'@import stats
245. #'@importFrom pls explvar explvar
246. #'@keywords predict drug sensitivity and phenotype
247. #'@export
248. calcPhenotype<-function (trainingExprData,
249. trainingPtype,
250. testExprData,
251. batchCorrect,
252. powerTransformPhenotype=TRUE,
253. removeLowVaryingGenes=0.2,
254. minNumSamples,
255. selection=1,
256. printOutput,
257. pcr=FALSE,
258. removeLowVaringGenesFrom,
259. report_pc=FALSE,
260. cc=FALSE,
261. percent=80,
262. rsq=FALSE)
263. {
265. #Initiate empty lists for each data type you'd like to collect.
266. #_______________________________________________________________
267. DrugPredictions<-list() #Collects drug predictions.
268. rsqs<-list() #Collects R^2 values.
269. cors<-list() #Collects correlation coefficient for each gene across all samples vs. each drug across all samples.
270. pvalues<-list() #Collects p-values for the correlation coefficients.
272. #vs=c()
274. drugs<-colnames(trainingPtype) #Store all the possible drugs in a vector.
276. #Check the supplied data and parameters.
277. #_______________________________________________________________
278. if (class(testExprData)[1] != "matrix")
279. stop("\nERROR: \"testExprData\" must be a matrix.")
280. if (class(trainingExprData)[1] != "matrix")
281. stop("\nERROR: \"trainingExprData\" must be a matrix.")
282. if (class(trainingPtype)[1] != "matrix")
283. stop("\nERROR: \"trainingPtype\" must be a matrix.")
285. if (report_pc)
286. if (pcr == FALSE)
287. stop("\nERROR: pcr must be TRUE if report_pc is TRUE")
289. if (pcr)
290. if (cc)
291. stop("\nERROR: pcr must be FALSE if cc is TRUE")
293. #Make sure training samples are equivalent in both matrices.
294. if (!any(colnames(trainingExprData) %in% rownames(trainingPtype)))
295. stop("\nERROR: No Cell Lines Found in Common: Sample names must be consistent in training matrices")
297. #Subset and order the training Expr and trainingPtype to the cell lines in common (and order them)
299. commonCellLines<-colnames(trainingExprData)[colnames(trainingExprData) %in% rownames(trainingPtype)]
301. trainingExprData <- trainingExprData[,commonCellLines]
302. trainingPtype <- trainingPtype[commonCellLines,]
304. #Check if an adequate number of training and test samples have been supplied.
305. #_______________________________________________________________
306. if ((nrow(trainingExprData) < minNumSamples) || (nrow(testExprData) < minNumSamples)) {
307. stop(paste("\nThere are less than", minNumSamples, "samples in your test or training set. It is strongly recommended that you use larger numbers of samples in order to (a) correct for batch effects and (b) fit a reliable model. To supress this message, change the \"minNumSamples\" parameter to this function."))
308. }
310. #Get the homogenized data.
311. #_______________________________________________________________
312. homData <- homogenizeData(testExprMat=testExprData, trainExprMat=trainingExprData, batchCorrect, selection, printOutput)
314. #Remove low varying genes.
315. #_______________________________________________________________
316. #Do variable selection if specified. By default, we remove 20% of least varying genes from the homogenized dataset.
317. #We can also remove the intersection of the lowest 20% from both training and test sets (for the gene ids remaining in the homogenized data).
318. #Otherwise, keep all genes.
320. #Check batchCorrect parameter.
321. if (!(removeLowVaringGenesFrom %in% c("homogenizeData", "rawData"))) {
322. stop("\nremoveLowVaringGenesFrom\" must be one of \"homogenizeData\", \"rawData\"")
323. }
325. keepRows <- seq(1:nrow(homData$train)) #By default we will keep all the genes.
326. if (removeLowVaryingGenes > 0 && removeLowVaryingGenes < 1) { #If the proportion of variability to keep is between 0 and 1.
327. if (removeLowVaringGenesFrom == "homogenizeData") { #If you're filtering based on homogenized data.
328. keepRows <- doVariableSelection(cbind(homData$test, homData$train), removeLowVaryingGenes = removeLowVaryingGenes)
330. numberGenesRemoved <- nrow(homData$test) - length(keepRows)
331. if (printOutput) message(paste("\n", numberGenesRemoved, "low variabilty genes filtered."));
332. }
333. else if (removeLowVaringGenesFrom == "rawData") { #If we are filtering based on the raw data i.e. the intersection of the things filtered from both datasets.
334. evaluabeGenes <- rownames(homData$test)
335. keepRowsTrain <- doVariableSelection(trainingExprData[evaluabeGenes,], removeLowVaryingGenes = removeLowVaryingGenes)
336. keepRowsTest <- doVariableSelection(testExprData[evaluabeGenes,], removeLowVaryingGenes = removeLowVaryingGenes)
337. keepRows <- intersect(keepRowsTrain, keepRowsTest)
338. numberGenesRemoved <- nrow(homData$test) - length(keepRows)
339. if (printOutput)
340. message(paste("\n", numberGenesRemoved, "low variabilty genes filtered."));
341. }
342. }
344. #Predict for each drug.
345. #_______________________________________________________________
346. for(a in 1:length(drugs)){ #For each drug...
348. #Modify trainingPtype and trainingExprData so that you only use cell lines for which you have expression and response data for.
349. #_______________________________________________________________
350. trainingPtype2<-trainingPtype[,a] #Obtain the response data for the compound of interest.
351. NonNAindex <- which(!is.na(trainingPtype2)) #Get the indices of the non NAs. You only want the cell lines/cosmic ids that you have drug info for.
353. samps<-rownames(trainingPtype)[NonNAindex] #Obtain cell lines you have expression and response data for.
355. if (length(samps) == 1){ #Make sure training data has more than just 1 sample. If the drug has one sample, it will be skipped.
356. drugs = drugs[-a]
357. message(paste("\n", drugs[a], "is skipped due to insufficient cell lines to fit the model."))
358. next
359. } else {
361. trainingPtype4<-as.numeric(trainingPtype2[NonNAindex]) #This makes sure you use all the response data for the drug without NA values.
363. #PowerTransform the phenotype if specified.
364. #_______________________________________________________________
365. offset = 0
366. if (powerTransformPhenotype){
367. if (min(trainingPtype4) < 0){ #All numbers must be positive for a powertransform to work, so make them positive.
368. offset <- -min(trainingPtype4) + 1
369. trainingPtype4 <- trainingPtype4 + offset
370. }
371. transForm <- powerTransform(trainingPtype4)[[6]]
372. trainingPtype4 <- trainingPtype4^transForm
373. }
375. #Create the ridge regression model on the training data using pcr (keeping the number of components required for 80% variance).
376. #_______________________________________________________________
377. if (pcr){
379. #There are many ways for pcr in R...here I will use pcr() (principal component regression).
381. #Use pcr to predict for testing data.
382. #_______________________________________________________________
383. train_x<-(t(homData$train)[samps,keepRows]) #samps represent the cell lines that have been filtered, keepRows represents the genes.
384. train_y<-trainingPtype4
386. test_x<-(t(homData$test)[,keepRows])
388. #Remove genes that you have no transcriptome data for (aka columns are filled with 0's).
389. #If you don't, it causes problems in linearRidge().
390. #This is weird, but there were 2 genes in CTRPv2 data that had 0's for one drug (I suppose this occurs depending on sample filtration).
391. #Make sure you remove those same genes from the test data...want to make sure the genes are the same in both train_x and test_x
392. x<-as.vector(colSums(train_x)) #Sum each column.
393. bad<-which(x == 0) #Column/genes indices that are filled with only 0's.
394. if(length(bad) != 0){
395. train_x<-train_x[,-bad]
396. test_x<-data.frame(test_x[,-bad])
397. }
398. #Remove genes with 0 variance. If you don't, this will also cause problems in linearRidge().
399. variance<-c()
400. for(i in 1:ncol(train_x)){
401. variance[i]<-var(as.vector(train_x[,i]))
402. }
403. bi<-which(variance %in% 0) #Bad index...gene has 0 variance.
404. if(length(bi) != 0){
405. train_x<-train_x[,-bi]
406. test_x<-data.frame(test_x[,-bi])
407. }
409. #Check to make sure you have more than 1 training sample for the drug's model.
410. #_______________________________________________________________
411. trainFrame<-try(data.frame(Resp=train_y, train_x), silent = TRUE)
412. if (dim(trainFrame)[1] == 1){ #Make sure you have more than 1 sample.
413. drugs = drugs[-a]
414. message(paste("\n", drugs[a], "is skipped due to insufficient cell lines to fit the model."))
415. next
416. } else {
418. pcr_model<-pcr(Resp~., data=trainFrame, validation='CV')
420. v=cumsum(explvar(pcr_model)) #A vector of all the pcs and their percent of variance.
421. ncomp=min(which(v > percent)) #Identify which pcs will represent 80% variance.
423. #vs[a]<-ncomp
425. if(printOutput) message("\nCalculating predicted phenotype using pcr...")
426. preds<-predict(pcr_model, newdata=test_x, ncomp=ncomp)
428. #You can compute an R^2 value for the data you train on from true and predicted values.
429. #The rsq value represents the percentage in which the optimal model accounts for the variance in the training data.
430. #_______________________________________________________________
431. if (rsq){
433. if (dim(train_x)[1] < 4){ #The code will result in an error if you have 3 samples (which is enough for the model fitting but not when you do a 70/30% split)...
434. message(paste("\n", drugs[a], 'is skipped for R^2 analysis'))
435. }else{
437. #trainFrame<-data.frame(Resp=train_y, train_x)
439. data<-(cbind(train_x, train_y)) #Rows are samples, columns are genes.
441. dt<-sort(sample(nrow(data), nrow(data)*.7)) #sample() randomly picks 70% of rows/samples from the dataset. It samples without replacement.
443. #Prepare the training data (70% of original training data)
444. train_x<-data[dt,]
445. ncol<-dim(train_x)[2]
446. train_y<-train_x[,ncol]
447. train_x<-train_x[,-ncol]
449. #Prepare the testing data (30% of original training data)
450. test_x<-data[-dt,]
451. ncol<-dim(test_x)[2]
452. test_y<-test_x[,ncol]
453. test_x<-test_x[,-ncol]
455. #Remove genes that you have no transcriptome data for aka columns are filled with 0's.
456. #If you don't, it causes problems in linearRidge()
457. #Make sure you remove these same genes from the test data...want to make sure the genes are the same in both.
458. x<-as.vector(colSums(train_x))
459. bad<-which(x == 0)
460. if (length(bad) != 0){
461. train_x<-train_x[,-bad]
462. test_x<-data.frame(test_x[,-bad])
463. }
464. #Remove genes with 0 variance. If you don't, this will also cause problems in linearRidge()
465. variance<-c()
466. for(i in 1:ncol(train_x)){
467. variance[i]<-var(as.vector(train_x[,i]))
468. }
469. bi<-which(variance %in% 0) #Bad index...gene has 0 variance.
470. if (length(bi) != 0){
471. train_x<-train_x[,-bi]
472. }
474. data<-data.frame(Resp=train_y, train_x)
476. pcr_model<-pcr(Resp~., data=data, validation='CV') #Set validation argument to CV.
478. v=cumsum(explvar(pcr_model)) #A vector of all the pcs and their percent of variance.
479. ncomp=min(which(v > percent)) #Identify which pcs will represent 80% variance.
481. pcr_pred<-predict(pcr_model, test_x, ncomp=ncomp)
483. if (printOutput) message("\nCalculating R^2...")
484. sst<-sum((test_y - mean(test_y))^2) #Compute the sum of squares total.
485. sse<-sum((pcr_pred - test_y)^2) #Compute the sum of squares error.
486. rsq_value<-1 - sse/sst #Compute the rsq value.
487. }
488. }
489. }
491. if (report_pc){
492. if (printOutput) message("\nObtaining principal components...")
493. pcs<-coef(pcr_model, comps = ncomp) #comps: numeric, which components to return.
494. dir.create("./calcPhenotype_Output")
495. path<-paste('./calcPhenotype_Output/', drugs[a], '.RData', sep="")
496. save(pcs, file=path)
497. }
499. } else {
501. #Create the ridge regression model on our training data to predict for our actual testing data without pcr.
502. #_______________________________________________________________
503. if(printOutput) message("\nFitting Ridge Regression model...");
505. expression<-(t(homData$train)[samps,keepRows]) #samps represent the cell lines that have been filtered, keepRows represents the genes.
507. test_x<-(t(homData$test)[,keepRows])
509. #Remove genes that you have no transcriptome data for (aka columns are filled with 0's).
510. #If you don't, it causes problems in linearRidge().
511. #This is weird, but there were 2 genes in CTRPv2 data that had 0's for one drug (I suppose this occurs depending on sample filtration).
512. #Make sure you remove those same genes from the test data...want to make sure the genes are the same in both train_x and test_x
513. x<-as.vector(colSums(expression))
514. bad<-which(x == 0) #Column/genes indices that are filled with only 0's.
515. if (length(bad) != 0){
516. expression<-expression[,-bad]
517. test_x<-data.frame(test_x[,-bad])
518. }
520. #Remove genes with 0 variance. If you don't, this will also cause problems in linearRidge().
521. variance<-c()
522. for (i in 1:ncol(expression)){
523. variance[i]<-var(as.vector(expression[,i]))
524. }
525. bi<-which(variance %in% 0) #Bad index...gene has 0 variance.
526. if (length(bi) != 0){ #If there are actually bad indices/genes that had 0 variance across all samples, remove them.
527. expression<-expression[,-bi]
528. test_x<-data.frame(test_x[,-bi])
529. }
531. #Check to make sure you have more than 1 training sample for the drug's model.
532. #_______________________________________________________________
533. trainFrame<-try(data.frame(Resp=trainingPtype4, expression), silent = TRUE)
534. if (dim(trainFrame)[1] == 1){ #Make sure you have more than 1 sample.
535. drugs = drugs[-a]
536. message(paste("\n", drugs[a], "is skipped due to insufficient cell lines to fit the model."))
537. next
538. } else {
540. if(printOutput) message("\nCalculating predicted phenotype...")
542. rrModel<-linearRidge(Resp ~., data=trainFrame)
544. preds<-predict(rrModel, newdata=data.frame(test_x))
546. #You can compute an R^2 value for the data you train on from true and predicted values.
547. #The R^2 value represents the percentage in which the optimal model accounts for the variance in the training data.
548. #_______________________________________________________________
549. if(rsq){
551. if (dim(expression)[1] < 4){ #The code will result in an error if you have 3 samples...this makes sure you have more than 3 samples.
552. #It results in an error because there is a 70/30% split of training data.
553. message(paste("\n", drugs[a], 'is skipped for R^2 analysis'))
554. } else {
555. expression<-(cbind(expression, trainingPtype4))
556. dt<-sort(sample(nrow(expression), nrow(expression)*.7)) #sample() randomly picks 70% of rows/samples from the dataset. It samples without replacement.
558. #Prepare the training data (70% of original training data)
559. train_x<-expression[dt,]
560. ncol<-dim(train_x)[2]
561. train_y<-train_x[,ncol]
562. train_x<-train_x[,-ncol]
564. #Prepare the testing data (30% of original training data)
565. test_x<-expression[-dt,]
566. ncol<-dim(test_x)[2]
567. test_y<-test_x[,ncol]
568. test_x<-test_x[,-ncol]
570. #Remove genes that you have no transcriptome data for (aka columns are filled with 0's).
571. #If you don't, it causes problems in linearRidge().
572. #Make sure you remove those same genes from the test data...want to make sure the genes are the same in both train_x and test_x
573. x<-as.vector(colSums(train_x))
574. bad<-which(x == 0) #Column/genes indices that are filled with only 0's.
575. if (length(bad) != 0){
576. train_x<-train_x[,-bad]
577. test_x<-data.frame(test_x[,-bad])
578. }
580. #Remove genes with 0 variance. If you don't, this will also cause problems in linearRidge().
581. variance<-c()
582. for (i in 1:ncol(train_x)){
583. variance[i]<-var(as.vector(train_x[,i]))
584. }
585. bi<-which(variance %in% 0) #Bad index...gene has 0 variance.
586. if (length(bi) != 0){ #If there are actually bad indices/genes that had 0 variance across all samples, remove them.
587. train_x<-train_x[,-bi]
588. test_x<-data.frame(test_x[,-bi])
589. }
591. trainFrame<-data.frame(Resp=train_y, train_x)
592. rrModel<-linearRidge(Resp ~., data=trainFrame)
594. testFrame<-data.frame(test_x)
596. pred<-predict(rrModel, newdata=testFrame)
598. if(printOutput) message("\nCalculating R^2...")
600. sst<-sum((test_y - mean(test_y))^2) #Compute the sum of squares total.
601. sse<-sum((pred - test_y)^2) #Compute the sum of squares error.
602. rsq_value<-1 - sse/sst #Compute the rsq value.
604. }
605. }
606. }
607. }
609. #If the response variable was transformed (aka powerTransformPhenotype=TRUE), untransform it here.
610. #_______________________________________________________________
611. if(powerTransformPhenotype) {
612. preds <- preds^(1/transForm)
613. preds <- preds - offset
614. }
616. #Find correlation between imputed response and expression of a gene.
617. #Each gene is corrected differently; therefore, it may not be ideal to determine weight or percentage or weight in which each gene contributes to the prediction.
618. #_______________________________________________________________
619. if(cc){ #You can only collect correlations if pcr=FALSE!
620. if(pcr){
621. stop('ERROR: pcr must equal FALSE in order to compute correlations') #It doesn't make sense to compute correlations when the features have changed from genes to pcs.
622. }
624. if(printOutput) message("\nCalculating correlation coefficients...") #This is only relevant if you aren't using pcr.
626. cors_vec<-c() #This vector will store the correlation coefficients.
627. cors_vec2<-c() #This vector will store the p values.
628. matrix<-homData$test[keepRows,] #Matrix of genes x cell lines/cosmic ids.
629. for(d in 1:nrow(matrix)){ #For each gene...
630. cors_vec[d]<-cor.test(as.vector(matrix[d,]), as.vector(preds))$estimate #Compute correlation coefficient for expression of a given gene across patients vs.imputed values for a given drug across patients.
631. cors_vec2[d]<-cor.test(as.vector(matrix[d,]), as.vector(preds))$p.value
632. }
633. #indices<-order(cor, decreasing=TRUE)
634. #ordered_cor<-cor[indices] #Order the correlation coefficients from big to small.
635. #genes<-rownames(matrix)
636. #ordered_genes<-genes[indices]
637. #names(ordered_cor)<-ordered_genes
638. }
640. if(printOutput) message(paste("\nDone making prediction for drug", a, "of", ncol(trainingPtype)))
642. #Store the data in your lists.
643. #_______________________________________________________________
644. DrugPredictions[[a]]<-preds
646. if(rsq){
647. rsqs[[a]]<-rsq_value
648. }
650. if(cc){
651. cors[[a]]<-cors_vec
652. pvalues[[a]]<-cors_vec2
653. }
654. }
655. }
657. #Time to save the data!
658. #_______________________________________________________________
659. #Save drug prediction data to your home directory as a .txt file.
660. names(DrugPredictions)<-drugs
661. DrugPredictions_mat<-do.call(cbind, DrugPredictions)
662. colnames(DrugPredictions_mat)<-drugs
663. rownames(DrugPredictions_mat)<-colnames(testExprData)
664. dir.create("./calcPhenotype_Output")
665. write.csv(DrugPredictions_mat, file="./calcPhenotype_Output/DrugPredictions.csv", row.names = TRUE, col.names = TRUE)
667. #If rsq=TRUE, save R^2 data.
668. if(rsq){
669. names(rsqs)<-drugs
670. rsqs_mat<-do.call(cbind, rsqs)
671. dir.create("./calcPhenotype_Output")
672. write.table(rsqs_mat, file="./calcPhenotype_Output/R^2.txt")
673. }
675. #If CC=TRUE, save correlation coefficient data.
676. if(cc){
677. names(cors)<-drugs
678. cor_mat<-do.call(cbind, cors)
679. rownames(cor_mat)<-rownames(homData$train[keepRows,NonNAindex])
680. colnames(cor_mat)<-drugs
681. dir.create("./calcPhenotype_Output")
682. write.table(cor_mat, file="./calcPhenotype_Output/cors.txt")
684. names(pvalues)<-drugs
685. p_mat<-do.call(cbind, pvalues)
686. rownames(p_mat)<-rownames(homData$train[keepRows, NonNAindex])
687. colnames(p_mat)<-drugs
688. dir.create("./calcPhenotype_Output")
689. write.table(p_mat, file="./calcPhenotype_Output/pvalues.txt")
690. }
692. #print(vs)
693. }

CALCPHENOTYPE.R at commit f7c1944, no license · at the source