由于太久沒有分析過bulk數(shù)據(jù)了麻裳，前幾天突然被要求分析一個文章的數(shù)據(jù)時有點(diǎn)亂了手腳狞换，所以在這里想總結(jié)一下關(guān)于芯片數(shù)據(jù)與seq數(shù)據(jù)的差異分析流程谭网。
假如你是小白，也剛剛了解NGS的冰山一角（我指的是轉(zhuǎn)錄組）坝撑，就不要隨隨便便拿到數(shù)據(jù)去跑你那所謂的流程代碼，因?yàn)楹苡锌赡苁清e的粮揉，還得意洋洋（我會鄙視你的）
所以這篇blog想要闡明的問題就是：芯片數(shù)據(jù)和seq數(shù)據(jù)不同的差異分析流程
首先聲明绍载，你得去了解什么是基因芯片測序和seq測序，這里我不講滔蝉，自行Google击儡！
一個必要的好習(xí)慣：當(dāng)別人給你一個bulk的轉(zhuǎn)錄組數(shù)據(jù)時，你首先需要做的事蝠引，問這個測序是什么方式測序（芯片or二代測序）阳谍，因?yàn)椴粏柲銦o法分析或者分析的也是大錯特錯！ｓΩ拧矫夯！
ok，咱們就直接開始
1 芯片數(shù)據(jù)
假如你拿到是芯片數(shù)據(jù)那么就用下面的分析流程吧
首先你得需要一個表達(dá)矩陣吧吊洼，若你想分析文章的數(shù)據(jù)是芯片數(shù)據(jù)训貌，你需要怎么下載呢，很簡單用一個包
1.1下載和處理數(shù)據(jù)
#數(shù)據(jù)下載
rm(list = ls())
library(GEOquery)
gse_number = "GSE56649" #這里你只需要改成你要下載文獻(xiàn)的GEO號
eSet <- getGEO(gse_number,
destdir = '.',
getGPL = F)。 #這個包會自動聯(lián)網(wǎng)去GEO數(shù)據(jù)庫下載該文章的數(shù)據(jù)
class(eSet)
length(eSet)
eSet = eSet[[1]]
#(1)提取表達(dá)矩陣exp
exp <- exprs(eSet) #表達(dá)矩陣在assay data中exprs
exp[1:4,1:4]
if(T){exp = log2(exp+1)}#避免有0值的存在 加1不影響大械莼Α（觀察表達(dá)矩陣中有無取過log）
boxplot(exp) #觀察是否差不多的箱長豺鼻，如果不一樣
if(T){exp = limma::normalizeBetweenArrays(exp)}
boxplot(exp) #這時候箱長就會一樣的
#(2)提取臨床信息
pd <- pData(eSet)
#(3)調(diào)整pd的行名順序與exp列名完全一致
p = identical(rownames(pd),colnames(exp));p
if(!p) exp = exp[,match(rownames(pd),colnames(exp))]
#(4)提取芯片平臺編號
gpl_number <- eSet@annotation
對上面代碼塊中關(guān)于表達(dá)矩陣的處理進(jìn)行解釋，
通過exp[1:4,1:4]判斷數(shù)據(jù)長啥樣款慨，判斷表達(dá)矩陣是否已經(jīng)取過log儒飒，假如表達(dá)矩陣的數(shù)值范圍在0-15，最大不超過20檩奠，那么就是取過log了桩了，不要運(yùn)行exp = log2(exp+1)，反之就要運(yùn)行埠戳，這里解釋一下問什么要取log井誉，假如沒有取log的話，差異logFC會大的離譜整胃，但是你取了多次log的話送悔，后面的logFC會很小，我們對表達(dá)值取log不僅是對logFC有作用爪模，而且對后面畫熱圖數(shù)值的歸一化也很有用 然后我們就要?dú)w一化處理limma::normalizeBetweenArrays(exp)欠啤，什么是歸一化，為什么要?dú)w一化自己查
到此我們得到了我們需要的數(shù)據(jù)了屋灌，并且對數(shù)據(jù)也處理好了洁段，下面就可以進(jìn)行差異分析了
1.2 設(shè)置分組 易出錯的步驟
差異分析顧名思義就是對比出差異，你需要告訴包共郭，那幾個樣本是control組哪些樣本是實(shí)驗(yàn)組祠丝，在這里為了大家都能用，所以只展示手動設(shè)置分組
Group = rep(c("treat","control"),times = c(13,9))
#為了大家分析不出錯除嘹，建議大家把實(shí)驗(yàn)組的sample調(diào)到前面写半，control組的sample調(diào)到后面，
#這樣就可以無腦運(yùn)行后面的代碼了尉咕，并且匹配設(shè)置group
#設(shè)置參考水平叠蝇，指定levels，對照組在前年缎，處理組在后
Group = factor(Group,
levels = c("control","treat")) #這里一定要這樣設(shè)置因子的位置
Group #因?yàn)楹竺娴拇a是按照這樣匹配的
1.3由于芯片測序的特殊性悔捶，我們需要拿到探針注釋，這樣才能將探針轉(zhuǎn)換為gene symbol 有很多種方法 這里只推薦為經(jīng)常用的三種 選一種即可
#第一種
#利用tinyarray的包
library(tinyarray)
find_anno(gpl_number) #假如你自己的數(shù)據(jù)单芜，請自行給gpl_number賦予平臺number
ids <- AnnoProbe::idmap('GPL570') #在這里就可以看到ids有兩列 一個是探針 一個是symbol
#方法2
#首先你需要知道 每一個平臺都有對應(yīng)的包蜕该。我們需要找到對應(yīng)的包，然后安裝洲鸠，利用包來轉(zhuǎn)換
#那我們怎么知道每一個平臺對應(yīng)的包呢堂淡，這里就要感謝jimmy大神了，大佬弄了一個list 里面有對應(yīng)關(guān)系
gpl_number
#我們知道我們的gpl_number后。
#http://www.bio-info-trainee.com/1399.html #打開這個鏈接绢淀。在里面找對應(yīng)
if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db")
library(hgu133plus2.db)
ls("package:hgu133plus2.db")
ids <- toTable(hgu133plus2SYMBOL)
head(ids)
# 方法3 #比較慢
#利用geo數(shù)據(jù)庫
#讀取GPL平臺的soft文件萤悴，按列取子集
##https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570
if(F){
library(GEOquery)
#注：soft文件列名不統(tǒng)一，活學(xué)活用更啄，有的表格里沒有symbol列稚疹，也有的GPL平臺沒有提供注釋
a = getGEO(gpl_number,destdir = ".")
b = a@dataTable@table
colnames(b)
ids2 = b[,c("ID","Gene Symbol")]
colnames(ids2) = c("probe_id","symbol")
ids2 = ids2[ids2$symbol!="" & !str_detect(ids2$symbol,"http:///"),]
}
jimmy大佬的list
ids
1.4 在這里強(qiáng)調(diào)居灯，芯片數(shù)據(jù)的差異分析只能用limma包來做＜牢瘛！怪嫌！
#差異分析义锥，用limma包來做
#需要表達(dá)矩陣和Group，不需要改
library(limma)
design=model.matrix(~Group) #根據(jù)group分組信息進(jìn)行貝葉斯檢驗(yàn)得出差異分析logFC以及p value
fit=lmFit(exp,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)
#為deg數(shù)據(jù)框添加幾列
#1.加probe_id列岩灭，把行名變成一列
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg)) #mutate函數(shù)給數(shù)據(jù)框添加一列probe_id 內(nèi)容是rownames(deg)
#2.加上探針注釋
ids = ids[!duplicated(ids$symbol),] #探針去重 自行了解
deg <- inner_join(deg,ids,by="probe_id") #將deg與ids合并 自動把deg中重復(fù)的probe_id給刪除
nrow(deg)
沒有將探針轉(zhuǎn)換為symbol的數(shù)據(jù)
轉(zhuǎn)換了的數(shù)據(jù)
這樣我們就有芯片測序結(jié)果的差異分析了拌倍。
總結(jié)：
1 判斷是否為芯片數(shù)據(jù) 是否取log
2需要測序平臺的number
3設(shè)置group，這里拋出一個問題噪径，就是多組比較怎么辦
4探針與symbol的轉(zhuǎn)換
2 seq數(shù)據(jù)
假如是用二代測序技術(shù)測的seq那么分析就跟芯片數(shù)據(jù)有很大的不同柱恤，為什么二者的分析的流程不能夠相互用這個你需要自己去google了https://zhuanlan.zhihu.com/p/375797068
在此特別強(qiáng)調(diào)！Ｕ野梗顺！
由于現(xiàn)在seq數(shù)據(jù)的三個差異分析的包對于輸入數(shù)據(jù)的要求都是要求count數(shù)據(jù)，假如從網(wǎng)上下載的數(shù)據(jù)不是count數(shù)據(jù)车摄，那么你的數(shù)據(jù)分析到此就很有可能截然而止了寺谤，不過有以下幾種辦法你可以嘗試：
1:發(fā)email問作者要（但是幾率很小，懂得都懂）
2:下載作者的fastq數(shù)據(jù)進(jìn)行上游分析（在你真的很想用這篇文章的數(shù)據(jù)的時候吮播，因?yàn)楹苡锌赡苣愕姆治鼋Y(jié)果得不到作者的結(jié)果）
3:硬來变屁，但是你的分析是不符合規(guī)范的！Ｒ夂荨粟关！
4:找到作者對上傳數(shù)據(jù)處理的手法，例如log环戈，那么你就log回去誊役，就會得到count數(shù)據(jù)了
ok，廢話不多說谷市，直接開始seq的數(shù)據(jù)的差異分析
再強(qiáng)調(diào)一遍蛔垢，需要輸入count數(shù)據(jù)！Ｆ扔啤Ｅ羝帷！
2.1 導(dǎo)入count數(shù)據(jù)和清洗數(shù)據(jù)
#seq數(shù)據(jù)無法用geoquery包，r直接聯(lián)網(wǎng)去下載數(shù)據(jù)
#你需要去手動下載supplement file里下載 并且記得看大小對不對得上
dat = read.table(”你的count數(shù)據(jù)“艺玲，check.names = F,row.names = 1,header = T)
#表達(dá)矩陣過濾
#需要過濾一下那些在很多樣本里表達(dá)量都為0或者表達(dá)量很低的基因括蝠。過濾標(biāo)準(zhǔn)不唯一。
nrow(dat) #過濾之前基因數(shù)量：
exp = dat[apply(dat, 1, function(x) sum(x > 0) > 0.5*ncol(exp)), ]#僅保留在一半以上樣本里表達(dá)的基因
nrow(exp)
最后你需要將數(shù)據(jù)處理為這個樣子
我們可以看到列名為sample名饭聚，行名為gene symbol忌警，再仔細(xì)觀察一下，會發(fā)現(xiàn)count數(shù)有的很大有幾千秒梳，而且都是整數(shù)法绵，都是整數(shù)，都是整數(shù)＠业狻Ｅ笃！
2.2 設(shè)置分組
跟芯片數(shù)據(jù)的分析一樣兴垦，要設(shè)置分組,我們還是手動設(shè)置吧
Group = rep(c("treat","control"),times = c(13,9)) #必須嚴(yán)格按照sample設(shè)置
#設(shè)置參考水平徙赢，指定levels，對照組在前探越，處理組在后
Group = factor(Group,
levels = c("control","treat")) #這里一定要這樣設(shè)置因子的位置
Group
這里得說一嘴狡赐，就是這里的代碼是需要你組內(nèi)有重復(fù)，假如你組內(nèi)沒有重復(fù)至少在這個代碼里會報(bào)錯（例如control組至少有兩個sample）
2.3 差異分析
對于seq數(shù)據(jù)差異分析有三個包DESeq2钦幔，edgeR枕屉，limma這里全部奉上三個包的代碼
library(DESeq2)
library(edgeR)
library(limma)
logFC_t = 2
pvalue_t = 0.05
#DESeq2做差異分析
colData <- data.frame(row.names =colnames(exp), #DESeq2要先創(chuàng)建一個數(shù)據(jù)框
condition=Group)
if(!file.exists(paste0(proj,"_dd.Rdata"))){
dds <- DESeqDataSetFromMatrix(
countData = exp,
colData = colData,
design = ~ condition) #第10行到第14行是創(chuàng)建一個有擬合分組的對象
dds <- DESeq(dds) #進(jìn)行差異分析
save(dds,file = paste0(proj,"_dd.Rdata")) #保存差異分析數(shù)據(jù)
}
load(file = paste0(proj,"_dd.Rdata"))
class(dds)
res <- results(dds, contrast = c("condition",rev(levels(Group)))) #提取做差異的結(jié)果
c("condition",rev(levels(Group)))
class(res)
DEG1 <- as.data.frame(res) #變成一個數(shù)據(jù)框
DEG1 <- DEG1[order(DEG1$pvalue),]
DEG1 = na.omit(DEG1)
head(DEG1)
#添加change列標(biāo)記基因上調(diào)下調(diào)
k1 = (DEG1$pvalue < pvalue_t)&(DEG1$log2FoldChange < -logFC_t);table(k1)
k2 = (DEG1$pvalue < pvalue_t)&(DEG1$log2FoldChange > logFC_t);table(k2)
DEG1$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
table(DEG1$change)
head(DEG1)
#edgeR----
dge <- DGEList(counts=exp,group=Group) #也是要一個表達(dá)矩陣和分組信息
dge$samples$lib.size <- colSums(dge$counts)
dge <- calcNormFactors(dge)
design <- model.matrix(~Group)
dge <- estimateGLMCommonDisp(dge, design)
dge <- estimateGLMTrendedDisp(dge, design)
dge <- estimateGLMTagwiseDisp(dge, design)
fit <- glmFit(dge, design)
fit <- glmLRT(fit)
DEG2=topTags(fit, n=Inf)
class(DEG2)
DEG2=as.data.frame(DEG2) #提取差異結(jié)果
head(DEG2)
k1 = (DEG2$PValue < pvalue_t)&(DEG2$logFC < -logFC_t);table(k1)
k2 = (DEG2$PValue < pvalue_t)&(DEG2$logFC > logFC_t);table(k2)
DEG2$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
table(DEG2$change)
###limma----
dge <- DGEList(counts=exp)
dge <- calcNormFactors(dge)
design <- model.matrix(~Group)
v <- voom(dge,design, normalize="quantile") #limma包對RNA-SEQ的話要VOOM標(biāo)準(zhǔn)化
fit <- lmFit(v, design)
fit= eBayes(fit)
DEG3 = topTable(fit, coef=2, n=Inf) #直接提取 不用as.data.frame
DEG3 = na.omit(DEG3)
k1 = (DEG3$P.Value < pvalue_t)&(DEG3$logFC < -logFC_t);table(k1)
k2 = (DEG3$P.Value < pvalue_t)&(DEG3$logFC > logFC_t);table(k2)
DEG3$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
table(DEG3$change)
head(DEG3)
上面的包選擇其中一個即可，或者你可以取他們的交集节槐，怎么判斷自己分析出來的差異結(jié)果是對的呢搀庶，這里說一句就是對于無處理count的seq數(shù)據(jù)logFC大概會在正負(fù)10左右，而且stable的基因遠(yuǎn)遠(yuǎn)大于上下調(diào)gene铜异，假如你想要更深程度知道每一個包的分析流程（更多個性話多流程）哥倔，你可以去學(xué)習(xí)每一個包的分析流程，大概都大差不差揍庄，按照要求輸入數(shù)據(jù)咆蒿，清洗數(shù)據(jù)，設(shè)置組別蚂子，差異分析沃测。
ok后面的就是簡簡單單的富集分析了.....
上面有寫的不對的地方，歡迎大家批評指正與討論
References：
https://zhuanlan.zhihu.com/p/375797068
https://blog.csdn.net/qazplm12_3/article/details/121134183
https://www.sohu.com/a/446981070_120380672
https://zhuanlan.zhihu.com/p/493445452
這里面有幾篇寫的很好Ｊ尘ァ５倨啤！
加油１鹩妗８矫浴惧互！
