学习及参考教程：
https://mp.weixin.qq.com/s/uLt-bYCdVcBH6mqoIG_r6w这个教程里面有个别地方不太合理，所以调整并备注了一些命令函数的意义，以备后面复习或调用。
https://www.jianshu.com/p/f994631684b0

文中实验检测了三组数据，WT，AI和Y19F，每组3个样本。所以一共有9个counts数据。直接下载counts数据，在windows系统的R里面进行下面分析。先整体分析，后分别把AI和Y19F，与WT比较。然后看AIvsWT与Y19Fvs WT两种比对的交集，共性特性。

一、读入及格式转换

#先安装一些之前没有的包：
install.packages("FactoMineR")   #用来分析PCA或MCA等
install.packages("factoextra")  #用来可视化
install.packages("tidyverse") #一个基础包，含有ggplot2、dplyr。。很常用
install.packages("pheatmap")
install.packages("VennDiagram")
#清理变量加载各种包
rm(list = ls())
setwd("D:/work/NGSanalysis/rnaseq")
library(edgeR) 
library(clusterProfiler)
library(DESeq2)
library(FactoMineR)
library(factoextra)
library(org.Hs.eg.db)
library(stringr)
library(stringi)
library(tidyverse)
library(ggplot2)
library(patchwork)
library(pheatmap)
library(VennDiagram)
library(RColorBrewer)

#读入数据
rawcount <- read.delim("counts.xls",header = T,sep="\t")  #读入制表符分隔的数据
rawcount <- rawcount[,1:10]   #保留1-10列
table(!duplicated(rawcount[,1]))   #table函数用来展示重复项的频数。duplicated函数用来表示是否重复数据，要不是重复的返回为FALSE。!duplicated是取duplicated的反相，原来是FALSE的，就变为TURE。
rownames(rawcount) <- rawcount[,1]  #设定行名
rawcount <- rawcount[,-1] #去掉第一列
#转换gene id
id2symbol <- bitr(rownames(rawcount),   
                  fromType = "ENSEMBL", 
                  toType = "SYMBOL", 
                  OrgDb = org.Hs.eg.db)
rawcount <- cbind(GeneID=rownames(rawcount),rawcount)  #cbind函数用来合并列
rawcount <- merge(id2symbol,rawcount,
                  by.x="ENSEMBL",by.y="GeneID",all.y=T) 
rawcount=rawcount[!duplicated(rawcount$SYMBOL),]  #去除重复symbol行
rawcount <- drop_na(rawcount)  #去除含有NA值的行
rawcount <- rawcount[,!colnames(rawcount)%in%c("ENSEMBL","Ensmble","median")]#正式获得基因rawcount矩阵
rownames(rawcount) <- rawcount[,1]
rawcount <- rawcount[,-1]
keep <- rowSums(rawcount>0) >= floor(0.5*ncol(rawcount)) #至少在50%的样本中都有表达的基因要保留下来
filter_count <- rawcount[keep,] #获得filter_count数据框
express_cpm<-log2(cpm(filter_count)+1)  #edgeR中的cpm函数将原始计数转换为CPM之后取log。得到的express_cpm是矩阵，而不是数据框。
save(express_cpm,filter_count,file="filterGeneCountCpm.Rdata")

说明
1，Y 叔的 clusterProfiler 包（v4.4.4）中有个函数 bitr() ，可以把Ensembl_ID 转 Symbol。原来的数据有58735行，生成id2symbol这个变量，只有35309行，也就是只有这么多基因找到了symbol名字。
2，用到cbind和merge这样的函数，可以对数据进行合并操作。cbind必须具有相同的行数，就可以把列合并了。R中的merge函数类似于Excel中的Vlookup，可以实现对两个数据表进行匹配和拼接的功能。by.x,by.y:用于制定进行合并的两个数据集的共同列，不必须行数相同。比如id2symbol的ENSEMBL行只有35309行，而rawcount的GeneID有58735行。all.y意思是y矩阵（即rawcount）的行应该全在输出文件里。因为y的行多于x的行，那么输出文件里就会有一些行里面有空的数据NA。
3，!a %in% b表示a不在b中为ture。也就是列名不是"ENSEMBL","Ensmble","median"这几个的都要保留。
4，rowSums(rawcount>0)获取count数大于0的样本个数。本例中有9个样本，所以每行rowSums后结果必定是0-9之间的整数。
5，cpm函数是edgeR中的一种基因定量方式，count-per-millon。原始的表达量除以该样本表达量的总和，在乘以一百万就得到了CPM值 。
6，express_cpm<-log2(cpm(filter_count)+1)这个命令里面，要用+1。因为log2(0)是无意义的，如果count里面有某个数是0，那么就没法计算了。所以这样正合适log2(1)=0，log2(2)=1, log2(4)=2....

二、WT和AI组两组间比较

#展示WT和AI组的整体表达情况
filter_count1 <- filter_count[,c(1,2,3,4,5,6)]
express_cpm1 <- express_cpm[,c(1,2,3,4,5,6)]
colnames(filter_count1) <- c("WT1","WT2",'WT3',"Ai1","Ai2",'Ai3') #给样本命名
colnames(express_cpm1) <- colnames(filter_count1)
group1=rep(c("WT","Ai"),each=3)
group_list1=factor(group1,levels = c("WT","Ai"))
table(group_list1) #检查一下组别数量
exprSet1 <- express_cpm1
data1 <- data.frame(expression=c(exprSet1),
                   sample=rep(colnames(exprSet1),each=nrow(exprSet1)))
p1.1 <- ggplot(data = data1,aes(x=sample,y=expression,fill=sample))
plot1.1 <- p1.1 + geom_boxplot() +xlab(NULL) + ylab("log2(CPM)")+theme_bw()+ 
         theme(panel.background = element_blank(), #去掉背景色
         panel.grid = element_blank(),   #去掉网格线
         axis.text.x = element_text(face="bold",angle = 45,hjust = 1,color = 'black'),
         axis.title.x = element_blank(),
         legend.direction = "vertical",
         legend.title =element_blank())
plot1.1 #展示出来plot

说明
1，用data.frame函数建立一个新的数据框，有两列分别为expression和sample。 expression是变量 exprSet1里的数据。对应sample列是把变量 exprSet1的列名标注出来，each变量为重复次数，即变量 exprSet1的行数。
2，ggplot画图,data参数指定绘图的数据，aes参数指定图的x轴，y轴是什么，fill是填充的颜色怎么区分。通常为fill=分组组名。
3，theme_bw() 类似默认背景，调整为白色背景和浅灰色网格线，无边框
4，theme参数里面有比较多的选项
axis.text.x = element_text(angle = 90)) 横坐标轴标签的方向
theme参数挺复杂的，但是可以使得图更美观
panel.grid = element_blank() # 网格线不要。如果想要这个参数就删掉。
panel.background = element_blank() #绘图区背景不要。其实有了theme_bw() 选项，这个background没啥用啊。
axis.title.x = element_blank() # x轴标题不要 ，如果要可以用element_text("xxx")来标注。
legend.title =element_blank() #图例标题不要，如果要可以用element_text("xxx")来标注。
这几个后面element_blank()，表示为空
5，xlab(NULL)表示x轴标签不显示。ylab("log2(CPM)")表示y轴标签为log2(CPM)

image.png

#WT和AI两组比较的PCA图
data1.2 <- data.frame(t(exprSet1)) #t()为转置，就是横竖对换
data1.2 <- na.omit(data1.2) #删除有NA的行
dat_pca1 <- PCA(data1.2[,], graph = FALSE)
plot1.2 <- fviz_pca_ind(dat_pca1,
                     geom.ind = "point", # 只显示点，不显示文字
                     mean.point=F, #去除中心的一个大点
                     col.ind = group_list1, # 用不同颜色表示分组
                     palette = c("#00AFBB", "#E7B800"),
                     addEllipses = T, # 是否圈起来，少于4个样圈不起来
                     legend.title = "Groups") + theme_bw()
plot1.2 #展示这个plot
ggsave(plot1.2, filename = "AIvsWT_PCA.jpg",width = 6,height = 4,units = "in",dpi = 300) 

image.png

说明
PCA函数是FactoMineR包里面的一个重要命令，用于分析。
fviz_pca_ind函数是R包factoextra里面的命令，用于作图。
上述两个R包相辅相成，一起完成主成分分析及作图。
mean.point=F这个选项如果不加，会为每个组自动添加一个大的中心点，不怎么好看。
palette选项是指定各个分组的颜色。如果有三组可以设成这样：palette = c("#00AFBB", "#E7B800", "#FC4E07"),#三个组设置三种颜色

#WT和AI两组表达的差异分析(用edgeR包做分析）
exprSet1 <- filter_count1
design1 <- model.matrix(~0+factor(group_list1))  #建立分组文件
rownames(design1) <- colnames(exprSet1) #行名
colnames(design1) <- levels(factor(group_list1))  #列名
DEG1 <- DGEList(counts=exprSet1,  
                group=factor(group_list1))
DEG1 <- calcNormFactors(DEG1)
# 计算线性模型的参数
DEG1 <- estimateGLMCommonDisp(DEG1,design1)
DEG1 <- estimateGLMTrendedDisp(DEG1, design1)
DEG1 <- estimateGLMTagwiseDisp(DEG1, design1)
# 拟合线性模型
fit1 <- glmFit(DEG1, design1)
lrt1 <- glmLRT(fit1, contrast=c(1,-1)) 
# 提取过滤差异分析结果
DEG_edgeR1 <- as.data.frame(topTags(lrt1, n=nrow(DEG1)))
head(DEG_edgeR1)
fc <- 2.0
p <- 0.05
DEG_edgeR1$regulated <- ifelse(DEG_edgeR1$logFC>log2(fc)&DEG_edgeR1$PValue<p,
                              "up",ifelse(DEG_edgeR1$logFC<(-log2(fc))&DEG_edgeR1$PValue<p,"down","normal"))
write.csv(DEG_edgeR1,"DEG1.csv")  

说明
DGEList函数是edgeR包中的用来构建一个含有组别标记的DGElist对象，是一个可以包含多种内容和统计的列表。用DGEList函数读取count matrix数据，需要提供一个现成的matrix数据，而不是指望它能读取单独的文件。前面已经有了exprSet1这个变量。用class()查看，它是matrix，可以直接调用。
counts 就是这个matrix，group是因子，即分组信息，一定要是factor才行。
前面命令是这样的：
group1=rep(c("WT","Ai"),each=3)
group_list1=factor(group1,levels = c("WT","Ai"))
class(group_list1)可以查看到它是“factor”，而class(group1)则是“character”
DGEList函数生成的对象DEG1有两个部分：counts（数据矩阵），sample（那些样本，分组情况），这个列表还可以添加更多内容，如下。
calcNormFactors这个函数，是计算每个样本的标准化参数，每个样本都用不同的参数，以消除由于测序深度以及有效文库大小不同等带来的影响。这一步之后每个样本norm.factors变得不同了。之前都是1。
estimateGLMCommonDisp(x,design)：为所有基因都计算同样的dispersion。
estimateGLMTrendedDisp(x,design)：根据不同基因的均值--方差之间的关系来计算dispersion，并拟合一个trended model。
estimateGLMTagwiseDisp(x,design)：计算每个基因的dispersion，并利用经验性贝叶斯方法（empirical bayes）压缩至Trend Didspersion。
上面三个参数依此计算后，会在DGE1这个对象后面添加好几个参数，都是离散相关的系数之类。
之后，glmFit函数是根据design进行拟合，glmLRT函数利用likelihood ratio test（LRT）进行统计检验。glmFit 和 glmLRT 函数是配对使用的，用于 likelihood ratio test (似然比检验)。 contrast这个参数里面1和-1的前后位置无所谓。因为都是把WT作为对照组，AI作为实验组。在分组设定design1里面WT在前面就是对照。
之后生成数据框DEG_edgeR1

image.png
最后进行筛选。

#火山图
data1.3 <- DEG_edgeR1
plot1.3 <- ggplot(data=data1, aes(x=logFC, y=-log10(PValue),color=AIvsWT)) + 
  geom_point(alpha=0.5, size=1.8) + 
  theme_set(theme_set(theme_bw(base_size=20))) + 
  xlab("log2FC") + ylab("-log10(Pvalue)") +
  scale_colour_manual(values = c('blue','black','red'))+theme_bw()
plot1.3
ggsave(plot1.3, filename = "AIvsWT_volcano.jpg",width = 8,height = 6,units = "in",dpi = 300)
#热图
edgeR1_sigGene1 <- DEG_edgeR1[DEG_edgeR1$AIvsWT!="normal",]
edgeR1_sigGene1 <-rownames(edgeR1_sigGene1)
data1.4 <- express_cpm1[match(edgeR1_sigGene1,rownames(express_cpm1)),]
data1.4[1:4,1:4]
group1 <- data.frame(group=group_list1)
rownames(group1)=colnames(data1.4)
library(pheatmap)
plot1.4 <- pheatmap(data1.4,scale = "row",show_colnames =T,show_rownames = F, 
                 cluster_cols = T, 
                 annotation_col=group1,
                 main = "edgeR's DEG")
library(ggplotify)
plot1.4 <- as.ggplot(plot1.4)
plot1.4
ggsave(plot1.4,filename = "AIvsWT_heatmap.jpg",width = 8,height = 6,units = "in",dpi = 300 )

用ggsave命令把前面做出来的图保存在工作目录下。

三、Y19FvsWT两组间比较

#Y19FvsWT 
#  boxplot
filter_count2 <- filter_count[,c(1,2,3,7,8,9)]
express_cpm2 <- express_cpm[,c(1,2,3,7,8,9)]
colnames(filter_count2) <- c("WT1","WT2",'WT3',"Y19F1","Y19F2","Y19F3")
colnames(express_cpm2) <- colnames(filter_count2)
group2=rep(c("WT","Y19F"),each=3)
group_list2=factor(group2,levels = c("WT","Y19F"))
exprSet2 <- express_cpm2
data2.1 <- data.frame(expression=c(exprSet2),
                      sample=rep(colnames(exprSet2),each=nrow(exprSet2)))
p2.1 <- ggplot(data = data2.1,aes(x=sample,y=expression,fill=sample))
plot2.1 <- p2.1 + geom_boxplot(outlier.shape = NA) +xlab("sample") + ylab("log2(CPM)")+theme_bw()+ 
  theme(panel.background = element_blank(),
        panel.grid = element_blank(),  
        axis.text.x = element_text(face="bold",angle = 45,hjust = 1,color = 'black'),
        axis.title.x = element_blank(),
        legend.direction = "vertical",
        legend.title =element_blank())
plot2.1
ggsave(plot2.1, filename = "AIvsY19F_boxplot.jpg",width = 6,height = 4,units = "in",dpi = 300)
#AI$WT PCA
data2.2 <- data.frame(t(exprSet2))
data2.2 <- na.omit(data2.2)
dat_pca2 <- PCA(data2.2[,], graph = FALSE)
plot2.2 <- fviz_pca_ind(dat_pca1,
                        geom.ind = "point", 
                        mean.point=F,
                        col.ind = group_list1, 
                        palette = c("#00AFBB", "#E7B800"),
                        addEllipses = T, 
                        legend.title = "Groups") + theme_bw()
plot2.2
ggsave(plot2.2, filename = "AIvsY19F_PCA.jpg",width = 6,height = 4,units = "in",dpi = 300)
# DEG using edgeR
exprSet2 <- filter_count2
design2 <- model.matrix(~0+factor(group_list2))
rownames(design2) <- colnames(exprSet2)
colnames(design2) <- levels(factor(group_list2))
DEG2 <- DGEList(counts=exprSet2,  
                group=factor(group_list2))
DEG2 <- calcNormFactors(DEG2)
DEG2 <- estimateGLMCommonDisp(DEG2,design2)
DEG2 <- estimateGLMTrendedDisp(DEG2, design2)
DEG2 <- estimateGLMTagwiseDisp(DEG2, design2)
fit2 <- glmFit(DEG2, design1try)
lrt2 <- glmLRT(fit2, contrast=c(-1,1)) 
DEG_edgeR2<- as.data.frame(topTags(lrt2, n=nrow(DEG2)))
fc <- 2.0
p <- 0.05
DEG_edgeR2$AIvsY19F <- ifelse(DEG_edgeR2$logFC>log2(fc)&DEG_edgeR2$PValue<p,
                              "up",ifelse(DEG_edgeR2$logFC<(-log2(fc))&DEG_edgeR2$PValue<p,"down","normal"))
table(DEG_edgeR2$AIvsWT)

write.csv(DEG_edgeR1,"DEG2.csv")

#volcano
data2.3 <- DEG_edgeR2
plot2.3 <- ggplot(data=data2.3, aes(x=logFC, y=-log10(PValue),color=AIvsY19F)) + 
  geom_point(alpha=0.5, size=1.8) + 
  theme_set(theme_set(theme_bw(base_size=20))) + 
  xlab("log2FC") + ylab("-log10(Pvalue)") +
  scale_colour_manual(values = c('blue','black','red'))+theme_bw()
plot2.3
ggsave(plot2.3, filename = "AIvsY10F_volcano.jpg",width = 8,height = 6,units = "in",dpi = 300)
#heatmap
edgeR2_sigGene2 <- DEG_edgeR2[DEG_edgeR2$AIvsY19F!="normal",]
edgeR2_sigGene2 <-rownames(edgeR2_sigGene2)
data2.4 <- express_cpm2[match(edgeR2_sigGene2,rownames(express_cpm2)),]
data2.4[1:4,1:4]
group2 <- data.frame(group=group_list2)
rownames(group2)=colnames(data2.4)
library(pheatmap)
plot2.4 <- pheatmap(data2.4,scale = "row",show_colnames =T,show_rownames = F, 
                    cluster_cols = T, 
                    annotation_col=group2,
                    main = "edgeR's DEG")
library(ggplotify)
plot2.4 <- as.ggplot(plot2.4)
plot2.4
ggsave(plot2.4,filename = "AIvsY19F_heatmap.jpg",width = 8,height = 6,units = "in",dpi = 300 )

四、三个样本，在两两比较后的整合韦恩图

# 整体venn图，包含上调下调的
deg1 <- rownames(DEG_edgeR1[DEG_edgeR1$AIvsWT!="normal",])
up1 <- rownames(DEG_edgeR1[DEG_edgeR1$AIvsWT=="up",])
down1 <- rownames(DEG_edgeR1[DEG_edgeR1$AIvsWT=="down",])
deg2 <- rownames(DEG_edgeR2[DEG_edgeR2$Y19FvsWT!="normal",])
up2 <- rownames(DEG_edgeR2[DEG_edgeR2$Y19FvsWT=="up",])
down2 <- rownames(DEG_edgeR2[DEG_edgeR2$Y19FvsWT=="down",])
library(ggvenn)
deg<-list(`AIvsWT`=deg1,
          `Y19FvsWT`=deg2)
plot5 <- ggvenn(deg,
             show_percentage = F,
             stroke_color = "white",
             fill_color = c("#b2d4ec","#d3c0e2"),
             set_name_color = c("#ff0000","#4a9b83"))
plot5
plot5 <- as.ggplot(plot5)
ggsave(plot5, filename = "venn.jpg",width = 6,height = 4,units = "in",dpi = 300)
#上调，下调的韦恩图
up<-list(`up1`=up1,
         `up2`=up2)

plot6 <- ggvenn(up,
             show_percentage = F,
             stroke_color = "white",
             fill_color = c("#ffb2b2","#b2e7cb"),
             set_name_color = c("#ff0000","#4a9b83"))
plot6
plot6 <- as.ggplot(plot6)
ggsave(plot6, filename = "venn-up.jpg",width = 6,height = 4,units = "in",dpi = 300)
down<-list(`down1`=down1,
           `down2`=down2)

plot7 <- ggvenn(down,
             show_percentage = F,
             stroke_color = "white",
             fill_color = c("#b2d4ec","08a4ad"),
             set_name_color = c("#ff0000","#4a9b83"))
plot7
plot7 <- as.ggplot(plot7)
ggsave(plot7, filename = "venn-down.jpg",width = 6,height = 4,units = "in",dpi = 300)

说明
ggvenn是基于ggplot2的专门绘制韦恩图的R包。韦恩图，Venn diagram是用来展示集合的共性和特性。
先分别读取各个集合，比如第一组和第二组比较的差异基因集合deg1和deg2,以及上调基因集合up1/up2，下调基因集合down1/down2。每个集合都是一个列表，其中包含基因名。
[DEG_edgeR1$AIvsWT!="normal",]是取AIvsWT这一列为非normal的行。再把这些行的行名赋值给deg1这个变量。
在用list把需要做图的集合（可以是两个，三个或更多）赋值给一个对象。这里可以对每个集合进行命名，之后在图中就会显示该集合的名字。
最后用ggvenn函数为这个对象做图。

image.png