由于太久没有分析过bulk数据了，前几天突然被要求分析一个文章的数据时有点乱了手脚，所以在这里想总结一下关于芯片数据与seq数据的差异分析流程。
假如你是小白，也刚刚了解NGS的冰山一角（我指的是转录组），就不要随随便便拿到数据去跑你那所谓的流程代码，因为很有可能是错的，还得意洋洋（我会鄙视你的）
所以这篇blog想要阐明的问题就是：芯片数据和seq数据不同的差异分析流程
首先声明，你得去了解什么是基因芯片测序和seq测序，这里我不讲，自行Google！
一个必要的好习惯：当别人给你一个bulk的转录组数据时，你首先需要做的事，问这个测序是什么方式测序（芯片or二代测序），因为不问你无法分析或者分析的也是大错特错！！！
ok，咱们就直接开始
1 芯片数据
假如你拿到是芯片数据那么就用下面的分析流程吧
首先你得需要一个表达矩阵吧，若你想分析文章的数据是芯片数据，你需要怎么下载呢，很简单用一个包
1.1下载和处理数据
#数据下载
rm(list = ls())
library(GEOquery)
gse_number = "GSE56649" #这里你只需要改成你要下载文献的GEO号
eSet <- getGEO(gse_number,
destdir = '.',
getGPL = F)。 #这个包会自动联网去GEO数据库下载该文章的数据
class(eSet)
length(eSet)
eSet = eSet[[1]]
#(1)提取表达矩阵exp
exp <- exprs(eSet) #表达矩阵在assay data中exprs
exp[1:4,1:4]
if(T){exp = log2(exp+1)}#避免有0值的存在 加1不影响大小（观察表达矩阵中有无取过log）
boxplot(exp) #观察是否差不多的箱长，如果不一样
if(T){exp = limma::normalizeBetweenArrays(exp)}
boxplot(exp) #这时候箱长就会一样的
#(2)提取临床信息
pd <- pData(eSet)
#(3)调整pd的行名顺序与exp列名完全一致
p = identical(rownames(pd),colnames(exp));p
if(!p) exp = exp[,match(rownames(pd),colnames(exp))]
#(4)提取芯片平台编号
gpl_number <- eSet@annotation
对上面代码块中关于表达矩阵的处理进行解释，
通过exp[1:4,1:4]判断数据长啥样，判断表达矩阵是否已经取过log，假如表达矩阵的数值范围在0-15，最大不超过20，那么就是取过log了，不要运行exp = log2(exp+1)，反之就要运行，这里解释一下问什么要取log，假如没有取log的话，差异logFC会大的离谱，但是你取了多次log的话，后面的logFC会很小，我们对表达值取log不仅是对logFC有作用，而且对后面画热图数值的归一化也很有用 然后我们就要归一化处理limma::normalizeBetweenArrays(exp)，什么是归一化，为什么要归一化自己查
到此我们得到了我们需要的数据了，并且对数据也处理好了，下面就可以进行差异分析了
1.2 设置分组 易出错的步骤
差异分析顾名思义就是对比出差异，你需要告诉包，那几个样本是control组哪些样本是实验组，在这里为了大家都能用，所以只展示手动设置分组
Group = rep(c("treat","control"),times = c(13,9))
#为了大家分析不出错，建议大家把实验组的sample调到前面，control组的sample调到后面，
#这样就可以无脑运行后面的代码了，并且匹配设置group
#设置参考水平，指定levels，对照组在前，处理组在后
Group = factor(Group,
levels = c("control","treat")) #这里一定要这样设置因子的位置
Group #因为后面的代码是按照这样匹配的
1.3由于芯片测序的特殊性，我们需要拿到探针注释，这样才能将探针转换为gene symbol 有很多种方法 这里只推荐为经常用的三种 选一种即可
#第一种
#利用tinyarray的包
library(tinyarray)
find_anno(gpl_number) #假如你自己的数据，请自行给gpl_number赋予平台number
ids <- AnnoProbe::idmap('GPL570') #在这里就可以看到ids有两列 一个是探针 一个是symbol
#方法2
#首先你需要知道 每一个平台都有对应的包。我们需要找到对应的包，然后安装，利用包来转换
#那我们怎么知道每一个平台对应的包呢，这里就要感谢jimmy大神了，大佬弄了一个list 里面有对应关系
gpl_number
#我们知道我们的gpl_number后。
#http://www.bio-info-trainee.com/1399.html #打开这个链接。在里面找对应
if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db")
library(hgu133plus2.db)
ls("package:hgu133plus2.db")
ids <- toTable(hgu133plus2SYMBOL)
head(ids)
# 方法3 #比较慢
#利用geo数据库
#读取GPL平台的soft文件，按列取子集
##https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570
if(F){
library(GEOquery)
#注：soft文件列名不统一，活学活用，有的表格里没有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,"///"),]
}
jimmy大佬的list
ids
1.4 在这里强调，芯片数据的差异分析只能用limma包来做！！！
#差异分析，用limma包来做
#需要表达矩阵和Group，不需要改
library(limma)
design=model.matrix(~Group) #根据group分组信息进行贝叶斯检验得出差异分析logFC以及p value
fit=lmFit(exp,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)
#为deg数据框添加几列
#1.加probe_id列，把行名变成一列
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg)) #mutate函数给数据框添加一列probe_id 内容是rownames(deg)
#2.加上探针注释
ids = ids[!duplicated(ids$symbol),] #探针去重 自行了解
deg <- inner_join(deg,ids,by="probe_id") #将deg与ids合并 自动把deg中重复的probe_id给删除
nrow(deg)
没有将探针转换为symbol的数据
转换了的数据
这样我们就有芯片测序结果的差异分析了。
总结：
1 判断是否为芯片数据 是否取log
2需要测序平台的number
3设置group，这里抛出一个问题，就是多组比较怎么办
4探针与symbol的转换
2 seq数据
假如是用二代测序技术测的seq那么分析就跟芯片数据有很大的不同，为什么二者的分析的流程不能够相互用这个你需要自己去google了https://zhuanlan.zhihu.com/p/375797068
在此特别强调！！！
由于现在seq数据的三个差异分析的包对于输入数据的要求都是要求count数据，假如从网上下载的数据不是count数据，那么你的数据分析到此就很有可能截然而止了，不过有以下几种办法你可以尝试：
1:发email问作者要（但是几率很小，懂得都懂）
2:下载作者的fastq数据进行上游分析（在你真的很想用这篇文章的数据的时候，因为很有可能你的分析结果得不到作者的结果）
3:硬来，但是你的分析是不符合规范的！！！
4:找到作者对上传数据处理的手法，例如log，那么你就log回去，就会得到count数据了
ok，废话不多说，直接开始seq的数据的差异分析
再强调一遍，需要输入count数据！！！！
2.1 导入count数据和清洗数据
#seq数据无法用geoquery包，r直接联网去下载数据
#你需要去手动下载supplement file里下载 并且记得看大小对不对得上
dat = read.table(”你的count数据“，check.names = F,row.names = 1,header = T)
#表达矩阵过滤
#需要过滤一下那些在很多样本里表达量都为0或者表达量很低的基因。过滤标准不唯一。
nrow(dat) #过滤之前基因数量：
exp = dat[apply(dat, 1, function(x) sum(x > 0) > 0.5*ncol(exp)), ]#仅保留在一半以上样本里表达的基因
nrow(exp)
最后你需要将数据处理为这个样子
我们可以看到列名为sample名，行名为gene symbol，再仔细观察一下，会发现count数有的很大有几千，而且都是整数，都是整数，都是整数！！！
2.2 设置分组
跟芯片数据的分析一样，要设置分组,我们还是手动设置吧
Group = rep(c("treat","control"),times = c(13,9)) #必须严格按照sample设置
#设置参考水平，指定levels，对照组在前，处理组在后
Group = factor(Group,
levels = c("control","treat")) #这里一定要这样设置因子的位置
Group
这里得说一嘴，就是这里的代码是需要你组内有重复，假如你组内没有重复至少在这个代码里会报错（例如control组至少有两个sample）
2.3 差异分析
对于seq数据差异分析有三个包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要先创建一个数据框
condition=Group)
if(!file.exists(paste0(proj,"_dd.Rdata"))){
dds <- DESeqDataSetFromMatrix(
countData = exp,
colData = colData,
design = ~ condition) #第10行到第14行是创建一个有拟合分组的对象
dds <- DESeq(dds) #进行差异分析
save(dds,file = paste0(proj,"_dd.Rdata")) #保存差异分析数据
}
load(file = paste0(proj,"_dd.Rdata"))
class(dds)
res <- results(dds, contrast = c("condition",rev(levels(Group)))) #提取做差异的结果
c("condition",rev(levels(Group)))
class(res)
DEG1 <- as.data.frame(res) #变成一个数据框
DEG1 <- DEG1[order(DEG1$pvalue),]
DEG1 = na.omit(DEG1)
head(DEG1)
#添加change列标记基因上调下调
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) #也是要一个表达矩阵和分组信息
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) #提取差异结果
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标准化
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)
上面的包选择其中一个即可，或者你可以取他们的交集，怎么判断自己分析出来的差异结果是对的呢，这里说一句就是对于无处理count的seq数据logFC大概会在正负10左右，而且stable的基因远远大于上下调gene，假如你想要更深程度知道每一个包的分析流程（更多个性话多流程），你可以去学习每一个包的分析流程，大概都大差不差，按照要求输入数据，清洗数据，设置组别，差异分析。
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
这里面有几篇写的很好！！！
加油！！！
