①高低风险组某基因表达量差异分析

输入文件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)

#病人风险值和表达数据合并
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")    #设置工作目录

#读取表达输入文件,并对数据进行处理
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()