在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)) <>
根据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 <>
差异结果需要先根据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]="">
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(="">
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 ="" ="" ="">
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<><0.01 &="" kk_gse$enrichmentscore="">< 0,];down_kegg$group="">
up_kegg<><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(="">
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 ="" ="" ="" ="" =" ">
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()
}
}
图略~~~~~~
联系客服
