看nature文章是如何设计和使用普通转录组数据

你好啊，我是Christine，这里是生信菜鸟团周一数据挖掘专栏，欢迎和我们一起学习 ~

本周我们学习的文章是Glucocorticoids promote breast cancer metastasis，2019年3月发表在Nature上，文章里面的NGS数据并不多，也就是17只小鼠的转录组测序而已！

肿瘤的转移受到多种因素的因素的影响，比如经常提到的肿瘤异质性，这篇文章提出糖皮质激素也能促进转移，提示临床上癌症的激素治疗需要慎重考虑。

文章框架：

1. 作者将17名来自肿瘤患者的异种移植物或细胞系植入免疫缺陷小鼠的乳腺中 ，从转录组数据发现，转移灶的糖皮质激素受体（GR）比原发肿瘤更多地被激活

2. 动物实验验证GR与转移的相关性

3. 转录组，蛋白质组学和磷酸化蛋白质组学等手段发现糖皮质激素受体激活多种转移相关通路和ROR1激酶的表达

4. 实验验证敲除ROR1后减少转移，延长生存

需要重复的图是：

Fig1b,c：MDA-MB 231模型中的肿瘤及配对转移样本的表达数据主成分分析；原发及转移样本的差异基因热图

Fig4a：乳腺癌转移样本中GR激活信号与ROR1的相关系数热图

转录组数据PCA

取网页下载GSE124817表达量信息，批量读入R，形成表达矩阵，PCA分析检查根据转录组数据能否将不同组样本分开，同组样本间一致性是否更高。

rm(list = ls())
options(stringsAsFactors = F)
#读入表达信息
dir <- "./GSE124817_RAW/"
files <- list.files(path = dir)
expr <- lapply(files,
               function(x){
                 # x <- files[1]
                 expr <- read.csv(file = file.path(dir,x), header = T, stringsAsFactors = F)[,c(3,6)]
                 return(expr)
               })
df <- do.call(cbind, expr)
df <- na.omit(df)
rownames(df) <- df[,1]
# 只保留一列样本名
df <- df[,seq(2,ncol(df))]
expr <- df
save(expr,file = "./GSE76250_expr.Rdata")

## 主成分分析
rm(list = ls())
load("./GSE76250_expr.Rdata")
# 提取MDAMB231模型的表达矩阵
expr <- expr[,grep("MDAMB231",colnames(expr))]
colnames(expr) <- substr(colnames(expr),10,nchar(colnames(expr)))
# 去掉Spleen样本
expr <- expr[,grep("Spleen",colnames(expr),invert = T)]
# 构建group_list
library(stringr)
group_list <- str_split(colnames(expr),"_",simplify = T)[,1]
save(expr,group_list,file = "./expr_processed.Rdata")

## 主成分分析
library(FactoMineR)
library(factoextra)
# 表达矩阵归一化
library(edgeR)
PCA_matrix <- cpm(expr,log = T)
# 取标准差最大的1000个基因
genes <- sort(apply(PCA_matrix,1,sd),decreasing = T)[1:1000]
PCA_matrix <- PCA_matrix[names(genes),]
df_pca <- PCA(t(PCA_matrix),graph = F)
fviz_pca_ind(df_pca,col.ind = group_list)

多组样本差异分析

原发肿瘤和不同转移样本间的差异表达分析

rm(list = ls())
load("./expr_processed.Rdata")

# 用edgeR做差异分析
library(edgeR)
d <- DGEList(counts=expr,group=factor(group_list))
keep <- rowSums(cpm(d)>1) >= 2 #筛选cpm值
table(keep)
d <- d[keep, , keep.lib.sizes=FALSE]
d$samples$lib.size <- colSums(d$counts)
d <- calcNormFactors(d)
#d$samples

# 构建分组矩阵
design <- model.matrix(~0+factor(group_list))
rownames(design) <- colnames(d)
colnames(design)<-levels(factor(group_list))
dge <- estimateGLMCommonDisp(d,design)
dge <- estimateGLMTrendedDisp(dge, design)
dge <- estimateGLMTagwiseDisp(dge, design)

fit <- glmFit(dge, design)

# 构建多组样本的比较矩阵
cont_matrix <- makeContrasts(CTC_Tumor=CTC-Tumor, Liver_Tumor=Liver-Tumor, Lungs_Tumor=Lungs-Tumor,
                             levels=design)
lrt <- glmLRT(fit,contrast = cont_matrix)
#提取结果
nrDEG <- topTags(lrt, n=nrow(dge))
nrDEG <- as.data.frame(nrDEG)
DEG_edgeR <- na.omit(nrDEG)

筛选差异表达基因，画热图：

DEG_edgeR$change <- ifelse(DEG_edgeR$FDR < 0.05 & 
                             (abs(DEG_edgeR$logFC.CTC_Tumor) >= 2 |
                             abs(DEG_edgeR$logFC.Liver_Tumor >= 2) |
                             abs(DEG_edgeR$logFC.Lungs_Tumor >= 2)),
                           "choose","Not")
table(DEG_edgeR$change) 
#340个差异基因，文章是280个，可能和我们的logCPM筛选不同相关吧

## 画热图
choose_genes <- rownames(DEG_edgeR[DEG_edgeR$change == "choose",]) 
library(pheatmap)
choose_matrix <- expr[choose_genes,]
#将表达矩阵归一化，中心化，标准化
choose_matrix <- t(scale(t(log2(choose_matrix+1))))   
choose_matrix[1:4,1:4]
boxplot(choose_matrix[,1:4])
###调整一些特殊的表达值
choose_matrix[choose_matrix < -2] = -2
choose_matrix[choose_matrix > 2] = 2

pheatmap(choose_matrix,show_rownames = F, 
         color = colorRampPalette(c("blue","black","yellow"))(50))

GR激活信号和ROR1基因表达的相关系数热图

为了验证BRCA转移样本中GR激活信号和ROR1表达量之间的相关性，作者画了它们的相关关系热图，根据文献中的描述，所用的转移样本来自https://met500.path.med.umich.edu/，然而我在下载页面只看到了突变和拷贝数变异数据。

这里先用TCGA-BRCA的转移样本数据代替一下：

rm(list = ls())
# 从UCSC-XENA下载BRCA表达矩阵
# 单位是log2(RSEM+1)
expr <- read.delim("./HiSeqV2.gz",header = T, stringsAsFactors = F)
rownames(expr) <- expr[,1]
expr_mat <- expr[,substr(colnames(expr),14,15)=="06"] # 7个转移样本

## 对表达矩阵归一化
# 根据文献方法的描述：先去掉编码线粒体和核糖体蛋白的基因，再归一化文库大小
# 可能是作者用的数据在建库时引入了很多核糖体和线粒体基因，过滤后只剩18,115个基因
# 预处理后的TCGA RNA-seq是不需要这一步的，这里就当学习一下过滤方法
expr_mat <- 2^expr_mat - 1 #先还原为count

grep("^MT-",rownames(expr),value = T) #没有线粒体蛋白基因

#以RPL*,RPS*,MRPL*,MRPS*开头的基因为核糖体蛋白基因
# Ribosomal Protein Gene Database：http://ribosome.med.miyazaki-u.ac.jp/
ribo_g <- grep("^M*RP[SL]",rownames(expr),value = T) # 199个
expr_mat <- expr_mat[!(rownames(expr_mat) %in% ribo_g),]

library(edgeR)
expr_mat <- cpm(expr_mat,log = T) #校正文库大小，归一化

# 挑选基因画相关系数热图
# 把GR基因表达量加起来，形成RG signature (SGR)
SGR_g <- c("FN1","KLF9", "ANKRD1", "MT2A", "VIM", "SNAI2", "POU5F1", "ID3")
SRG <- colSums(expr_mat[SGR_g,])
choose_matrix <- rbind(expr_mat[c(SGR_g,"ROR1","NR3C1","WNT5A"),],SRG)

cor_mat <- cor(t(choose_matrix))
library(pheatmap)
pheatmap(cor_mat,main = "All metastasis")

虽然数据不同，但趋势和文章是类似的。

总结一下：

1. 对于熟悉我们专栏的朋友来说，本周的数据分析内容是比较简单的，所有的图表之前都做过。同时对多组样本进行差异分析时，构建比较矩阵和平时有一点点的不同。

2. 作者在对样本进行文库校正前先去除了编码核糖体和线粒体蛋白的基因，难道是这个数据本身建库时的问题？单细胞转录组经常会做这样的过滤。TCGA预处理好的RNA-seq数据是没有线粒体蛋白基因的，只有少量的核糖体蛋白基因。

3. 我又一次没有找到作者提到的原数据，不过事实证明，Nature文章的结论即便换个数据集也是同样很靠谱的。

如果你也想完成这样的数据分析，欢迎选择我们全国巡讲的的生物信息学爆款入门线下3天课程！

号外：中秋节广州3天入门课程报名马上截止：（中秋节一起来学习！）全国巡讲第16站-广州（生信入门课加量不加价）