三阴性乳腺癌表达数据分析笔记之PAM50

生信技能树

发布于 2020-10-26 10:45:54

2.5K0

发布于 2020-10-26 10:45:54

文章被收录于专栏：生信技能树生信技能树

学徒和学员已经陆续出师，是时候把生信技能树的舞台交给后辈了！

下面是学徒写的《GEO数据挖掘课程》的配套笔记（第6篇）

PAM50的介绍

在临床实践中，就需要HR阳性，HER2阴性乳腺癌的预后和预测模型，为了实现这些目的，基因检测，如OncotypeDx和PAM50, 已经被开发并在多个临床实验中得以验证

基于基因在癌症中的表达，乳腺癌可以分为5个固有的分子亚型

Luminal A：荷尔蒙受体阳性（雌激素受体(ER)或者/以及雄激素受体(PR)阳性），HER2阴性，以及低水平的蛋白Ki-67，它可以帮助控制癌细胞生长的快慢。Luminal A癌症是低级的，趋向于生长缓慢，而且预后最好
Luminal B：荷尔蒙受体阳（性雌激素受体(ER)或者/以及雄激素受体(PR)阳性），HER2阳性或者HER2阴性伴随着较高水平的Ki-67。Luminal B这种类型的癌症一般比Luminal A型的癌症生长较快而且预后也较差。
Triple-negative/basal-like: 荷尔蒙受体阴性（性雌激素受体(ER)以及雄激素受体(PR)阴性）而且HER2阴性。这种类型的癌症在BRCA1基因突变的女性中是最常见的。研究中并不确定原因，但是这种类型的癌症也是更常见于年轻以及非洲裔美国人。
HER2-enriched: 荷尔蒙受体阴性（性雌激素受体(ER)以及雄激素受体(PR)阴性），HER2阳性。HER2-enriched 的癌症趋向于比luminal性癌症生长地快，而有较差的预后，但是它们通常可以用靶向HER2蛋白的治疗方法成功治疗。
Normal-like: 类似于Luminal A型癌症：荷尔蒙受体阳（性雌激素受体(ER)或者/以及雄激素受体(PR)阳性），HER2阴性伴随着较低水平的Ki-67。尽管normal-like型乳腺癌有一个较好的预后，它的预后依然要比Luminal A型癌症的预后要差。

R 语言实现PAM50分类：genefu

load(file = 'step1-output.Rdata')
# 每次都要检测数据
dat[1:4,1:4]  
dim(dat)
class(dat)
str(dat)
library(hgu133plus2.db)
ids=toTable(hgu133plus2SYMBOL)
head(ids)
dat=dat[ids$probe_id,]
dat[1:4,1:4] 

ids$median=apply(dat,1,median)
ids=ids[order(ids$symbol,ids$median,decreasing = T),]
ids=ids[!duplicated(ids$symbol),]
dat=dat[ids$probe_id,]
rownames(dat)=ids$symbol
dat[1:4,1:4]  
dim(dat)
dat['GAPDH',]
boxplot(dat[,1])

if(T){
  ddata=t(dat) #首先要进行转置，列名为基因，行名为样本名
  ddata[1:4,1:4]
  s=colnames(ddata);head(s)
  library(org.Hs.eg.db)
  s2g=toTable(org.Hs.egSYMBOL)
  g=s2g[match(s,s2g$symbol),1];head(g)
  #  probe Gene.symbol Gene.ID
  dannot=data.frame(probe=s,
                    "Gene.Symbol" =s, 
                    "EntrezGene.ID"=g) #将symbol转换成ENTREZID是下游分析必须的
  ddata=ddata[,!is.na(dannot$EntrezGene.ID)]#去除列的NA值
  dannot=dannot[!is.na(dannot$EntrezGene.ID),] #去除行的NA值
  head(dannot)
  library(genefu)
  # c("scmgene", "scmod1", "scmod2","pam50", "ssp2006", "ssp2003", "intClust", "AIMS","claudinLow")
  
  s<-molecular.subtyping(sbt.model = "pam50",data=ddata,
                         annot=dannot,do.mapping=TRUE) #挑选模型pam50
  table(s$subtype) #这就是一个分组依据，和我们之前根据是否是TNBC分组一样
  tmp=as.data.frame(s$subtype)
  subtypes=as.character(s$subtype)
}

取出PAM50基因，根据这些基因的表达了绘制热图，并添加分组信息，与原始分组（TNBC，noTNBC）进行对比。

library(genefu)
pam50genes=pam50$centroids.map[c(1,3)]
#我们表达数据中的这三个基因名与PAM50中的基因名不一样，所以需要先进行转换
pam50genes[pam50genes$probe=='CDCA1',1]='NUF2' 
pam50genes[pam50genes$probe=='KNTC2',1]='NDC80'
pam50genes[pam50genes$probe=='ORC6L',1]='ORC6'
x=dat
x=x[pam50genes$probe[pam50genes$probe %in% rownames(x)] ,] #判断pam50基因是否存在于我们的基因中，存在才取出对应的表达信息

table(group_list)
tmp=data.frame(group=group_list,
                  subtypes=subtypes) #制作分组信息，分组信息有两个，第一个是原始的分组信息，也就是是否是TNBC，第二个分组信息为PAM50给出的分组信息
rownames(tmp)=colnames(x)

library(pheatmap)
#原始表达数据
pheatmap(x,show_rownames = T,show_colnames = F,
         annotation_col = tmp,fontsize = 5)

#表达值进行归一化，只对基因的表达值进行样本之间的比较，不考虑基因之间的差异
x=t(scale(t(x)))
x[x>1.6]=1.6
x[x< -1.6]= -1.6
pheatmap(x,show_rownames = T,show_colnames = F,
         annotation_col = tmp,fontsize = 5)

未进行归一化之前的热图