一、數(shù)據(jù)集介紹
GEO鏈接： https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62452
芯片平臺(tái)： GPL6244 [HuGene-1_0-st] Affymetrix Human Gene 1.0 ST Array
平臺(tái)鏈接： https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL6244
樣品信息： 61個(gè)正常樣本與69個(gè)胰腺導(dǎo)管腺癌（PDAC）樣本
文章鏈接：A Novel MIF Signaling Pathway Drives the Malignant Character of Pancreatic Cancer by
Targeting NR3C2. Cancer Res 2016 Jul 1;76(13)
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2755578/
二藐守、核心步驟
1.運(yùn)行準(zhǔn)備與獲取輸入數(shù)據(jù)及輸入數(shù)據(jù)檢查
- 01步R包加載、加載數(shù)據(jù)處理腳本與加載繪圖腳本
清空環(huán)境胜宇，加載R包與腳本----
source("scripts/step0-load-packages.R") #加載R包
source("scripts/functions.R") #加載數(shù)據(jù)處理腳本
source('scripts/draw-figures.R')#加載繪圖腳本
- 02步獲取輸入數(shù)據(jù)
- 主要獲得三方面輸入數(shù)據(jù)：（表達(dá)矩陣乏沸、臨床信息與芯片注釋包）
# 獲取表達(dá)矩陣袜茧、臨床信息與芯片注釋包----
getOption('timeout')
options(timeout=10000)
gse_number = "GSEGSE62452"#需要指定gse_number
f="step1-data.Rdata"
source("scripts/step1-download-data.R") #下載三方面輸入數(shù)據(jù)
- 03步輸入數(shù)據(jù)檢查（兩方面）
- 主要檢查兩方面：1.檢查exp列名與pd的行名順序是否完全一致 2.檢查是否需要取log（箱線圖）
# 數(shù)據(jù)兩方面檢查----
## 01-檢查1：exp列名與pd的行名順序是否完全一致
p = identical(rownames(pd),colnames(probeM));p
if(!p) probeM = probeM[,match(rownames(pd),colnames(probeM))]
## 02-檢查2：是否需要取log
probeM[1:4,1:4]#檢測(cè)整體探針是否需要取log，小于20則不需要取
#probeM =log2(probeM+1) #不取log則不要運(yùn)行此步倍踪，不運(yùn)行可以加個(gè)注釋
## 03-開(kāi)始繪圖
draw_p1_boxplot(probeM)#箱線圖初看組內(nèi)與組間測(cè)序差異 p1
2.整理數(shù)據(jù)
- 整理兩方面：1.整理獲得去除冗余探針的表達(dá)矩陣 2.整理獲得分組信息
- 04步整理探針矩陣獲取基因表達(dá)矩陣系宫，去除冗余探針
# 整理探針矩陣獲取基因表達(dá)矩陣，去除冗余探針----
## 01將探針表達(dá)矩陣轉(zhuǎn)換為基因表達(dá)矩陣惭适，探針I(yè)D與基因symbol一一對(duì)應(yīng)
geneM = probeM2geneM(ids,probeM)
## 02檢查轉(zhuǎn)換情況笙瑟，看兩個(gè)管家基因表達(dá)范圍（正常10-15之間）
fivenum(geneM['ACTB',])
fivenum(geneM['GAPDH',])
- 05步整理獲得分組信息(根據(jù)生物學(xué)背景及研究目的人為分組)
# 獲取樣品分組信息并進(jìn)行分組----
## 01查找分組信息所在列，通過(guò)使用ifelse與str_detect進(jìn)行分組
pd$source_name_ch1 #找列癞志，看有無(wú)分組信息
Group=ifelse(str_detect(pd$source_name_ch1,"Non"),"Normal","PDAC") #進(jìn)行分組
## 02將Group轉(zhuǎn)換成因子往枷，指定levels；對(duì)照組在前，處理組在后
Group = factor(Group,levels = c("Normal","PDAC"))
table(Group)
## 03保存基因表達(dá)矩陣與樣品分組信息
save(geneM,Group,file = "step2-genemGro.Rdata")
- 3個(gè)質(zhì)控圖及其拼圖(樣本PCA圖错洁、高變基因熱圖與樣品相關(guān)性圖)
- 3個(gè)圖的作用：從樣本間的差異與樣本內(nèi)的差異來(lái)看測(cè)序效果與分組是否有差異
# 繪制三個(gè)質(zhì)控圖----
## 01PCA看組內(nèi)與組間整體差異 p2
p2=draw_p2_PCA(probeM,gse_number,Group)
## 02高變基因熱圖（畫(huà)sd排名前1000基因的熱圖）p3
p3=draw_p3_pheatmap(geneM,Group,gse_number)
## 03樣品相關(guān)性熱圖 p4
p4=draw_p4_colsample(Group,exp,gse_number)
## 04拼圖與保存
##注熱圖得轉(zhuǎn)換為ggplot格式才能進(jìn)行拼圖
p2+as.ggplot(p3)+ as.ggplot(p4)
ggsave("png/p234.pdf",width = 12,height = 5)
p2 |(as.ggplot(p3)/ as.ggplot(p4) )
- 4.差異基因-limma分析（火山圖與熱圖秉宿，兩圖）
# 差異分析與繪制火山圖與熱圖----
## 01差異分析-limma分析（得到6列差異基因表達(dá)矩陣）
deg=DEG_limma(Group,probeM)
## 02第7列：加上下調(diào)基因列（為繪制火山圖與差異基因熱圖做準(zhǔn)備）
logFC_t=1;P.Value_t = 0.05 #設(shè)置顯著差異過(guò)濾條件
k1 = (deg$P.Value < P.Value_t)&(deg$logFC < -logFC_t)
k2 = (deg$P.Value < P.Value_t)&(deg$logFC > logFC_t)
deg <- mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable")))
## 03第8列：加ENTREZID列及其余列,為富集分析
s2e <- bitr(deg$symbol,
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)#人類(lèi)
#其他物種http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
## 04去重ENTREZEID對(duì)應(yīng)多個(gè)基因的現(xiàn)象
deg=deg[order(deg$symbol,abs(deg$logFC),decreasing = T),] #先按絕對(duì)值排
deg = deg[!duplicated(deg$symbol),]
table(deg)
## 05差異基因好了，繪制火山圖 p5
p5=draw_p5_volcano(Group,deg,gse_number)
## 06前十與后十個(gè)差異基因繪制熱圖 p6
### 001挑選前十與后十變化最為顯著的差異基因
dat1=filter(deg,change!="stable") #去組別間不發(fā)生顯著變化的基因
dat2=arrange(dat1,logFC)#排序
cg = c(head(dat2$symbol,10),tail(dat2$symbol,10)) #取前十與后十
n=geneM[cg,] ##差異表達(dá)矩陣?yán)锶∏笆c后十
dim(n)
## 002 20個(gè)差異表達(dá)矩陣好了屯碴，繪制熱圖
p6 <- draw_p6_geneheatmap(Group,n)
p6
## 07拼圖與保存圖片
plot_grid(p5, as.ggplot(p6))
ggsave("png/p56.pdf",width=15, height=15)
plot_grid(p5, as.ggplot(p6))
- 5利用GO與KEGG富集分析差異功能與差異通路
# GO富集分析----
## 01先對(duì)差異基因進(jìn)行富集描睦，包括上調(diào)、下調(diào)與整體差異基因
ego=GO_enrich(deg)
## 02對(duì)富集后的結(jié)果可視化
p7=draw_p7_GObarplot(ego)
# KEGG富集分析---
## 01先對(duì)上下調(diào)差異基因進(jìn)行KEGG富集
load(paste0(gse_number,"_GO.Rdata")) #輸入已有的上下調(diào)參數(shù)
KEGG_enrich(deg)
## 02按照q值分別挑選10條上調(diào)與下調(diào)最顯著的KEGG通路
load(paste0(gse_number,"_KEGG.Rdata"))
top10updownKEGG(kk.down,kk.up)
## 03繪制top10down與top10up的KEGG條形圖
load("top10updown.Rdata")
p8 <- draw_p8_KEGGBarplot (up.data,down.data)
## 04拼圖（GO與KEGG結(jié)果拼圖）
plot_grid(p7,p8)
ggsave("png/p78.pdf",width = 15,height=15)
plot_grid(p7,p8)
- 6 GSEA富集分析
# KEGG的GSEA富集分析---
## 01KEGG的GSEA富集
KeggGSEA_enrich(deg)
## 02挑選前6個(gè)上調(diào)與后6個(gè)下調(diào)通路
load('Rdata/gsea_kk.Rdata')
dat=top6downupGSEA(deg)
## 03繪制KEGG的GSEA富集分析條形圖
p9 <- draw_p9_gsea(dat)
7.KEGG的GSEA富集分析具體通路可視化（折線圖）
# 針對(duì)上面獲得的12條通路進(jìn)行具體GSEA可視化
## 01先看下調(diào)的6個(gè)通路并繪折線圖
load("top6downup.Rdata")
p10 <- gseaplot2(kk, geneSetID = rownames(down_k))
p10
ggsave("png/p10.png")
## 02再看上調(diào)的6個(gè)通路并繪折線圖
p11=gseaplot2(kk, geneSetID = rownames(up_k))
ggsave("png/p11.png")
## 03拼圖并保存
plot_grid(p9,p10,p11)
ggsave("png/p91011.pdf",width = 15,height=15)
plot_grid(p9,p10,p11)
以上腳本均來(lái)自生信技能樹(shù)导而，生信技能樹(shù)目前正在進(jìn)行萬(wàn)能芯片數(shù)據(jù)的處理忱叭，Jimmy老師現(xiàn)場(chǎng)演示。詳情請(qǐng)見(jiàn)生信技能樹(shù)公眾號(hào)今艺，以上腳本以及各種芯片數(shù)據(jù)的處理韵丑，你也可以手到擒來(lái)，輕松復(fù)現(xiàn)虚缎。
