友情提示：一定要有R语言的基础，建议学习生信技能树的《生信人这样这样学R语言》之后（最好再把中级20题做完），再来学习本系列视频！
总的来说，理解了R语言中级20题之后，这一块的内容学起来就轻松了很多。
1.通读文献阅读及规律
这一块儿主要就是在你感兴趣的文章里，找到其测序信息的GSE号，在GEO数据库中直接检索该文件号。
2.了解GEO数据库
比如我们对GSE42872 这个测序结果比较感兴趣，我们现在GEO Database中搜索该文件号，点击进入。
我们可以通过原文章了解实验设计及摘要，也可以在GEO Database中浏览。原文章是比较简单的一个实验设计思路，直接用药物vemurafenib或空白对照处理了黑色素瘤细胞后去做RNA seq，每组3个对照。
再向下看，就是测序信息了，这里要说明三个号：
数据集：GSE(GEO series)，也就是本次查询的号，是一个实验设计下的测序结果的集合，一个数据集下有多个样本，一篇文章也可能会出现多个数据集。
平台：GPL(GEO platform)，测序公司所提供的的测序系统，在这次实验中运用的是 Affymetrix Human Gene 1.0 ST Array。
样本：GSM(GEO sample)，每个样本的测序结果都会有一个GSM号。在这个实验中一共有6个样本，3个实验组，3个对照组。
3.数据下载的3种方法
①RAW data处下载，不推荐，不好分析。旧版本芯片用readaffy包来处理，illumina用iumiR.batch来处理。
②下载matrix表达矩阵，可直接导入R。（我比较喜欢这种方法）
gset=read.table('Downloads/GSE42872_series_matrix.txt.gz',sep='\t',comment.char='!',header = T)#读入，txt文件的分隔符为\t，前面会有一些以！开头的描述信息，删去这些，表头要有
class(gset)
str(gset)
rownames(gset)=gset[,1]#原本的第一列为探针名，现在将行名改为探针名
gset=gset[,-1]#删去原来的第一列
③在R中直接用GSE号下载：
source("[http://www.bioconductor.org/biocLite.R)](http://www.bioconductor.org/biocLite.R))
biocLite("GEOquery")
library(GEOquery)
gset<-getGEO("GEO number",GSEMatrix=TURE,AnnotGPL=TURE，destdir=".")
需要注意，如果是第二种方法下载，导入的直接是dataframe，第三种方法下载的是对象，里面还包含了描述信息，但是getGEO会帮助提取其中的表达矩阵。
4.ID转换技巧大全
同R中级20题中的一样，不过需要注意的是该套测序的探针用的是哪套注释信息，大家可以从platforms这里查看：
可以在网上找注释信息对应合集，如下述链接，在其中找GPL6244 [HuGene-1_0-st]对应的R包是hugene10sttranscriptcluster，下载该注释包。
https://blog.csdn.net/weixin_40739969/article/details/103186027
#探针注释信息下载
BiocManager::install('hugene10sttranscriptcluster.db') #下载，注意在包的后面手动加上.db
library('hugene10sttranscriptcluster.db')#载入
ls("package:hugene10sttranscriptcluster.db")#查看，找到后面是SYMBOL的，复制名字
ids=toTable(hugene10sttranscriptclusterSYMBOL)#提取注释信息
table(table(ids$symbol)>1)#统计有多少个基因对应一个以上的探针
table(sort(table(ids$symbol)))#统计有1、2、3等个探针数的基因分别有多少个，大多数基因只有一个探针
下载好之后，就开始将探针名更改为对应的基因名。先判断探针有无对应的基因名，然后对有多个探针的基因保留平均数或中位数或其他指标作为该基因的表达量。
#更改ID
rownames(exprSet)%in% ids$probe_id#a %in% table a是否位于table中，返回T or F。这里判断exprSet的行名是否位于ids的probe_id这一列中
table(rownames(exprSet)%in% ids$probe_id)#统计T、F的个数，有13483个探针不存在对应的基因名
n_exprSet=exprSet[!(rownames(exprSet)%in% ids$probe_id)==F,]#将包含在探针包里的探针测序数据新建一个数据框，这里很奇怪，因该是满足条件的包括进来，这里却填F的时候会包括进来
tmp = by(n_exprSet,ids$symbol,
function(x)
rownames(x)[which.max(rowMeans(x))] )#by函数，将对象按行或某种方式分为一个个小的子集，对每个子集进行操作。这里将对n_exprSet的行按与ids中的symbol的对应关系，将同一个symbol的探针放在一起，然后对子集中的每一行取平均值，返回每个子集中平均值大的那一行的行名（探针名）
probes = as.character(tmp)#定义为字符型
View(probes)#这里就是样本中相同基因平均值最大的探针
exprSet=exprSet[rownames(exprSet) %in% probes ,]#对exprSet的行进行操作，若该行的行名在probes中出现了则保留下来
dim(exprSet)#现在为18821个探针在6个样本中的表达量
View(exprSet)
rownames(exprSet)=ids[match(rownames(exprSet),ids$probe_id),2]#将exprSet的行名与ids的probe_id相互匹配，若相同，则exprSet的行名为ids的第二列
View(exprSet)
5.了解你的表达矩阵
通过一些简单的分析查看自己的表达矩阵是否合乎逻辑、前面的处理有无错误。
①查看管家基因的表达量，和表达量的中位数（boxplot）作比较，管家基因表达量一般比较高。
GAPDH的表达量均在14以上，β-ACTIN的表达量均在13以上，而表达量的中位数是在8到9之间，关键基因的表达量在所有表达量中处于较高的水平，这一项是合格的。
exprSet['GAPDH',]#查看GAPDH在六个样本中的表达量
exprSet['ACTB',]#查看β-actin在六个样本中的表达量
boxplot(exprSet)#画六个样本表达量的箱型图
②看表达矩阵分布图，查看所有样本表达量的中位数，应该是处于差不多水平的才可以进行比较。
六个样本的表达量中位数基本处于同一水平线上。
③画hclust图，查看同一分组是否在一起。
首先处理数据：
#添加样本分组信息
group_list=c(rep('control',3),rep('case',3))#新建list，根据GEO描述得知前三个为对照组、后3个为实验组
group_list
#将宽数据变为长数据
library(reshape2)#使用该包中的melt函数
a=rownames(exprSet)#这里需要操作一波，在实践中出现的，我单独写一篇专门讲述melt函数
a=as.data.frame(a)
exprSet=cbind(a,exprSet)
exprSet=exprSet[,-1]
exprSet1=cbind(a,exprSet)
exprSet_L=melt(exprSet1)#处理数据，将多维数据转变为一维数据，即宽数据变为长数据
colnames(exprSet_L)=c('probe','sample','value')#每一列依次为探针（基因）、样本来源、表达量
exprSet_L$group=rep(group_list,each=nrow(exprSet))
head(exprSet_L)
聚类分析，相同处理的样本应该是聚在一起的：
#聚类分析
nodePar <- list(lab.cex = 0.6, pch = c(NA, 19),
cex = 0.7, col = "blue")
hc=hclust(dist(t(exprSet)))
par(mar=c(5,5,5,10))
plot(as.dendrogram(hc), nodePar = nodePar, horiz = TRUE)
④PCA主成分分析，相同处理的样本也是聚在一起的：
#主成分分析
pc <- prcomp(t(exprSet),scale=TRUE)
pcx=data.frame(pc$x)
pcr=cbind(samples=rownames(pcx),group_list, pcx)
p=ggplot(pcr, aes(PC1, PC2))+geom_point(size=5, aes(color=group_list)) +#能最大反映样本差异性的两个成分（PC1、PC2）
geom_text(aes(label=samples),hjust=-0.1, vjust=-0.3)#label=samples可以加上样本名称
print(p)
6.差异分析（limma包），火山图、热图的绘制及美化
limma包分析时需要提前准备好三个矩阵：表达矩阵、分组矩阵、比较矩阵。表达矩阵已经处理好了，下面只需制作分组矩阵（用0、1表示实验组和对照组的样本）和比较矩阵（谁和谁比，实验组比上对照组）。
library(limma)#就是给表达矩阵、分组矩阵、比较矩阵，limma给你分析差异结果
design=model.matrix(~0+factor(group_list))#分组矩阵
design#分组矩阵制作好了，用0和1代替
colnames(design)=levels(factor(group_list))#给分组矩阵的列名用分组信息来表示
rownames(design)=colnames(exprSet)#行名用样本名来表示
contrast.matrix=makeContrasts(case-control,levels = design)
contrast.matrix #比较矩阵，makeContrasts用来告诉谁比谁
fit=lmFit(exprSet,design)#lmFit函数为线性模型拟合，给表达矩阵和分组矩阵，输出
fit2 <- contrasts.fit(fit, contrast.matrix) #把上一步的线性拟合矩阵矩阵和比较矩阵给它，计算差异和标准误
fit2 <- eBayes(fit2)#通过eBayes将标准误调整到一个相同的水平再来计算t、logFC等值
tempOutput = topTable(fit2, coef=1, n=Inf)#用topTable提取
nrDEG = na.omit(tempOutput) #该函数是删去缺少某些值的对象
head(nrDEG)
write.csv(nrDEG,"limma_notrend.results.csv",quote = F)#将差异分析结果保存到本地
检查一下差异分析的结果：
exprSet['CD36',]#查看表达情况
#热图画一下前25个基因，检验是否有差异
choose_gene=head(rownames(nrDEG),25)
choose_matrix=exprSet[choose_gene,]
choose_matrix=t(scale(t(choose_matrix)))
pheatmap::pheatmap(choose_matrix)#可以看到这些基因的表达量在试验和对照组中差异还是很大的
画火山图（和R20题的操作一样）：
#带标题、上下调基因的火山图
logFC_cutoff <- with(nrDEG,mean(abs( logFC)) + 2*sd(abs( logFC)) )
nrDEG$change = as.factor(ifelse(nrDEG$P.Value < 0.05 & abs(nrDEG$logFC) > logFC_cutoff,
ifelse(nrDEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
)#设置阈值
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',])
)#图片标题
library(ggplot2)
g = ggplot(data=nrDEG, aes(x=logFC, y=-log10(P.Value), color=change)) +
geom_point(alpha=0.4, size=1.75) +
theme_set(theme_set(theme_bw(base_size=20)))+
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')) ## corresponding to the levels(res$change)
print(g)
7.基因集的富集分析（注释）KEGG
富集分析或注释等通俗来讲就是这些差异基因处于哪些通路或者是功能中，即造成了哪些通路、功能发生改变。这里用到的是超几何分布的原理，简单介绍一下：假设一共有20000个基因，差异基因有1000个，也就是说每个基因被抽中的概率为1/20。现在有一个通路有100个基因，若随机抽，则差异基因中属于该通路的有5个，但这1000个基因中有20个基因都是属于这个通路，这就说明该通路是有问题的。
在用clusterProfiler包进行分析时，输入的ID类型为ENTREZID，但之前我们用的是SYMBOL，需要转换一下。ENTREZID即gene id。
#ID转换
gene=head(rownames(nrDEG),1000)#我们在这里只取前1000个差异基因来做，具体数字看自己的需求
library(clusterProfiler)#之前的gene ID是symbol，现在要用ENTREZID，所以需要进行转换
library(org.Hs.eg.db)
gene.df <- bitr(gene, fromType = "SYMBOL",
toType = c("ENSEMBL","ENTREZID"),
OrgDb = org.Hs.eg.db)#从SYMBOL转换为ENTREZID和ENSEMBL，可根据自己的需要转换
head(gene.df)#检查几个看看做得对不对
#分析
kk <- enrichKEGG(gene = gene.df$ENTREZID,
organism = 'hsa',
pvalueCutoff = 0.05)
head(kk)[,1:6]
如第一个通路是和DNA复制相关的，这条通路中一共有36个基因，有18个基因都出现在差异基因中了，说明该项实验中，药物处理后DNA的复制发生了很大很大的变化。
超几何分布检验的弊端就是只知道哪些通路发生了改变，但不知道是哪些基因在其中起到了比较重要的作用，这个时候可以用GSE进行分析（后续单独写一篇）。
