下载乳腺癌的芯片表达数据进行差异分析
下载乳腺癌的芯片表达数据进行差异分析
在TCGA数据挖掘如此流行的现在，小编也来插一脚，我就先拿最简单的差异分析做例子吧。
比如我下载乳腺癌的所有芯片得到的表达矩阵，然后根据样本的分组，比如正常组织的表达矩阵和疾病组织的表达矩阵就可以做差异分析。其实和我们之前推出的GEO数据挖掘系列相差无几，只不过芯片表达矩阵的下载地址换成了TCGA的数据库。
假如你成功地学习了我们之前的GEO系列教程，那么本教程对你而言，最重要的知识点其实就是数据如何下载。
那么我们就开始吧！！
尤其需要注意的是，TCGA数据库里对病人来说，量化他们的基因的表达值，经历了两个阶段。起初都是用表达芯片的手段，对乳腺癌来说有接近600个表达芯片的数据，后来随着NGS技术的大行其道，乳腺癌全部的1000多个病人都进行了RNA-seq测序，考虑到有些病人既测量了他们的疾病组织也测量了他们的正常组织，所以我们其实是可以下载到1300个样本的数据。
芯片数据
1
通过firehose下载......
数据存放网址：
http://gdac.broadinstitute.org/runs/stddata__2016_01_28/data/BRCA/20160128/
下载命令：
wget http://gdac.broadinstitute.org/runs/stddata__2016_01_28/data/BRCA/20160128/gdac.broadinstitute.org_BRCA.mRNA_Preprocess_Median.Level_3.2016012800.0.0.tar.gz
2
通过UCSC Xena下载......
数据存放网址：
https://xenabrowser.net/datapages/?dataset=TCGA.BRCA.sampleMap%2FAgilentG4502A_07_3&host=https%3A%2F%2Ftcga.xenahubs.net&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443
下载命令：
wget https://tcga.xenahubs.net/download/TCGA.BRCA.sampleMap/AgilentG4502A_07_3.gz
UCSC Xena的测序数据位于：
https://xenabrowser.net/datapages/?cohort=TCGA%20Breast%20Cancer%20(BRCA)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443
所以我们在做芯片差异分析表达数据的时候，一定要下载正确的数据哦~
数据准备
rm(list = ls())
options(stringsAsFactors = F)
if(F){
array_brca=read.table('BRCA.medianexp.txt.gz',header = T,sep=' ',fill=T,quote = '')
array_brca[1:4,1:4]
array_brca=array_brca[-1,]
rownames(array_brca)=array_brca[,1]
array_brca=array_brca[,-1]
exprSet=array_brca
exprSet[1:4,1:4]
group_list=ifelse(as.numeric(substr(colnames(array_brca),14,15)) < 10,'tumor','normal')
根据TCGA样本的命名可以区分正常组织和肿瘤样本的测序结果，详情阅读最后的原文。
exprSet=as.data.frame(lapply(exprSet,as.numeric))
rownames(exprSet)=rownames(array_brca)
exprSet=na.omit(exprSet)
exprSet[1:4,1:4]
dim(exprSet)
save(exprSet,group_list,file = "tcga_brca_array_input.Rdata")
}
load(file = "tcga_brca_array_input.Rdata")
差异分析
library( "limma" )
{
design <- model.matrix( ~0 + factor( group_list ) )
colnames( design ) = levels( factor( group_list ) )
rownames( design ) = colnames( exprSet )
}
design
contrast.matrix <- makeContrasts( "tumor-normal", levels = design )
contrast.matrix
{
fit <- lmFit( exprSet, design )
fit2 <- contrasts.fit( fit, contrast.matrix )
fit2 <- eBayes( fit2 )
nrDEG = topTable( fit2, coef = 1, n = Inf )
write.table( nrDEG, file = "nrDEG_BRCA_medianexp.out")
}
head(nrDEG)
绘制热图
library( "pheatmap" )
{
tmp = nrDEG[nrDEG$P.Value < 0.05,]
差异结果需要先根据p值挑选
nrDEG_Z = tmp[ order( tmp$logFC ), ]
nrDEG_F = tmp[ order( -tmp$logFC ), ]
choose_gene = c( rownames( nrDEG_Z )[1:100], rownames( nrDEG_F )[1:100] )
choose_matrix = exprSet[ choose_gene, ]
choose_matrix = t( scale( t( choose_matrix ) ) )
choose_matrix[choose_matrix > 2] = 2
choose_matrix[choose_matrix < -2] = -2
annotation_col = data.frame( CellType = factor( group_list ) )
rownames( annotation_col ) = colnames( exprSet )
pheatmap( fontsize = 2, choose_matrix, annotation_col = annotation_col, show_rownames = F, annotation_legend = F, filename = "heatmap_BRCA_medianexp.png")
}
👁🗨瞄一眼👁🗨
■ ■ ■
绘制火山图
library( "ggplot2" )
logFC_cutoff <- with( nrDEG, mean( abs( logFC ) ) + 2 * sd( abs( logFC ) ) )
logFC_cutoff
logFC_cutoff = 1.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_array_medianexp.Rdata" )
this_tile <- paste0( 'Cutoff for logFC is ', round( logFC_cutoff, 3 ),
' The number of up gene is ', nrow(nrDEG[ nrDEG$change =='UP', ] ),
' The 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_BRCA_medianexp.png' )
dev.off()
}
👁🗨瞄一眼👁🗨
■ ■ ■
KEGG注释
library( "clusterProfiler" )
library( "org.Hs.eg.db" )
df <- bitr( rownames( nrDEG ), fromType = "SYMBOL", toType = c( "ENTREZID" ), OrgDb = org.Hs.eg.db )
head( df )
{
nrDEG$SYMBOL = rownames( nrDEG )
nrDEG = merge( nrDEG, df, by='SYMBOL' )
}
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 )
}
library( "ggplot2" )
# kegg enrich
{
{
## KEGG pathway analysis
kk.up <- enrichKEGG( gene = gene_up ,
organism = 'hsa' ,
universe = gene_all ,
pvalueCutoff = 0.99 ,
qvalueCutoff = 0.99 )
kk.down <- enrichKEGG( gene = gene_down ,
organism = 'hsa' ,
universe = gene_all ,
pvalueCutoff = 0.99 ,
qvalueCutoff = 0.99 )
}
head( kk.up )[ ,1:6 ]
head( kk.down )[ ,1:6 ]
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.png' )
}
👁🗨瞄一眼👁🗨
■ ■ ■
GSEA注释
{
### GSEA
kk_gse <- gseKEGG(geneList = geneList,
organism = 'hsa',
nPerm = 1000,
minGSSize = 30,
pvalueCutoff = 0.9,
verbose = FALSE)
head(kk_gse)[,1:6]
gseaplot(kk_gse, geneSetID = rownames(kk_gse[1,]))
down_kegg<-kk_gse[kk_gse$pvalue<0.01 & kk_gse$enrichmentScore < 0,];down_kegg$group=-1
up_kegg<-kk_gse[kk_gse$pvalue<0.01 & kk_gse$enrichmentScore > 0,];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_gsea.png')
}
👁🗨瞄一眼👁🗨
■ ■ ■
GO注释
g_list = list( gene_up = gene_up, gene_down = gene_down, gene_diff = gene_diff)
go_enrich_results <- lapply( g_list, function( gene ) {
lapply( c( 'BP', 'MF', 'CC' ) , function( ont ) {
cat( paste( 'Now process', ont ) )
ego <- enrichGO( gene = gene,
universe = gene_all,
OrgDb = org.Hs.eg.db,
ont = ont ,
pAdjustMethod = "BH",
pvalueCutoff = 0.99,
qvalueCutoff = 0.99,
readable = TRUE)
print( head( ego ) )
return( ego )
})
})
save( go_enrich_results, file = 'go_enrich_results.Rdata' )
n1 = c( 'gene_up', 'gene_down', 'gene_diff' )
n2 = c( 'BP', 'MF', 'CC' )
for ( i in 1:3 ){
for ( j in 1:3 ){
fn = paste0( 'dotplot_', n1[i], '_', n2[j], '.png' )
cat( paste0( fn, ' ' ) )
png( fn, res = 150, width = 1080 )
print( dotplot( go_enrich_results[[i]][[j]] ) )
dev.off()
}
}
图略~~~~~~
-END-
GEO数据挖掘系列文-第四期-肝细胞癌 （WGCNA）
GEO数据挖掘系列文-第五期-肝细胞癌 （多组差异分析）