我的GEO練習(xí)

GEO數(shù)據(jù)挖掘練習(xí)

搜索文獻(xiàn)，找到GSE號(hào)

? ? ? 成骨細(xì)胞礦化對(duì)基因的影響（Matrix mineralization controls gene expression in osteoblasts）

GSE114237

實(shí)驗(yàn)設(shè)計(jì)：4個(gè)非礦化樣本 4個(gè)礦化樣本

R代碼參考

https://mp.weixin.qq.com/s/Z4fK6RObUEfjEyY_2VS4Nw

安裝或載入所需各種包

library(stringi)

library(GEOquery)

library(limma)

library(ggfortify)

library(ggstatsplot)

library(VennDiagram)

? ? 如果有沒安裝的包? 先設(shè)置鏡像

r[ "CRAN" ] <- "- "https://mirrors.tuna.tsinghua.edu.cn/CRAN/";

/";

options( repos = r )

BioC <- getOption( "BioC_mirror" );

BioC[ "BioC_mirror" ] <- "- "https://mirrors.ustc.edu.cn/bioc/";

/";

options( BioC_mirror = BioC )

下載GEO數(shù)據(jù)

library( "GEOquery" )

GSE_name = 'GSE114237'? #這填要下的號(hào)

options( 'download.file.method.GEOquery' = 'libcurl' )

gset <- getGEO( GSE_name, getGPL = T )? #getGPL = T意思就是下載這個(gè)GSE用的平臺(tái)

save( gset, file = 'gset.Rdata' )

制作數(shù)據(jù)集

load( './'./gset.Rdata' )

View（gset）

library( "GEOquery" )

gset = gset[[1]]? #這時(shí)候gset是個(gè)list 提取list中的第一個(gè)元素 就是表達(dá)相關(guān)信息

exprSet = exprs( gset )? #Geoquery包函數(shù)exprs用來提取表達(dá)矩陣

pdata = pData( gset )? #Geoquery包函數(shù)pData用來提取樣本信息

group_list = as.character( pdata[, 24] )? #grouplist是個(gè)分組信息 作圖用的 該取哪列取哪列 進(jìn)去看一眼改個(gè)數(shù) 折騰半天? 他給分成min組合demin組（礦化非礦化）

dim( exprSet )

exprSet[ 1:3, 1:5 ]? #看一眼

n_expr = exprSet[ , grep( "^min", group_list )]? #^這個(gè)符號(hào)表示取開頭是什么什么的東西

View(n_expr)

g_expr = exprSet[ , grep( "^demin", group_list )]? #n和g這倆命名照代碼上抄的 應(yīng)該自己改一個(gè)能看懂的

exprSet = cbind( n_expr, g_expr )

View(exprSet)

group_list = c(rep( 'min', ncol( n_expr ) ),

rep( 'demin', ncol( g_expr ) ) )

dim( exprSet )? #就是讓min和demin重復(fù)多少列就多少遍完寫grouplist里頭

exprSet[ 1:8, 1:8 ]? #瞅一眼

table( group_list )

save( exprSet, group_list, file = 'exprSet_by_group.Rdata')? ? #保存Rdata的寫法

篩選探針? 去除沒有注釋的 或者沒測(cè)出來東西的探針

GPL = gset@featureData@data? #之前下GSE同時(shí)下了平臺(tái) GPL就是平臺(tái)對(duì)探針的注釋（意思就是每個(gè)探針都對(duì)應(yīng)啥玩意）

colnames( GPL )

view( GPL )

看完感覺不對(duì)? 因?yàn)镚PL里頭探針號(hào)沒有對(duì)應(yīng)的Gene symbol? 得下個(gè)包找找這個(gè)平臺(tái)的對(duì)應(yīng)信息 先去GEO上找這個(gè)平臺(tái)對(duì)應(yīng)什么GPL號(hào)? 然后根據(jù)GPL號(hào)

? 上曾老師博客里找這個(gè)平臺(tái)對(duì)應(yīng)什么探針? 找著了就把包下下來

BiocInstaller::biocLite('hugene10sttranscriptcluster.db')? #這包里有信息

library(hugene10sttranscriptcluster.db)

toTable(hugene10sttranscriptclusterSYMBOL)? #toTable可以查看包里的對(duì)應(yīng)信息

hgnc_id=toTable(hugene10sttranscriptclusterSYMBOL)? #把這個(gè)信息賦值

View(hgnc_id)

合并提取出來的矩陣和對(duì)應(yīng)信息

exprSet = exprSet[ rownames(exprSet) %in% hgnc_id[ , 1 ], ]

hgnc_id = hgnc_id[ match(rownames(exprSet), hgnc_id[ , 1 ] ), ]

dim( exprSet )

dim( hgnc_id )

tail( sort( table( hgnc_id[ , 2 ] ) ), n = 12L )? #n=12L啥意思沒看懂? 回頭查查說明書 先往下進(jìn)行

取出現(xiàn)頻率最大的交集

MAX = by( exprSet, hgnc_id[ , 2 ],

? ? ? ? ? function(x) rownames(x)[ which.max( rowMeans(x) ) ] )

MAX = as.character(MAX)

exprSet = exprSet[ rownames(exprSet) %in% MAX , ]

rownames( exprSet ) = hgnc_id[ match( rownames( exprSet ), hgnc_id[ , 1 ] ), 2 ]

exprSet = log(exprSet)

dim(exprSet)

exprSet[1:5,1:5]

save(exprSet, group_list, file = 'final_exprSet.Rdata')

聚類分析

install.packages("dplyr")

library(ggfortify)

library(stringi)

install.packages("stringi")

library(ggfortify)

colnames( exprSet ) = paste( group_list, 1:ncol( exprSet ), sep = '_' )? #sep = '_'分隔符默認(rèn)

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圖

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


出完圖發(fā)現(xiàn)有兩個(gè)樣本數(shù)據(jù)偏差跟別的比特別大? 決定給他倆刪了

剔除不好的數(shù)據(jù)樣本

exprSet= exprSet[,1:2&4:7]

exprSet = exprSet[,-3]

exprSet = exprSet[,-7]

完事又發(fā)現(xiàn)group list又不對(duì)了? 好像還得改一下group list 可能有不太蠢的辦法 但我用了下面這個(gè)

#重新定義一下group list

table(group_list)

n_expr=n_expr[,-3]

g_expr=g_expr[,-4]

group_list = c(rep( 'min', ncol( n_expr ) ),

? ? ? ? ? ? ? rep( 'demin',? ? ncol( g_expr ) ) )

table(group_list)

再做一次PCA看看? 結(jié)果還湊合吧


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

差異分析

library( "limma" )

design <- model.matrix( ~0 + factor( group_list ) )

colnames( design ) = levels( factor( group_list ) )

rownames( design ) = colnames( exprSet )

design

contrast.matrix <- makeContrasts( "demin-min", levels = design )? #這塊得改引號(hào)里的名 該改啥改啥

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.out")

head(nrDEG)

熱圖

library( "pheatmap" )

choose_gene = head( rownames( nrDEG ), 50 )? #取了前50個(gè)差異最大的基因? 這數(shù)可以自己定

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 = 5, choose_matrix, annotation_col = annotation_col, show_rownames = T,

? ? ? ? ? annotation_legend = T, filename = "heatmap.png")


火山圖

library( "ggplot2" )

logFC_cutoff <- with( nrDEG, mean( abs( logFC ) ) + 2 * sd( abs( logFC ) ) )

logFC_cutoff

logFC_cutoff = 1? #這1就是火山圖橫坐標(biāo) 把1里頭的都算成沒啥變異? 這個(gè)也可以自己定 根據(jù)圖定吧 可以調(diào)調(diào)

nrDEG$change = as.factor( ifelse( nrDEG$P.Value < 0.05 & abs(nrDEG$logFC) > logFC_cutoff,

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ifelse( nrDEG$logFC > logFC_cutoff , 'UP', 'DOWN' ), 'STABLE' ) )? #這down up stable都是自己寫的? 可以寫別的 我感覺寫這個(gè)比較好

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) +? #這里頭這個(gè)size是火山圖上那個(gè)小點(diǎn)點(diǎn)的大小 可以改 別的也能改 這里頭帶數(shù)的都可以改改包括下面的? 我還沒試

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

save.image(file = 'project1.Rdata') #保存了一下? 有點(diǎn)干太多了


葫蘆圖

隨便挑了4個(gè)基因? 做了一下

special_gene = c( 'GRB14', 'KDM7A', 'DPP4', 'MYPN' )

for( gene in special_gene ){

# gene='GRB14' # 先做一個(gè)看看好使不? 要好使就注釋了

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,ylab = gene,xlab = 'g')

? ggsave( p, filename = filename)

}

save.image(file = 'project1.Rdata')


最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

  • 十三屆全國人大一次會(huì)議今日將舉行第五次全體會(huì)議，選舉新一屆國家機(jī)構(gòu)領(lǐng)導(dǎo)人，表決關(guān)于國務(wù)院機(jī)構(gòu)改革方案的決定草案等。...
    羊羊sasa閱讀 129評(píng)論 0贊 0
  • 當(dāng)現(xiàn)實(shí)不盡如人意的時(shí)候，我們總是會(huì)幻想，如果…… 今天麻婆要推薦的影片是—— 白日夢(mèng)想家 The Secret L...
    麻婆電影閱讀 361評(píng)論 0贊 1
  • 曾經(jīng)我有一只叫做大煩的貓，其實(shí)到現(xiàn)在我也不愿意說“曾經(jīng)”這兩個(gè)字，因?yàn)槲乙恢庇X得有一天，它還會(huì)回來的。 大煩是...
    小渦魚閱讀 407評(píng)論 3贊 2

友情鏈接更多精彩內(nèi)容