①高低風險組某基因表達量差異分析

輸入文件symbol.png

輸入文件risk.png

#引用包
library(limma)
library(ggpubr)
gene="CXCL11"              #基因名稱日杈，需修改
expFile="symbol.txt"     #表達輸入文件
riskFile="risk.txt"      #風險輸入文件
setwd("E:\\research")     #設置工作目錄

#讀取表達輸入文件
rt=read.table(expFile, header=T, sep="\t", check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)

#提取目標基因表達量
data=data[gene,,drop=F]
data=t(data)
rownames(data)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*)", "\\1\\-\\2\\-\\3", rownames(data))
data=avereps(data)

#讀取風險輸入文件
risk=read.table(riskFile, header=T, sep="\t", check.names=F, row.names=1)

#病人風險值和表達數(shù)據(jù)合并
sameSample=intersect(row.names(risk), row.names(data))
risk=risk[sameSample, "risk",drop=F]
data=data[sameSample, , drop=F]
rt=cbind(risk, data)

#設置比較組
#rt[,gene][rt[,gene]>quantile(rt[,gene],0.99)]=quantile(rt[,gene],0.99)
rt$risk=factor(rt$risk, levels=c("low", "high"))
type=levels(factor(rt[,"risk"]))
comp=combn(type, 2)
my_comparisons=list()
for(i in 1:ncol(comp)){my_comparisons[[i]]<-comp[,i]}

#繪制小提琴圖
pdf(file=paste0(gene, ".vioplot.pdf"), width=6, height=5)
ggviolin(rt, x="risk", y=gene, fill = "risk", 
         xlab="Risk", ylab=paste0(gene, " expression"),
         legend.title="Risk",
         palette=c("green", "red"),
         add = "boxplot", add.params = list(fill="white"))+ 
         stat_compare_means(comparisons = my_comparisons,symnum.args=list(cutpoints = c(0, 0.001, 0.01, 0.05, 1), symbols = c("***", "**", "*", "ns")),label = "p.signif")
dev.off()

②藥物敏感性分析

安裝 pRRophetic 包蘸泻，首先要安裝它的依賴包

if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install(c("car", "ridge", "preprocessCore", "genefilter", "sva"))

依賴包安裝成功后碉纳，從本地安裝pRRophetic壓縮包

從本地安裝pRRophetic 包.png

安裝好之后我們需要測試下包能不能正常使用

> library(pRRophetic)
Warning message:
replacing previous import ‘car::Anova’ by ‘genefilter::Anova’ when loading ‘pRRophetic’ 
> set.seed(12345)

輸入文件symbol.png

輸入文件risk.png

#####藥物列表######
#A.443654, A.770041, ABT.263, ABT.888, AG.014699, AICAR, AKT.inhibitor.VIII, AMG.706, AP.24534, 
#AS601245, ATRA, AUY922, Axitinib, AZ628, AZD.0530, AZD.2281, AZD6244, AZD6482, AZD7762, AZD8055, 
#BAY.61.3606, Bexarotene, BI.2536, BIBW2992, Bicalutamide, BI.D1870, BIRB.0796, Bleomycin, 
#BMS.509744, BMS.536924, BMS.708163, BMS.754807, Bortezomib, Bosutinib, Bryostatin.1, BX.795, 
#Camptothecin, CCT007093, CCT018159, CEP.701, CGP.082996, CGP.60474, CHIR.99021, CI.1040, Cisplatin, 
#CMK, Cyclopamine, Cytarabine, Dasatinib, DMOG, Docetaxel, Doxorubicin, EHT.1864, Elesclomol, 
#Embelin, Epothilone.B, Erlotinib, Etoposide, FH535, FTI.277, GDC.0449, GDC0941, Gefitinib, 
#Gemcitabine, GNF.2, GSK269962A, GSK.650394, GW.441756, GW843682X, Imatinib, IPA.3, 
#JNJ.26854165, JNK.9L, JNK.Inhibitor.VIII, JW.7.52.1, KIN001.135, KU.55933, Lapatinib, 
#Lenalidomide, LFM.A13, Metformin, Methotrexate, MG.132, Midostaurin, Mitomycin.C, MK.2206, 
#MS.275, Nilotinib, NSC.87877, NU.7441, Nutlin.3a, NVP.BEZ235, NVP.TAE684, Obatoclax.Mesylate, 
#OSI.906, PAC.1, Paclitaxel, Parthenolide, Pazopanib, PD.0325901, PD.0332991, PD.173074, PF.02341066, 
#PF.4708671, PF.562271, PHA.665752, PLX4720, Pyrimethamine, QS11, Rapamycin, RDEA119, RO.3306, 
#Roscovitine, Salubrinal, SB.216763, SB590885, Shikonin, SL.0101.1, Sorafenib, S.Trityl.L.cysteine, 
#Sunitinib, Temsirolimus, Thapsigargin, Tipifarnib, TW.37, Vinblastine, Vinorelbine, Vorinostat, 
#VX.680, VX.702, WH.4.023, WO2009093972, WZ.1.84, X17.AAG, X681640, XMD8.85, Z.LLNle.CHO, ZM.447439
#####藥物列表######

#引用包
library(limma)
library(ggpubr)
library(pRRophetic)
library(ggplot2)
set.seed(12345)
expFile="symbol.txt"     #表達輸入文件
riskFile="risk.txt"      #風險輸入文件
drug="Rapamycin"         #藥物名稱裙戏，需要修改
setwd("E:\\research")    #設置工作目錄

#讀取表達輸入文件,并對數(shù)據(jù)進行處理
rt = read.table(expFile, header=T, sep="\t", check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
data=data[rowMeans(data)>0.5,]

#刪掉正常樣品
group=sapply(strsplit(colnames(data),"\\-"), "[", 4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2","1",group)
data=data[,group==0]
data=t(data)
rownames(data)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*)", "\\1\\-\\2\\-\\3", rownames(data))
data=avereps(data)
data=t(data)

#預測藥物敏感性
senstivity=pRRopheticPredict(data, drug, selection=1)
senstivity=senstivity[senstivity!="NaN"]
#senstivity[senstivity>quantile(senstivity,0.99)]=quantile(senstivity,0.99)

#讀取風險輸入文件
risk=read.table(riskFile, header=T, sep="\t", check.names=F, row.names=1)

#風險文件和藥物敏感性合并
sameSample=intersect(row.names(risk), names(senstivity))
risk=risk[sameSample, "risk",drop=F]
senstivity=senstivity[sameSample]
rt=cbind(risk, senstivity)

#設置比較組
rt$risk=factor(rt$risk, levels=c("low", "high"))
type=levels(factor(rt[,"risk"]))
comp=combn(type, 2)
my_comparisons=list()
for(i in 1:ncol(comp)){my_comparisons[[i]]<-comp[,i]}

#繪制箱線圖
boxplot=ggboxplot(rt, x="risk", y="senstivity", fill="risk",
                  xlab="Risk",
                  ylab=paste0(drug, " senstivity (IC50)"),
                  legend.title="Risk",
                  palette=c("green", "red")
                 )+ 
    stat_compare_means(comparisons=my_comparisons)
pdf(file=paste0(drug, ".pdf"), width=5, height=4.5)
print(boxplot)
dev.off()