一欧募、舉例介紹

本節(jié)下載GSE1009數(shù)據(jù)集粮坞，使用limma包進(jìn)行差異分析舉例晦雨。

GSE1009??

樣本量：共6個(gè)樣本，其中3個(gè)為糖尿病腎病(DN)腎小球樣本魄揉，另3個(gè)為正常腎小球樣本剪侮。

使用芯片：Affymetrix Human Genome U95 Version 2 Array。

平臺(tái)：GPL8300洛退。

二瓣俯、差異分析舉例：

###第一步：GEO數(shù)據(jù)下載

>setwd("D:\\Rfile")??

#設(shè)置工作目錄為D 盤的Rfile文件夾，大家可根據(jù)自己需要設(shè)置工作目錄兵怯。

>rm(list = ls())?

#清除所有變量

>options(stringsAsFactors=F)

#R導(dǎo)入數(shù)據(jù)時(shí)不自動(dòng)轉(zhuǎn)換為因子變量

##下面需要安裝GEOquery包彩匕，用于下載GEO數(shù)據(jù)庫的文件，已經(jīng)安裝GEOquery包不需要重復(fù)安裝媒区，直接調(diào)用GEOquery包驼仪。

#安裝命令如下(本人已經(jīng)安裝了，所以只給出安裝命令袜漩，去掉前面的注釋符號(hào)即可安裝)：

#if (!requireNamespace("BiocManager", quietly = TRUE))

#install.packages("BiocManager")

#BiocManager::install("GEOquery")

#bioconductor中的包的安裝命令經(jīng)常更新绪爸，如果大家怕安裝命令已經(jīng)更新了，現(xiàn)在的命令不能運(yùn)行宙攻，可以百度“GEOquerybioconductor”進(jìn)入主頁奠货，下拉到安裝命令處，將網(wǎng)站上的最新安裝命令復(fù)制即可座掘，這就是最新的安裝命令递惋，其他bioconductor中的包的安裝是一樣的道理柔滔。具體如下：

復(fù)制上面紅色這段即可。

##下面是用GEOquery包下載數(shù)據(jù)

>library(GEOquery)?

#加載GEOquery包

>gset <- getGEO("GSE1009",destdir = ".",AnnotGPL = F,getGPL = F) ???????????????

#下載GSE1009表達(dá)矩陣文件并賦值給gset變量,destdir = "."表示下載后的文件放在什么地方萍虽，默認(rèn)為當(dāng)前工作目錄睛廊。AnnotGPL = F表示不下載注釋文件，getGPL = F表示不下載平臺(tái)文件杉编。我們這里為了下載速度超全，就設(shè)置成不下載平臺(tái)文件和注釋文件。

#下載好后打開D盤,可看到Rfile中有一個(gè)“GSE1009_series_matrix.txt.gz”文件王财。

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

#將gset(也就是下載后的GSE1009矩陣文件)保存為R文件卵迂，文件名為“GSE1009.gset”方便以后調(diào)用。

>gset[["GSE1009_series_matrix.txt.gz"]]@annotation

#查看GSE1009_series_matrix矩陣的平臺(tái)文件

>gpl <- getGEO('GPL8300', destdir=".")

#下載這個(gè)平臺(tái)文件绒净。用于后面基因symbol注釋见咒。

總體先看看上面數(shù)據(jù)下載代碼：

代碼運(yùn)行過程及結(jié)果：

運(yùn)行代碼后Rfile中的文件：

###第二步：數(shù)據(jù)整理（將探針I(yè)D給為基因symbol）

>a=gset[[1]]

#提取gset中第一個(gè)元素（包含基因表達(dá)矩陣和注釋信息），并賦值給變量a

>dat=exprs(a)

#提取a中的表達(dá)矩陣并賦值給dat挂疆，這時(shí)的表達(dá)矩陣的基因名字為探針I(yè)D

>head(dat)

#展示表達(dá)矩陣的前6行改览，看下數(shù)據(jù)是否經(jīng)過log轉(zhuǎn)換，一般數(shù)據(jù)在20以內(nèi)缤言，是經(jīng)過log轉(zhuǎn)換的宝当，若有成百上千的數(shù)據(jù)，表示這個(gè)數(shù)據(jù)還需要進(jìn)一步log處理胆萧。

>ex <- dat

qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))

LogC <- (qx[5] > 100) ||

??(qx[6]-qx[1] > 50 && qx[2] > 0) ||

??(qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)

if (LogC) { ex[which(ex <= 0)] <- NA

dat <- log2(ex)

print("log2 transform finished")}else{print("log2 transform not needed")}

#也可以輸入上面這段代碼庆揩，如果數(shù)據(jù)已經(jīng)經(jīng)log后，顯示log2 transform not needed跌穗；如果數(shù)據(jù)沒有盡行l(wèi)og订晌，需要log程序按照代碼盡行l(wèi)og轉(zhuǎn)換，完成后顯示log2 transform finished蚌吸。

再看看轉(zhuǎn)換后的數(shù)據(jù)：

不是成千上百的啦锈拨。

>pd=pData(a)

#查看a的臨床信息，為后面選擇用于分組的變量做準(zhǔn)備

>View(pd)

#查看pd

本例中pd共有26個(gè)變量(未完全展示)羹唠，我們可以看到只有title下的字段可以分為正常奕枢、疾病，其余變量下的字段都是一樣的佩微，所以我們選擇title字段用于后續(xù)分組缝彬。

Title變量中的記錄為control 1a....Diabetes 1a.......，中間用空格分割哺眯，我們需要提取出Control和Diabetes作為分組跌造，所以需要用字符分割包來實(shí)現(xiàn)。

##安裝字符分割包stringr包,命令如下，已經(jīng)安裝過的可以直接調(diào)用壳贪，不用重復(fù)安裝。

#install.packages("stringr")

>library(stringr)

#調(diào)用stringr包

>group_list=str_split(pd$title,' ',simplify = T)[,1]

#把pd中title變量的字段寝杖，按照空格分割违施，為了給大家展示，我們先運(yùn)行分割代碼瑟幕，分割后如下圖所示磕蒲，所以我們選擇分割后的第一列即為分組狀態(tài)

>table(group_list)

#計(jì)數(shù)每一組個(gè)數(shù)

>gpl1<-Table(gpl)

>save(gpl1,file = 'gpl1.Rdata')

#我們將平臺(tái)文件轉(zhuǎn)為list格式，并賦值給gpl1只盹，將gpl1保存為R文件辣往，方便以后調(diào)用。

>colnames(Table(gpl))

#查看平臺(tái)文件的列名殖卑，我們看到有ID和gene symbol,記住gene symbol在第幾列站削，我們這里在第11列。

>View(gpl1)

#再了解一下平臺(tái)文件的數(shù)據(jù)孵稽，當(dāng)然這里大家可以直接選擇gene symbol字段也行许起，主要了解一下symbol中基因symbol值是否只有一個(gè)。

我們這里的gene symbol字段中的symbol菩鲜，有的基因就不止一個(gè)名稱园细，///后面有重名，我們需要第一個(gè)名字接校，所以需要用字符分割包再處理一下,原理同上面處理title一樣猛频。

>gpl1$`Gene Symbol`=str_split(gpl1[,11],'///',simplify = T)[,1]

現(xiàn)在是不是變成唯一啦。

>probe2symbol_df<-gpl1[,c(1,11)]

#提取平臺(tái)文件gpl1中的ID和gene symbol蛛勉，并賦值給probe2symbol_df

>colnames(probe2symbol_df)=c('probe_id','symbol')

#將列名改為probe_id和symbol

#這兩步是因?yàn)槲覒新寡埃幌朐僬{(diào)整代碼，為了方便后面代碼運(yùn)行董习，我改的烈和，大家也可以不改。

>length(unique(probe2symbol_df$symbol))

#查看symbol為唯一值的個(gè)數(shù)

>table(sort(table(probe2symbol_df$symbol)))

#查看每個(gè)基因出現(xiàn)n次的個(gè)數(shù),我們可以看到皿淋，symbol出現(xiàn)一次的基因有7050個(gè)招刹，出現(xiàn)2次有1493個(gè)。窝趣。疯暑。

>ids=probe2symbol_df[probe2symbol_df$symbol != '',]

#去掉沒有注釋symbol的探針（其實(shí)這里沒有注釋的探針數(shù)量即為上面出現(xiàn)次數(shù)最多的基因440，也就是說有440個(gè)探針沒有symbol）

>ids=probe2symbol_df[probe2symbol_df$probe_id %in% ?rownames(dat),] ??

##%in%用于判斷是否匹配哑舒，

#注意這里dat是gset表達(dá)矩陣的數(shù)據(jù)妇拯，這一步就是把平臺(tái)文件的ID和矩陣中的ID匹配。

>dat=dat[ids$probe_id,]

#取表達(dá)矩陣中可以與探針名匹配的那些，去掉無法匹配的表達(dá)數(shù)據(jù)

>ids$mean=apply(dat,1,mean) ?

#ids新建mean這一列越锈，列名為mean仗嗦，同時(shí)對(duì)dat這個(gè)矩陣按行操作，取每一行（即每個(gè)樣本在這個(gè)探針上的）的均值甘凭，將結(jié)果給到mean這一列的每一行稀拐，這里也可以改為中位值，median.

>ids=ids[order(ids$symbol,ids$mean,decreasing = T),] ???

#即先按symbol排序丹弱，相同的symbol再按照均值從大到小排列德撬，方便后續(xù)保留第一個(gè)值。

>ids=ids[!duplicated(ids$symbol),] ?

#將symbol這一列取出重復(fù)項(xiàng)躲胳，'!'為否蜓洪，即取出不重復(fù)的項(xiàng)，去除重復(fù)的gene

#取median的話坯苹，這樣就保留每個(gè)基因最大表達(dá)量結(jié)果.最后得到n個(gè)基因隆檀。

>dat=dat[ids$probe_id,]

#新的ids取出probe_id這一列，將dat按照取出的這一列北滥，用他們的每一行組成一個(gè)新的dat

>rownames(dat)=ids$symbol ?

#把ids的symbol這一列中的每一行給dat作為dat的行名

>View(dat)

現(xiàn)在就得到行名為基因名刚操，列名為樣本名字的矩陣?yán)病?/div>

>dat<-dat[-9173,]

但是我們看到矩陣最后一行，是沒有symbol名字的再芋，我們把他去掉菊霜，數(shù)字自己更改。

>save(dat,group_list,file = 'GSE1009.Rdata')

>write.csv(dat,file="GSE1009expressionmetrix_GSE.csv")

#最后我們把結(jié)果保存济赎。下一節(jié)會(huì)講解差異分析鉴逞。

現(xiàn)在看看，整體代碼如下：

運(yùn)行后的Rfile中的文件如下：

整理好的數(shù)據(jù)如下：