TCGA数据库中三阴性乳腺癌在亚洲人群中的差异表达

还是已经挖掘过很多次的TCGA数据库的BRCA数据，在哪里下载就不多说了。这期是根据种族来对数据进行分类，我们就只以亚洲人群为例，其他数据大家可以线下进行探索。

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

## 建议大家先处理表达矩阵，处理过程中会发现表达矩阵中的样本数量少于样本信息中的样本数量，
## 所以，以表达矩阵中的样本数为准
if(!file.exists('TCGA-BRCA.htseq_counts.Rdata')){
  library(data.table)
  ## fread函数读取大数据的速度较快（20s），但是读进来的数据不是data.frame的格式，需要转换
  raw_data <- fread('TCGA-BRCA.htseq_counts.tsv', sep = '\t', header = T)
  raw_data <- as.data.frame(raw_data)
  raw_data[1:4,1:4] ## 以ENTREZ ID为行名
  

  rownames(raw_data) <- raw_data[,1]
  raw_data <- raw_data[,-1]
  raw_data[1:4,1:4]
  ## 在数据的介绍页面上我们已经得知了数据的计算方法现在我们只要把它还原回去就可以了
  raw_data <- 2^raw_data-1
  raw_data=ceiling(raw_data)
  raw_data[1:4,1:4]     
    

  save(raw_data, file = 'TCGA-BRCA.htseq_counts.Rdata')
}else{
  load('TCGA-BRCA.htseq_counts.Rdata')
}

## 下载样本信息
phenotype_file <- read.table('TCGA-BRCA.GDC_phenotype.tsv.gz',header = T,sep = '\t',quote = '')
## 只保留有表达矩阵的样本信息，1217个样本
phenotype_file <- phenotype_file[phenotype_file[,1] %in% colnames(raw_data),] 

## 取出含有'race.demographic'信息的列，包含了种族信息
colnames_num_race <- grep('race.demographic',colnames(phenotype_file))
rph <- phenotype_file[,colnames_num_race]
table(rph)
## asian     black or african american       not reported        white 
##  62                            189                 95          869

## 根据种族将样本分组
col_asian <- grep("asian", rph) # 62个亚洲人
col_black <- grep("black or african american", rph) # 189个黑人
col_white <- grep("white", rph) # 869个白人

asian_sample <- phenotype_file[col_asian, 1]
black_sample <- phenotype_file[col_black, 1]
white_sample <- phenotype_file[col_white, 1]
save(asian_sample, black_sample, white_sample, file = 'race_sample.Rdata')

## 取出含有'receptor_status'信息的列，按照tnbc分类
colnames_num_tnbc <- grep('receptor_status',colnames(phenotype_file))
colnames(phenotype_file)[colnames_num_tnbc]
eph <- phenotype_file[,colnames_num_tnbc[1:3]]
## 根据ER,PR,HER2将样本分组
tnbc_row_num <- apply(eph, 1, function(x) sum(x =='Negative')) ## 均为阴性的为三阴性乳腺癌
colnames_num_tnbc <- phenotype_file[tnbc_row_num == 3, 1] # 111个三阴性乳腺癌样本

## 将表达矩阵拆分为tnbc样本和非tnbc样本分组
tnbc_exprSet=raw_data[,colnames_num_tnbc]
non_tnbcE_exprSet=raw_data[, setdiff(colnames(raw_data), colnames_num_tnbc)]

## 计算每个种族中的tnbc样本
asian_tnbc_sample <- intersect(colnames_num_tnbc, asian_sample) # 8个
black_tnbc_sample <- intersect(colnames_num_tnbc, black_sample) # 32个
white_tnbc_sample <- intersect(colnames_num_tnbc, white_sample) # 64个

save(asian_tnbc_sample, black_tnbc_sample, white_tnbc_sample, file = 'tnbc_sample.Rdata')

save(tnbc_exprSet, non_tnbcE_exprSet, file = 'group_exprset.Rdata')

到这里样本就已经处理好了，接下来用LIMMA包做差异分析。

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

load('race_sample.Rdata')
load('tnbc_sample.Rdata')
load('group_exprset.Rdata')

## source("http://bioconductor.org/biocLite.R")
## biocLite("edgeR")
library("edgeR")
library('ggplot2')
library('limma')

## 挑出亚洲人群中的tnbc和非tnbc样本，这样的方法可以确定矩阵和group_list是一一对应的
tnbc_asian = tnbc_exprSet[,asian_tnbc_sample]
non_tnbcE_asian = non_tnbcE_exprSet[, setdiff(asian_sample,asian_tnbc_sample)]
exprSet = cbind(tnbc_asian, non_tnbcE_asian)
group_list = c(rep('tnbc',ncol(tnbc_asian)),rep('nontnbc',ncol(non_tnbcE_asian)))

## 接下来就是寻常的limma差异分析
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
dge <- DGEList(counts=exprSet)
dge <- calcNormFactors(dge)
logCPM <- cpm(dge, log=TRUE, prior.count=3)
v <- voom(dge,design,plot=TRUE, normalize="quantile")
fit <- lmFit(v, design)
cont.matrix=makeContrasts(contrasts=c('tnbc-nontnbc'),levels = design)
fit2=contrasts.fit(fit,cont.matrix)
fit2=eBayes(fit2)
tempOutput = topTable(fit2, coef='tnbc-nontnbc', n=Inf)
DEG_limma_voom = na.omit(tempOutput)
head(DEG_limma_voom)
save(exprSet, group_list, DEG_limma_voom, file = "DEG.Rdata")

到这里差异结果也有了，就可以根据挑选出的差异基因做火山图

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

load(file = 'DEG.Rdate')

## 画火山图的一般流程
nrDEG <- DEG_limma_voom
library( "ggplot2" )
logFC_cutoff = 2
{
  nrDEG$change = as.factor( ifelse( nrDEG$P.Value < 0.05 & abs(nrDEG$logFC) > logFC_cutoff,
                                    ifelse( nrDEG$logFC > logFC_cutoff , 'UP', 'DOWN' ), 'NOT' ) )

  save( nrDEG, file = "nrDEG.Rdata" )

  this_tile <- paste0( 'Cutoff for logFC is ', round( logFC_cutoff, 3 ),
                       '\nThe number of up gene is ', nrow(nrDEG[ nrDEG$change =='UP', ] ),
                       '\nThe number of down gene is ', nrow(nrDEG[ nrDEG$change =='DOWN', ] ) )

  volcano = ggplot(data = nrDEG, aes( x = logFC, y = -log10(P.Value), color = change)) +
    geom_point( alpha = 0.4, size = 1.75) +
    theme_set( theme_set( theme_bw( base_size = 15 ) ) ) +
    xlab( "log2 fold change" ) + ylab( "-log10 p-value" ) +
    ggtitle( this_tile ) + theme( plot.title = element_text( size = 15, hjust = 0.5)) +
    scale_colour_manual( values = c('blue','black','red') )
  print( volcano )
  ggsave( volcano, filename = 'volcano.png' )
  dev.off()
}

下面我们就可以做通路注释了

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

load(file = 'DEG.Rdate')

## 注释
library( "clusterProfiler" )
library( "org.Hs.eg.db" )

nrDEG = DEG_limma_voom
logFC_cutoff = 2
nrDEG$change = as.factor( ifelse( nrDEG$P.Value < 0.05 & abs(nrDEG$logFC) > logFC_cutoff,
                                  ifelse( nrDEG$logFC > logFC_cutoff , 'UP', 'DOWN' ), 'NOT' ) )
table(nrDEG$change)
keytypes(org.Hs.eg.db)
library("stringr")
rownames( nrDEG ) <- str_sub(rownames( nrDEG ), start = 1, end = 15)
nrDEG$ENSEMBL <- rownames( nrDEG )
df <- bitr( rownames( nrDEG ), fromType = "ENSEMBL", toType = c( "ENTREZID" ), OrgDb = org.Hs.eg.db )
head( df )
{
  nrDEG$SYMBOL = rownames( nrDEG )
  nrDEG = merge( nrDEG, df, by='ENSEMBL' )
}
head( nrDEG )


{
  gene_up = nrDEG[ nrDEG$change == 'UP', 'ENTREZID' ] 
  gene_down = nrDEG[ nrDEG$change == 'DOWN', 'ENTREZID' ]
  gene_diff = c( gene_up, gene_down )
  gene_all = as.character(nrDEG[ ,'ENTREZID'] )
}

{
  geneList = nrDEG$logFC
  names( geneList ) = nrDEG$ENTREZID
  geneList = sort( geneList, decreasing = T )
}

{
  ## KEGG pathway analysis
  kk.up <- enrichKEGG(   gene          =  gene_up    ,
                         organism      =  'hsa'      ,
                         universe      =  gene_all   ,
                         pvalueCutoff  =  0.8        ,
                         qvalueCutoff  =  0.8        )
  kk.down <- enrichKEGG( gene          =  gene_down  ,
                         organism      =  'hsa'      ,
                         universe      =  gene_all   ,
                         pvalueCutoff  =  0.8        ,
                         qvalueCutoff  =  0.8       )
}

head( kk.up )[ ,1:6 ]
head( kk.down )[ ,1:6 ]

library( "ggplot2" )
{
  kegg_down_dt <- as.data.frame( kk.down )
  kegg_up_dt <- as.data.frame( kk.up )
  down_kegg <- kegg_down_dt[ kegg_down_dt$pvalue < 0.05, ]
  down_kegg$group = -1
  up_kegg <- kegg_up_dt[ kegg_up_dt$pvalue < 0.05, ]
  up_kegg$group = 1

  dat = rbind( up_kegg, down_kegg )
  dat$pvalue = -log10( dat$pvalue )
  dat$pvalue = dat$pvalue * dat$group

  dat = dat[ order( dat$pvalue, decreasing = F ), ]

  g_kegg <- ggplot( dat, 
    aes(x = reorder( Description, order( pvalue, decreasing=F ) ), y = pvalue, fill = group)) + 
    geom_bar( stat = "identity" ) + 
    scale_fill_gradient( low = "blue", high = "red", guide = FALSE ) + 
    scale_x_discrete( name = "Pathway names" ) +
    scale_y_continuous( name = "log10P-value" ) +
    coord_flip() + theme_bw() + 
    theme( plot.title = element_text( hjust = 0.5 ) ) +
    ggtitle( "Pathway Enrichment" ) 
  print( g_kegg )
  ggsave( g_kegg, filename = 'kegg_up_down_edgeR.png' )
}

好了，到这里就结束了，大家可以也做一下其他两个种族的差异比较