在TCGA數(shù)據(jù)挖掘如此流行的現(xiàn)在，小編也來(lái)插一腳，我就先拿最簡(jiǎn)單的差異分析做例子吧。

比如我下載乳腺癌的所有芯片得到的表達(dá)矩陣，然后根據(jù)樣本的分組，比如正常組織的表達(dá)矩陣和疾病組織的表達(dá)矩陣就可以做差異分析。其實(shí)和我們之前推出的GEO數(shù)據(jù)挖掘系列相差無(wú)幾，只不過(guò)芯片表達(dá)矩陣的下載地址換成了TCGA的數(shù)據(jù)庫(kù)。

假如你成功地學(xué)習(xí)了我們之前的GEO系列教程，那么本教程對(duì)你而言，最重要的知識(shí)點(diǎn)其實(shí)就是數(shù)據(jù)如何下載。

那么我們就開始吧！！

下載

地址

尤其需要注意的是，TCGA數(shù)據(jù)庫(kù)里對(duì)病人來(lái)說(shuō)，量化他們的基因的表達(dá)值，經(jīng)歷了兩個(gè)階段。起初都是用表達(dá)芯片的手段，對(duì)乳腺癌來(lái)說(shuō)有接近600個(gè)表達(dá)芯片的數(shù)據(jù)，后來(lái)隨著NGS技術(shù)的大行其道，乳腺癌全部的1000多個(gè)病人都進(jìn)行了RNA-seq測(cè)序，考慮到有些病人既測(cè)量了他們的疾病組織也測(cè)量了他們的正常組織，所以我們其實(shí)是可以下載到1300個(gè)樣本的數(shù)據(jù)。

芯片數(shù)據(jù)

通過(guò)firehose下載......

數(shù)據(jù)存放網(wǎng)址：

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

通過(guò)UCSC Xena下載......

數(shù)據(jù)存放網(wǎng)址：

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的測(cè)序數(shù)據(jù)位于：

https://xenabrowser.net/datapages/?cohort=TCGA%20Breast%20Cancer%20(BRCA)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443

所以我們?cè)谧鲂酒町惙治霰磉_(dá)數(shù)據(jù)的時(shí)候，一定要下載正確的數(shù)據(jù)哦~

數(shù)據(jù)準(zhǔn)備

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)) <>

根據(jù)TCGA樣本的命名可以區(qū)分正常組織和腫瘤樣本的測(cè)序結(jié)果，詳情閱讀最后的原文。

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 <>

差異結(jié)果需要先根據(jù)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()

}

圖略~~~~~~

本站僅提供存儲(chǔ)服務(wù)，所有內(nèi)容均由用戶發(fā)布，如發(fā)現(xiàn)有害或侵權(quán)內(nèi)容，請(qǐng)點(diǎn)擊舉報(bào)。

免费视频淫片aa毛片_日韩高清在线亚洲专区vr_日韩大片免费观看视频播放_亚洲欧美国产精品完整版