①免疫细胞浸润与风险值相关性分析

下载infiltration_estimation_for_tcga

下载界面.png
输入文件infiltration_estimation_for_tcga.csv.png
输入文件risk.png

#引用包
library(limma)
library(scales)
library(ggplot2)
library(ggtext)
riskFile="risk.txt"      #风险输入文件
immFile="infiltration_estimation_for_tcga.csv"     #免疫细胞浸润文件
setwd("D:\\biowolf")     #设置工作目录

#读取风险输入文件
risk=read.table(riskFile, header=T, sep="\t", check.names=F, row.names=1)

#读取免疫细胞浸润文件
immune=read.csv(immFile, header=T, sep=",", check.names=F, row.names=1)
immune=as.matrix(immune)
rownames(immune)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*)", "\\1\\-\\2\\-\\3", rownames(immune))
immune=avereps(immune)

#对风险文件和免疫细胞浸润文件取交集，得到交集样品
sameSample=intersect(row.names(risk), row.names(immune))
risk=risk[sameSample, "riskScore"]
immune=immune[sameSample,]

#对风险打分和免疫细胞进行相关性分析
x=as.numeric(risk)
outTab=data.frame()
for(i in colnames(immune)){
    y=as.numeric(immune[,i])
    corT=cor.test(x, y, method="spearman")
    cor=corT$estimate
    pvalue=corT$p.value
    if(pvalue<0.05){
        outTab=rbind(outTab,cbind(immune=i, cor, pvalue))
    }
}
#输出相关性结果
write.table(file="corResult.txt", outTab, sep="\t", quote=F, row.names=F)

#绘制气泡图
corResult=read.table("corResult.txt", head=T, sep="\t")
corResult$Software=sapply(strsplit(corResult[,1],"_"), '[', 2)
corResult$Software=factor(corResult$Software,level=as.character(unique(corResult$Software[rev(order(as.character(corResult$Software)))])))
b=corResult[order(corResult$Software),]
b$immune=factor(b$immune,levels=rev(as.character(b$immune)))
colslabels=rep(hue_pal()(length(levels(b$Software))),table(b$Software))     #定义颜色
pdf(file="cor.pdf", width=10, height=6)       #保存图片
ggplot(data=b, aes(x=cor, y=immune, color=Software))+
    labs(x="Correlation coefficient",y="Immune cell")+
    geom_point(size=4.1)+
    theme(panel.background=element_rect(fill="white",size=1,color="black"),
          panel.grid=element_line(color="grey75",size=0.5),
          axis.ticks = element_line(size=0.5),
          axis.text.y = ggtext::element_markdown(colour=rev(colslabels)))
dev.off()

②高低风险组免疫细胞浸润差异分析

输入文件infiltration_estimation_for_tcga.csv.png
输入文件risk.png

#引用包
library(limma)
library(ggpubr)
riskFile="risk.txt"      #风险输入文件
immFile="infiltration_estimation_for_tcga.csv"     #免疫细胞浸润文件
setwd("E:\\research")     #设置工作目录

#读取风险输入文件
risk=read.table(riskFile, header=T, sep="\t", check.names=F, row.names=1)

#读取免疫浸润文件
immune=read.csv(immFile, header=T, sep=",", check.names=F, row.names=1)
immune=as.matrix(immune)
rownames(immune)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*)","\\1\\-\\2\\-\\3",rownames(immune))
immune=avereps(immune)

#病人风险值和免疫细胞合并
sameSample=intersect(row.names(risk), row.names(immune))
risk=risk[sameSample, "risk", drop=F]
immune=immune[sameSample,]
data=cbind(risk, immune)

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

#高低风险组免疫差异分析
for(i in colnames(data)[2:ncol(data)]){
    #绘制箱线图
    boxplot=ggboxplot(data, x="risk", y=i, fill="risk",
                      xlab="Risk",
                      ylab=i,
                      legend.title="Risk",
                      palette=c("green", "red")
                      )+ 
            stat_compare_means(comparisons=my_comparisons)
    wilcoxTest=wilcox.test(data[,i] ~ data[,"risk"])
    if(wilcoxTest$p.value<0.05){
        j=gsub("/", "-", i)
        pdf(file=paste0(j, ".pdf"), width=5, height=4.5)
        print(boxplot)
        dev.off()
    }
}