公共数据辅助乳腺癌的免疫治疗机制研究

自我介绍：

大家好啊，我是菜鸟团的新成员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与癌细胞干性相关，有望成为乳腺癌免疫治疗的靶标，最后用到了公共数据验证，给结论锦上添花。

文章的大致思路如下：

1. 肿瘤微环境中巨噬细胞LSECtin高表达，增强肿瘤细胞干性

2. 筛选LSECtin的受体，得到BTN3A3并验证

3. 体内和体外证明LSECtin-BTN3A3 axis及其对乳腺癌细胞干性的影响

4. 公共数据库中验证，也就是接下来我将重复的Fig7a, c, e

Fig.7

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组间差异的：

备注：这里的图，其实是Y坐标轴跟作者不一样，所以导致图并没有原文作者那样明显的差异表达直观感受，大家可以思考一下如何修改绘图代码，达到作者的高度！

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"))


画的最好的应该就是这幅图了，而且基本上跟作者差不多的相关系数。

一点心得：

1. 在实验为主的文章中用公共数据验证已经快成为“标配”了，多学习一点还是很有必要的 ~

2. 公共数据生存分析的样本选择似乎是个很有学问的事情… …

最后，谢谢你的支持，欢迎和我一起学习！