公共数据辅助乳腺癌的免疫治疗机制研究
自我介绍：
大家好啊，我是菜鸟团的新成员Christine，目前算是生信的业余学习者，名副其实的菜鸟一枚，尝试接手最优秀学员的周一数据挖掘专场，珠玉在前，激励我好好努力，也请大家多多指教 ~
菜鸟团（周一数据挖掘专栏）成果展 （如果大家感兴趣前面学徒的教程，可以收藏学习）
文章简介
本次学习的文章是2019年3月发表于cell research的LSECtin on tumor-associated macrophages enhances breast cancer stemness via interaction with its receptor BTN3A3
LSECtin：一种跨膜蛋白，在肿瘤相关的巨噬细胞中高表达，能增强乳腺癌的细胞干性
BTN3A3：LSECtin的受体
这是一篇很传统的实验生物学文章，证明了乳腺癌中LSECtin-BTN3A3 axis与癌细胞干性相关，有望成为乳腺癌免疫治疗的靶标，最后用到了公共数据验证，给结论锦上添花。
文章的大致思路如下：
肿瘤微环境中巨噬细胞LSECtin高表达，增强肿瘤细胞干性
筛选LSECtin的受体，得到BTN3A3并验证
体内和体外证明LSECtin-BTN3A3 axis及其对乳腺癌细胞干性的影响
公共数据库中验证，也就是接下来我将重复的Fig7a, c, e
PartI. TCGA-BRCA数据
BTN3A3基因在TNBC和Non-TNBC中的表达（图A）
Step1.下载TCGA-BRCA的数据
UCSC Xena网站上直接下载：
基因表达数据：已经是log2(norm_count+1)处理好的表达矩阵了
临床信息
前面学徒多次讲解了下载方法，这里不再赘述，自行查看教程：菜鸟团（周一数据挖掘专栏）成果展
Step2. 提取TNBC样本
TNBC：triple-negative breast cancer，是指癌组织免疫组织化学检查结果为雌激素受体（ER）、孕激素受体（PR）和原癌基因Her-2均为阴性的乳腺癌，预后较其它类型更差。
备注：这里使用3个受体的阴性来判断乳腺癌患者为TNBC，可能并不是最精确的
rm(list = ls())
options(stringsAsFactors = F)
# 读取临床表型数据
phe <- read.table("./raw_data/BRCA_clinicalMatrix.gz", sep = "\t", header = T, stringsAsFactors = F)
colnames(phe) #挑选需要的表型信息，这里选择了ER,PR,HER2状态，生存状况及生存时间
phe <- phe[,c(1,8,10,22,32,33)]
phe$sampleID <- gsub("-",".",phe$sampleID) #使样本名与之后的表达矩阵保持一致
# 筛选TNBC样本
TNBC_phe <- phe[phe$ER_Status_nature2012 == "Negative" &
phe$HER2_Final_Status_nature2012 == "Negative" &
phe$PR_Status_nature2012 == "Negative",] #123个
Non_TNBC_phe <- phe[phe$ER_Status_nature2012 == "Positive" |
phe$HER2_Final_Status_nature2012 == "Positive" |
phe$PR_Status_nature2012 == "Positive",] #648个
# 只保留具有ER,PR,HER2信息的样本
phe <- phe[phe$sampleID %in% c(TNBC_phe$sampleID,Non_TNBC_phe$sampleID),] #771个
phe$TNBCtype <- ifelse(phe$sampleID %in% TNBC_phe$sampleID,"TNBC","Non_TNBC") #添加TNBC信息
Step3. 读取基因表达信息
###读入表达数据
expr <- read.table("./raw_data/HiSeqV2.gz",sep = "\t",header = T, stringsAsFactors = F)
expr[1:4,1:4]
rownames(expr) <- expr[,1]
expr <- expr[,-1]
grep("btn3a3",rownames(expr),ignore.case = T)#该基因位于20274行
df <- as.data.frame(t(expr[20274,]))
Step4. 画图：BTN3A3基因在TNBC和Non-TNBC中的表达量
# 文章中barplot
phe$BTN3A3 <- df[phe$sampleID,1] #将基因表达量添加到表型数据中
library(ggpubr)
ggbarplot(phe,x="TNBCtype",y="BTN3A3",fill = "TNBCtype",
add = "mean_sd",
error.plot = "errorbar") +
stat_compare_means(method = "t.test")
#boxplot结果更清楚
ggboxplot(phe,x="TNBCtype",y="BTN3A3",color = "TNBCtype",
add = "jitter") +
stat_compare_means(method = "t.test")
结果和文章的图不大一样，我们的样本筛选结果是有出入的，而我也无法得知作者是怎么选的，不过从boxplot还是能看出2组间差异的：
BTN3A3基因表达的生存分析（图C）
library(survival)
library(survminer)
#TNBC
#按文章的描述，表达量超过60%的记为High
TNBC_phe$group <- ifelse(TNBC_phe$BTN3A3 > quantile(TNBC_phe$BTN3A3,0.6),'high',"low")
TNBC_phe$OS.time <- TNBC_phe$OS.time/30 #单位改为月
# 生存曲线图
sfit1 <- survfit(Surv(OS.time, OS) ~ group, data=TNBC_phe)
ggsurvplot(sfit1, conf.int=F, pval=TRUE,risk.table =TRUE) +
labs(title = "TNBC")
#Non_TNBC
Non_TNBC_phe$group <- ifelse(Non_TNBC_phe$BTN3A3 > quantile(Non_TNBC_phe$BTN3A3,0.6),'high',"low")
Non_TNBC_phe$OS.time <- Non_TNBC_phe$OS.time/30 #单位改为月
sfit2 <- survfit(Surv(OS.time, OS) ~ group, data=Non_TNBC_phe)
ggsurvplot(sfit2, conf.int=F, pval=TRUE,risk.table =TRUE) +
labs(title = "Non-TNBC")
呃，结果和文章的图完全不同，都没有显著影响生存，主要应该是样本的选择不同，我很奇怪的是，文章中的样本在120个月时居然几乎全部出现结局？
生存分析就是一坨浆糊，让我很生气！
PartII. GEO数据
验证LSECtine与TAM相关基因表达量间的关系（图E）
TAMs: tumor-associated macrophages，作者选了CD68 ,CSF1R ,IL10这3个基因
文章用到了GSE76250，165个TNBC样本，33个配对正常样本，检测肿瘤特异的mRNA和lncRNA
Step1. 下载芯片数据
rm(list = ls()) ## 魔幻操作，一键清空~
options(stringsAsFactors = F)#在调用as.data.frame的时，将stringsAsFactors设置为FALSE可以避免character类型自动转化为factor类型
### 下载芯片数据
# https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE76250
f='./Rdata/GSE76250_eSet.Rdata'
library(GEOquery)
if(!file.exists(f)){
gset <- getGEO('GSE76250', destdir=".",
AnnotGPL = F, ## 注释文件
getGPL = F) ## 平台文件
save(gset,file=f) ## 保存到本地
}
load('./Rdata/GSE76250_eSet.Rdata') ## 载入数据
### 查看一下数据
class(gset) #查看数据类型
length(gset) #此处只有1个芯片平台
class(gset[[1]])
gset #70523 features, 198 samples，芯片平台是GPL17586
Step2.提取TNBC样本
a=gset[[1]]
dat=exprs(a) #a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
dim(dat)#看一下dat这个矩阵的维度
dat[1:4,1:4] #查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列
pd=pData(a) #通过查看说明书知道取对象a里的临床信息用pData
## 挑出TNBC样本
colnames(pd)
table(pd$source_name_ch1) #找到breast_TNBC样本
samples <- rownames(pd[pd$source_name_ch1 == "breast_TNBC",])
Step3.提取4个基因的表达矩阵
### 先要找到基因对应的探针ID
# 下载注释文件
if(F){
library(GEOquery)
gpl <- getGEO('GPL17586', destdir="./raw_data")
colnames(Table(gpl)) # 找到需要的列
probe2gene=Table(gpl)[,c(1,8)]
head(probe2gene)
#提取gene symbol
library(stringr)
probe2gene$gene_assignment <- str_split(probe2gene$gene_assignment,"//",simplify = T)[,2]
probe2gene$gene_assignment <- trimws(probe2gene$gene_assignment) #去掉字符串两边的空格
head(probe2gene)
save(probe2gene,file='probe2gene_GPL17586.Rdata')
}
# 提取四个基因
# LSECtin的基因名是Clec4g
ids <-probe2gene[probe2gene$gene_assignment %in% c("CD68","CSF1R","IL10","CLEC4G"),]
#提取4个基因的表达矩阵
expr <- dat[ids$ID,samples]
rownames(expr) <- ids$gene_assignment
Step4.画散点图
df <- as.data.frame(t(expr))
colnames(df)[4] <- "LSECtin" #再把基因名改回LSECtin
#画图
library(ggstatsplot)
p1 <- ggscatterstats(df,x=LSECtin,y=CD68,marginal = F)
p2 <- ggscatterstats(df,x=LSECtin,y=CSF1R,marginal = F)
p3 <- ggscatterstats(df,x=LSECtin,y=IL10,marginal = F)
library(ggpubr)
ggarrange(p1,p2,p3,ncol=3,nrow=1,labels=c("A","B","C"))
一点心得：
在实验为主的文章中用公共数据验证已经快成为“标配”了，多学习一点还是很有必要的 ~
公共数据生存分析的样本选择似乎是个很有学问的事情… …
最后，谢谢你的支持，欢迎和我一起学习！