TCGA正常血液样本中PD-L1与BRCA1和NBS1的表达量相关性（没找到原数据，但重复出了结果）

你好啊，欢迎来到周一数据挖掘专栏，我是新学徒Christine，很高兴和你一起学习！

以后将由我负责带领大家进行数据挖掘实战训练。

菜鸟团（周一数据挖掘专栏）成果展 （如果大家感兴趣前面学徒的教程，可以收藏学习）

本周我们学习一篇与肿瘤免疫治疗明星分子PD-L1相关的文章：PD-L1 (B7-H1) Competes with the RNA Exosome to Regulate the DNA Damage Response and Can Be Targeted to Sensitize to Radiation or Chemotherapy，2019年4月30号于Molecular Cell上online。

背景介绍

PD-L1(Programmed death ligand 1，又称B7-H1)是一个免疫检验点蛋白，在多种免疫细胞和肿瘤细胞中表达。PD-L1与其受体PD-1的结合后会抑制T细胞的活性，是正常细胞避免自体免疫的保护机制，也是肿瘤细胞发生免疫逃逸的重要原因。目前，FDA已经批准了Durvalumab, Atezolizumab 等抑制PD-L1与PD-1结合的抗癌药。

文章思路

PD-L1与PD-1的结合发生在细胞外，作者想知道细胞内大量存在的PD-L1有什么意义？

• 首先，敲减PD-L1，发现DNA损伤增加，表现在NBS1和BRCA1的表达降低，究其原因是mRNA稳定性降低了

• 接着，证明PD-L1通过与RNA外切体竞争结合保护mRNA免受降解

• 然后，探究PD-L1在全基因组范围内对RNA稳定性的调控：用到了RNA免疫共沉淀（RIP-seq）和RNA-seq

• 最后，找到一个抗体H1A，能抑制PD-L1对DNA损伤的保护作用，增强放疗敏感性。

img

RIP-seq：RNA免疫共沉淀与高通量测序结合的技术，通过免疫沉淀目标蛋白来捕获蛋白体内结合的RNA，将捕获的RNA进行高通量测序，可获得体内RNA结合蛋白与RNA靶标的结合模式，帮助阐述目的蛋白对基因表达调控的分子机制。

PD-L1与BRCA1和NBS1表达量的相关性

这次我们复现文中一个很简单的图Fig2H：用TCGA中407个血液来源的正常样本证明PD-L1与BRCA1和NBS1的表达量相关性很高。

Fig2H

Step1. 下载pan-cancer的基因表达数据

从UCSCXena上下载校正了批次效应且标准化后的pan-cancer表达数据：EB++AdjustPANCAN_IlluminaHiSeq_RNASeqV2.geneExp.xena

Step2. 提取3个基因的表达矩阵

因为pan-cancer的表达矩阵很大，就不往R里读了，直接用linux命令提取：

PD-L1对应的基因名是CD274，NBS1对应的基因名是NBN

zcat EB++AdjustPANCAN_IlluminaHiSeq_RNASeqV2.geneExp.xena.gz |grep -E 'sample|BRCA1|NBN|CD274' > expr.txt

Step3. 找到数据中血液来源的正常样本，画图

TCGA样本的barcode中第14，15位数字代表样本类型，具体可以看官方的说明文档(https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables/sample-type-codes)比较常用的有：

• 01 原发样本

• 02 复发样本

• 06 转移样本

• 10 血液来源的正常样本

• 11 正常组织样本

然而，这个表达矩阵中并没有编号10的样本！！！

首先猜测是不是UCSC在整理数据时删掉了？

我用GDC API没有查到“Blood Derived Normal”样本的"RNA-seq"文件。难道是表达量还有别的实验手段？

我又查询了全部365,463个文件的实验方法，发现表达量的检测方法真的只有“RNA-seq”这一种。

难道这个数据就没有存放在GDC中？

总之最终我也没有找到。

如果你知道该怎么获得这个数据，麻烦告诉我一下吧 ~

先用正常样本画个图看看：

rm(list = ls())
options(stringsAsFactors = F)

# 读入表达数据
expr <- read.table("./expr.txt",sep = "\t",header = T, stringsAsFactors = F)
expr[1:3,1:4]
rownames(expr) <- expr$sample
expr <-expr[,-1]
# 查看样本的组成
table(substr(colnames(expr),14,15)) #发现没有血液来源的正常样本10 ?!
expr <- expr[,substr(colnames(expr),14,15)==11] #取正常样本
#画图
library(ggstatsplot)
df <- as.data.frame(t(expr))
colnames(df) <- c("BRCA1","PD_L1","NBS1")
p1 <- ggscatterstats(df, x="PD_L1",y="NBS1",marginal = F)
p2 <-  ggscatterstats(df, x="PD_L1",y="BRCA1",marginal = F)
library(ggpubr)
ggarrange(p1,p2,ncol=2,nrow=1,labels=c("A","B"))

因为样本用得不对，图中并未看到PD-L1与BRCA1和NBS1的表达相关。

GTEx数据

我觉得自己已经没招了，只能求助Jimmy老师了，正在西安“逛&吃”的老师给了很简单的一句话，“既然用的正常样本，考虑一下GTEx数据”。所以，以下是GTEx的表达数据，操作步骤和上面是一样的，但结果很惊喜哦~

Genotype-Tissue Expression (GTEx)：收集了大量健康人体器官mRNA测序的数据库，根据基因型-组织表达可进行跨器官的eQTL研究。

Step1.GTEx网站下载数据

网址： https://gtexportal.org/

基因表达矩阵：GTEx_Analysis_2016-01-15_v7_RNASeQCv1.1.8_gene_tpm.gct.gz

样本信息：GTEx_v7_Annotations_SampleAttributesDS.txt

Step2. 提取3个基因的表达矩阵

 zcat GTEx_Analysis_2016-01-15_v7_RNASeQCv1.1.8_gene_tpm.gct.gz |grep -E 'Name|BRCA1|NBN|CD274' > GTEx_expr.txt

Step3. 找到血液样本，画图

rm(list = ls())
options(stringsAsFactors = F)

## 读取表达矩阵
GTEx_expr <- read.delim("C:/Users/lx/Desktop/WEEK2/GTEx_expr.txt",header=T,stringsAsFactors = F)
GTEx_expr[1:3,1:5]
rownames(GTEx_expr) <- GTEx_expr$Description

## 读入GTEx的表型信息，筛选血液样本
phe <- read.delim("C:/Users/lx/Desktop/WEEK2/GTEx_v7_Annotations_SampleAttributesDS.txt",
                  header = T,stringsAsFactors = F)
table(phe$SMTS) 
s <- phe[phe$SMTS == "Blood",1] #筛选出血液样本
s <- gsub("-",".",s)

#画图
s <- intersect(s,colnames(GTEx_expr))
expr <- GTEx_expr[,s] #537个样本
expr <- log2(expr)
library(ggstatsplot)
df <- as.data.frame(t(expr))
colnames(df) <- c("BRCA1","PD_L1","NBS1")
p1 <- ggscatterstats(df, x="PD_L1",y="NBS1",marginal = F)
p2 <-  ggscatterstats(df, x="PD_L1",y="BRCA1",marginal = F)
library(ggpubr)
ggarrange(p1,p2,ncol=2,nrow=1)

果然，在正常血液样本中，PD-L1与BRCA1和NBS1表达的相关性确实很高，相关系数甚至高于原文。应该是因为正常样本中PD-L1主要在淋巴细胞中表达，只有测血液样本才能看出来意义。

另外，文中的RIP-seq和RNA-seq数据都已经上传到了GEO数据库（GSE128613和GSE128742），但是目前还无法访问，预计2019年6月1日开放，拿到表达数据后，参照学习菜鸟团（周一数据挖掘专栏）成果展，下图A和B的热图也是可以做的，也可以自己做差异分析、GO富集分析，感兴趣的同学到时候可以试一试~

Fig4

一点心得：

1. 数据分析中生物背景还是非常重要的，比如这个正常血液样本中的PD-L1表达，就不能随便拿正常组织样本替代。

2. 有些生信分析的虽然看起来很简单，但自己做起来还是会遇到各种意想不到的问题，这个时候需要换个思路，当然最好还是有“大神老师”指点一下啦！

3. 有哪位小伙伴知道TCGA中”Blood Derived Normal“样本的基因表达数据到底该去哪里找吗？期待你的答案~

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