|
GEO數(shù)據(jù)挖掘系列文-第一期-膠質(zhì)母細胞瘤 文章標題 lncRNAs PVT1 and HAR1A are prognosis biomarkers and indicate
therapy outcome for diffuse glioma patients 關(guān)鍵詞 膠質(zhì)母細胞瘤,，lncRNA 疾?。?/span>神經(jīng)膠質(zhì)瘤 1. 星形膠質(zhì)細胞（astrocytomas） · Grade I · Grade II：彌漫性星形細胞瘤 · Grade III：anaplastic astrocytoma · Grade IV：膠質(zhì)母細胞瘤（glioblastoma,，GBM） 2. 少突膠質(zhì)細胞（oligodendroglial） · Grade II · Grade III GEO數(shù)據(jù)庫編號：GSE4290 研究對象：lncRNA 實驗設(shè)計 結(jié)論：在神經(jīng)膠質(zhì)母細胞瘤中PVT1和CYTOR基因表達顯著上調(diào), HAR1A和MIAT基因表達顯著下調(diào)。 國外鏡像下載速度很慢,，所以下載之前一定要設(shè)置好國內(nèi)鏡像 bioPackages <-c(
"stringi", # 處理字符串
"GEOquery", # 下載GEO數(shù)據(jù)
"limma", # 差異分析
"ggfortify", "ggplot2", "pheatmap", "ggstatsplot", "VennDiagram", # 作圖
"clusterProfiler", "org.Hs.eg.db" # 注釋
)
## 設(shè)置軟件包的下載鏡像
local({
# CRAN的鏡像設(shè)置
r <- getOption( "repos" );
r[ "CRAN" ] <- "https://mirrors.tuna./CRAN/";
options( repos = r )
# bioconductor的鏡像設(shè)置
BioC <- getOption( "BioC_mirror" );
BioC[ "BioC_mirror" ] <- "https://mirrors.ustc.edu.cn/bioc/";
options( BioC_mirror = BioC )
})
## 安裝未安裝的軟件包
source( "https:///biocLite.R" ) #第一次使用bioconductor需要運行
lapply( bioPackages,
function( bioPackage ){
if( !require( bioPackage, character.only = T ) ){
CRANpackages <- available.packages()
if( bioPackage %in% rownames( CRANpackages) ){
install.packages( bioPackage )
}else{
BiocInstaller::biocLite( bioPackage, suppressUpdates = FALSE, ask = FALSE)
}
}
}
)
第二步 得到geneID和gene類型的對應(yīng)關(guān)系 這個關(guān)系可以從人類的GTF文件中提取出來,，GTF可以在這里下載：https://asia./info/data/ftp/index.html 下載后的文件經(jīng)過下面的shell腳本處理，即可以得到基因與基因類型的對應(yīng)關(guān)系 awk '{if(!NF || /^#/){next}}1' /public/reference/gtf/gencode/gencode.v25lift37.annotation.gtf | cut -f9 | sed 's/\"//g'| sed 's/;//g' | awk '{ print $4"\t"$8 }' |awk '{if(/^E/){next}}1'| awk '{ print $2"\t"$1 }' | sort -k 1 | uniq > gencode.v25lift37.annotation.gtf.gene2type
將處理好的文件導入R中進行后續(xù)處理,，因為這篇文章只研究lncRNA,，所以要去除編碼蛋白的基因ID {
gene2type = read.table( 'gencode.v25lift37.annotation.gtf.gene2type' )
colnames( gene2type ) = c( "gene", "type" )
}
dim( gene2type )
## [1] 58298 2
## 說明一共有五萬多個基因的對應(yīng)關(guān)系
head( gene2type )
## gene type
## 1 5S_rRNA rRNA
## 2 7SK misc_RNA
## 3 A1BG-AS1 antisense
sort( table( gene2type$type ) )
## lincRNA processed_pseudogene protein_coding
## 7601 10197 20087
## 說明絕大多數(shù)的基因是編碼蛋白的
## 剔除基因類型為“protein_coding”的對應(yīng)關(guān)系
gene2type = gene2type[ gene2type[,2] != 'protein_coding', ]
length( unique( gene2type$gene ) )
save( gene2type, file = 'Relationship_others_gene.Rdata' )
如果getGEO下載失敗，可以手動在項目主頁進行下載 library( "GEOquery" )
GSE_name = 'GSE4290'
options( 'download.file.method.GEOquery' = 'libcurl' ) #windows系統(tǒng)
gset <- getGEO( GSE_name, getGPL = T )
## getGEO函數(shù)會下載GSE項目下的所有子集,，得到的gset對象是一個list,，'GSE4290’只有一個項目，之后的實戰(zhàn)會遇到多子集的情況
## 'getGPL = T’會直接下載注釋平臺,，如果報錯,，本文最后會附上，其他進行平臺注釋的方法
save( gset, file = 'gset.Rdata' )對樣本進行不同分組,，以及探針的選取對之后的差異分析結(jié)果都會有影響。 作者研究的是GBM樣本和非腫瘤樣本在lncRNA表達上的差異,，所以先取出這180個樣本中的77個GBM樣本和23個非腫瘤樣本 options( stringsAsFactors = F )
load( './gset.Rdata' )
library( "GEOquery" )
## 取表達矩陣和樣本信息表
{
gset = gset[[1]]
exprSet = exprs( gset ) ## “GEOquery”包中的exprs函數(shù)用來取出表達矩陣
pdata = pData( gset ) ## “GEOquery”包中的pData函數(shù)用來取出樣本信息
group_list = as.character( pdata[, 35] )
}
dim( exprSet )
## [1] 54613 180
exprSet[ 1:3, 1:5 ]
## GSM97793 GSM97794 GSM97795 GSM97796 GSM97797
## 1007_s_at 10178.1 10122.9 7826.6 11098.4 8668.9
## 1053_at 388.2 517.5 352.4 609.9 430.1
## 現(xiàn)在就應(yīng)該對得到的矩陣有這樣一個印象
## 這個矩陣有180列,，列名是樣本編號，54613行,，行名是探針編號
## 矩陣的內(nèi)容就是各個樣本對應(yīng)探針檢測到的表達量
## 探針本身是沒有意義的,，所以我們下一步就要將探針轉(zhuǎn)為基因名
table( group_list )
## astrocytoma, grade 2 astrocytoma, grade 3 glioblastoma, grade 4
## 7 19 77
## non-tumor oligodendroglioma, grade 2 oligodendroglioma, grade 3
## 23 38 12
## 我們不能通過上面表達矩陣的樣本名知道這個樣本是腫瘤樣本還是非腫瘤樣本
## 而pdata記錄了各個樣本的臨床信息，就可以通過pdata得到樣本名和樣本類型的對應(yīng)關(guān)系
## 根據(jù)上面的group_list,，取出研究樣本
## 總結(jié)一下就是,，這180個病人中astrocytoma分為了二三期病人
## glioblastoma是膠質(zhì)母細胞瘤，屬于惡性細胞瘤,，只有第四期
## oligodendroglioma也分為二三期病人
## 其余樣本為對照
{
n_expr = exprSet[ , grep( "non-tumor", group_list )]
g_expr = exprSet[ , grep( "glioblastoma", group_list )]
a_expr = exprSet[ , grep( "astrocytoma", group_list )]
o_expr = exprSet[ , grep( "oligodendroglioma", group_list )]
}
## 樣本分組,，新的表達矩陣只有normal和gbm樣本
{
exprSet = cbind( n_expr, g_expr )
group_list = c(rep( 'normal', ncol( n_expr ) ),
rep( 'gbm', ncol( g_expr ) ) )
}
dim( exprSet )
exprSet[ 1:5, 1:5 ]
table( group_list )
save( exprSet, group_list, file = 'exprSet_by_group.Rdata')分組這一步就完成了，下面要去除沒有注釋的探針 并不是每一個探針都是對應(yīng)lncRNA的,，所以我們要用之前的基因和基因類型關(guān)系篩選一下,，之后再去除沒有注釋的探針 ## 篩選探針
GPL = gset@featureData@data ## 第三步getGEO函數(shù)下載數(shù)據(jù)時，直接下載了平臺,，GPL就是注釋矩陣的平臺數(shù)據(jù)
## 也就是探針和基因的對應(yīng)關(guān)系
colnames( GPL )
view( GPL )
## GPL的“ID”列是探針,，'Gene Symbol”列是基因名，將這兩列取出來
ids = GPL[ ,c( 1, 11 ) ]
## 'Gene Symbol”列格式有些特殊
## 一個探針對應(yīng)了多個基因
a<-strsplit(as.character(ids[,2]), " /// ")
tmp <- mapply( cbind, ids[,1], a )
ID2gene <- as.data.frame( tmp )
colnames( ID2gene ) = c( "id", "gene" )
load( "./Relationship_others_gene.Rdata" )
## 下面一步就是根據(jù)gene2type去除平臺中編碼蛋白的基因
ID2gene$type = gene2type[ match( ID2gene[ , 2 ], gene2type[ , 1 ] ), 2 ]
dim(ID2gene)
## [1] 59255 3
save(ID2gene, file = 'ID2gene.Rdata')
load('./ID2gene.Rdata')
load('./exprSet_by_group.Rdata')
## 這一步根據(jù)ID2gene去除沒有注釋的探針
{
exprSet = exprSet[ rownames(exprSet) %in% ID2gene[ , 1 ], ]
ID2gene = ID2gene[ match(rownames(exprSet), ID2gene[ , 1 ] ), ]
}
## 因為有些探針對應(yīng)同一個基因,，要將exprSet,，ID2gene的探針一致
## 即exprSet有的探針，ID2gene也一定有
## 即ID2gene有的探針,，exprSet也一定有
dim( exprSet )
dim( ID2gene )
tail( sort( table( ID2gene[ , 2 ] ) ), n = 12L )
## 相同基因的表達數(shù)據(jù)取最大值,，五萬多個探針,，這一步相對會運行較長時間
{
MAX = by( exprSet, ID2gene[ , 2 ],
function(x) rownames(x)[ which.max( rowMeans(x) ) ] )
MAX = as.character(MAX)
exprSet = exprSet[ rownames(exprSet) %in% MAX , ]
rownames( exprSet ) = ID2gene[ match( rownames( exprSet ), ID2gene[ , 1 ] ), 2 ]
}
exprSet = log(exprSet)
dim(exprSet)
exprSet[1:5,1:5]
save(exprSet, group_list, file = 'final_exprSet.Rdata')進行上訴操作后，基本我們就得到了之后要進行差異分析的數(shù)據(jù)集了,，我們可以先通過聚類分析看一下數(shù)據(jù)的聚類情況,，是否與分組一致 library( "ggfortify" )
## 聚類
{
colnames( exprSet ) = paste( group_list, 1:ncol( exprSet ), sep = '_' )
nodePar <- list( lab.cex = 0.3, pch = c( NA, 19 ), cex = 0.3, col = "red" )
hc = hclust( dist( t( exprSet ) ) )
png('hclust.png', res = 250, height = 1800)
plot( as.dendrogram( hc ), nodePar = nodePar, horiz = TRUE )
dev.off()
}
## 樣本過多，圖就不放了畫個PCA圖看一下,，基本是分開的,，少數(shù)樣本界限不清。 如果樣本量過少,，對于主成分分析結(jié)果與分組不一致的樣本,，可以選擇舍去 ## PCA
data = as.data.frame( t( exprSet ) )
data$group = group_list
png( 'pca_plot.png', res=80 )
autoplot( prcomp( data[ , 1:( ncol( data ) - 1 ) ] ), data = data, colour = 'group',
label =T, frame = T) + theme_bw()
dev.off()
load( "./final_exprSet.Rdata" ) library( "limma" )
{
design <- model.matrix( ~0 + factor( group_list ) )
colnames( design ) = levels( factor( group_list ) )
rownames( design ) = colnames( exprSet )
}
design
contrast.matrix <- makeContrasts( "gbm-normal", levels = design )
contrast.matrix
## design和contrast.matrix都是為了下面的差異分析制作的
load( "./ID2gene.Rdata" )
{
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.out")
}
head(nrDEG)到這里差異分析就結(jié)束了，我們畫個熱圖看一下效果 library( "pheatmap" )
{
choose_gene = head( rownames( nrDEG ), 50 )
choose_matrix = exprSet[ choose_gene, ]
choose_matrix = t( scale( t( exprSet ) ) )
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.png")
}
library( "ggplot2" )
logFC_cutoff <- with( nrDEG, mean( abs( logFC ) ) + 2 * sd( abs( logFC ) ) )
logFC_cutoff
logFC_cutoff = 1 ## 文章中設(shè)置為1
{
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()
}
作者挑出PVT1 CYTOR HAR1A MIAT這四個基因,，畫出他們在不同的腫瘤中的表達量 library( "ggstatsplot" )
load( './final_exprSet.Rdata' )
special_gene = c( 'PVT1', 'LINC00152', 'HAR1A', 'MIAT' )
## CYTOR的別名是LINC00152
for( gene in special_gene ){
filename <- paste( gene, '.png', sep = '' )
TMP = exprSet[ rownames( exprSet ) == gene, ]
data = as.data.frame(TMP)
data$group = group_list
p <- ggbetweenstats(data = data, x = group, y = TMP )
ggsave( p, filename = filename)
}
|
