怎样迈出精准医学分析第一步？Xeva帮你打好药物靶向的基础！

#导入包source("imp_functions.R")library(foreach)library(doParallel)library(psych)##plotViolin：使用ggplot2包生成小提琴图。plotViolin <- function(dfx, pltCol, ylim=NULL){  p <- ggplot(data=dfx, aes(x=type, y=ICC , fill=type,color=type))  p <- p + geom_boxplot(width = .001, outlier.size = 2.0, outlier.color = "#525252",                        outlier.shape = 95)  p <- p + geom_violin(trim=F)  p <- p + theme_classic()  p <- p + scale_fill_manual(values=pltCol)  p <- p + scale_color_manual(values=pltCol)      p <- p + stat_summary(fun.data=data_summary, colour = "white", size=0.4)  p <- p + theme(legend.position="none")  p <- p + scale_y_continuous(breaks = seq(-1,1,by = 0.25))  p <- p + theme(axis.text.x = element_text(colour="grey20",size=12, angle=0,face="bold"),                 axis.text.y = element_text(colour="grey20",size=11, face="bold"),                 axis.title.y = element_text(colour="grey20",size=11, face="bold"),                 axis.ticks.length=unit(0.25,"cm"))  p <- p + theme(axis.line.x = element_line(size = 0, colour = "white"),                 axis.title.x=element_blank(),                 axis.ticks.x=element_blank())if(!is.null(ylim))  { p <- p + coord_cartesian(ylim = ylim, expand = F) }return(p)}##computeRankMatFor1Gene函数：计算给定基因在样本中的排名矩阵    computeRankMatFor1Gene <- function(gn, pt, marExp){  findRankOfGene <- function(gn, sampleIds, marExp)  {    rtx <- rep(NA, length(sampleIds)); names(rtx)<- names(sampleIds)for(n in names(rtx))    {s <- sampleIds[n]if(!is.na(s))      {        gval <- sort(marExp[, s])        rtx[n] <- which(names(gval)==gn)      }    }return(rtx)  }  gnRankMat <- matrix(data = NA, nrow = nrow(pt), ncol = ncol(pt))  rownames(gnRankMat) <- rownames(pt); colnames(gnRankMat)<- colnames(pt)for(patientID in rownames(pt))  {        sampleIds <- pt[patientID, ]    gnRankMat[patientID, ] <- findRankOfGene(gn, sampleIds, marExp)  }  gnRankMat <- gnRankMat[!apply(gnRankMat, 1, function(x)all(is.na(x))), ]return(gnRankMat)}##getICCforAllGen函数：用于计算所有基因的ICC（Intraclass Correlation Coefficient）。getICCforAllGenes <- function(pt, allGenes, marExp, NCPU=2, verbose=F){if(verbose==TRUE) { cl<-makeCluster(NCPU, outfile="") } else  { cl<-makeCluster(NCPU)}  registerDoParallel(cl)  allICC <- foreach(i=1:length(allGenes), .export=c("computeRankMatFor1Gene"),                    .packages = c("psych") ) %dopar%                    {                      gn <- allGenes[i]                      gnRankMat <- computeRankMatFor1Gene(gn, pt, marExp)                          v <- NULLif(nrow(gnRankMat)>2)                      {                         v <- psych::ICC(gnRankMat, missing = F)                        v$geneName <- gn                      }                      v                    }  stopCluster(cl)  rtx <- data.frame()for(i in 1:length(allICC))  {if(!is.null(allICC[[i]]))    {      v <- allICC[[i]]$results["Single_raters_absolute",]      v["gene"] <- allICC[[i]]$geneName      rtx <- rbind(rtx, v[c("ICC", "p", "gene")])    }  }  rownames(rtx) <- NULLreturn(rtx)    }## randomizeMatGetICC函数：用于对数据进行随机化处理，并计算随机化后的基因ICC值。randomizeMatGetICC <- function(passMatx, allGenes, marExp, NCPU, verbose=F){  passMatRand <- passMatxfor(i in 1:ncol(passMatRand))  { passMatRand[,i] <- sample(passMatRand[,i]) }for(i in 1:nrow(passMatRand))  { passMatRand[i,] <- sample(passMatRand[i,]) }  passMatRand <- passMatRand[!(apply(passMatRand, 1, function(x) all(is.na(x)))),]  rownames(passMatRand) <- paste0("r", 1:nrow(passMatRand))  allICCRand <- getICCforAllGenes(passMatRand, allGenes, marExp, NCPU=NCPU,                                  verbose=verbose)return(allICCRand)}## 数据加载和准备NCPU = max(1, detectCores()-1 )cat(sprintf("n##---- number of CPU %d ------n", NCPU))passage <- readRDS("../data/passage_data.Rda")dp <- pData(passage$data)mar <- passage$dataallGenes <- rownames(mar)marExp <- Biobase::exprs(mar)pheno <- Biobase::pData(mar)pheno$passage <- paste0("P", pheno$passage)passMat <- passage$passageif(recompute==TRUE){  cat(sprintf("recomputing values.... it may take long timen"))  ICCv <- list()##-------------------------------------------------------##------------------For all tissue type -----------------##-------------------------------------------------------  allICC <- getICCforAllGenes(passMat, allGenes, marExp, NCPU=NCPU, verbose=F)  allICC$type <- "All tissue"  allICCRand <- randomizeMatGetICC(passMat, allGenes, marExp, NCPU=NCPU, verbose=F)  allICCRand$type <- "Null (All tissue)"  ICCv[["All tissue"]] <- rbind(allICC, allICCRand)##-------------------------------------------------------##------------------By Tissue Type ----------------------  tissueName <- c("breast", "lung",  "pancreas", "skin", "large_intestine")for(tt in tissueName)  {      ptt <- passMat[pheno[pheno$tumor.type==tt, "patient.id"], , drop=F]if(nrow(ptt)>8)      {        cat(sprintf("##--- running for %s -----n", tt))        allICCtt <- getICCforAllGenes(ptt, allGenes, marExp, NCPU=NCPU, verbose=F)        allICCtt$type <- tt        allICCttRand <- randomizeMatGetICC(ptt,allGenes,marExp,NCPU=NCPU,verbose=F)            allICCttRand$type <- sprintf("Null (%s)", tt)        ICCv[[tt]] <- rbind(allICCtt, allICCttRand)      }##储存结果  saveRDS(ICCv, "../results/ICC_for_genes.Rda")  }} else  {ICCv <- readRDS("../data/ICC_for_genes.Rda")}##绘制图片pltColo <- c('#e41a1c', '#8dd3c7','#bebada','#fb8072','#80b1d3','#fdb462')names(pltColo)<-c("All tissue","breast","lung","pancreas","skin","large_intestine")nullCol <- "#4d4d4d"pdf("../results/figure-3.pdf", width = 5.15, height = 2.15)for(tt in names(pltColo)){      dfx <- ICCv[[tt]]  pltCol <- c(pltColo[tt], nullCol); names(pltCol) <- c(tt, sprintf("Null (%s)", tt))  dfx$type <- factor(as.character(dfx$type), levels = names(pltCol))  ylim <- c(-0.25,1)  p <- plotViolin(dfx, pltCol, ylim = ylim)print(p)}dev.off()