學習及參考教程：
https://mp.weixin.qq.com/s/uLt-bYCdVcBH6mqoIG_r6w這個教程里面有個別地方不太合理幅骄，所以調整并備注了一些命令函數(shù)的意義程腹，以備后面復習或調用具壮。
http://www.reibang.com/p/f994631684b0
文中實驗檢測了三組數(shù)據(jù)吩愧，WT第美，AI和Y19F橙喘，每組3個樣本。所以一共有9個counts數(shù)據(jù)亚侠。直接下載counts數(shù)據(jù)沐寺，在windows系統(tǒng)的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)
#讀入數(shù)據(jù)
rawcount <- read.delim("counts.xls",header = T,sep="\t") #讀入制表符分隔的數(shù)據(jù)
rawcount <- rawcount[,1:10] #保留1-10列
table(!duplicated(rawcount[,1])) #table函數(shù)用來展示重復項的頻數(shù)。duplicated函數(shù)用來表示是否重復數(shù)據(jù)，要不是重復的返回為FALSE颤陶。!duplicated是取duplicated的反相颗管，原來是FALSE的，就變?yōu)門URE滓走。
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函數(shù)用來合并列
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數(shù)據(jù)框
express_cpm<-log2(cpm(filter_count)+1) #edgeR中的cpm函數(shù)將原始計數(shù)轉換為CPM之后取log垦江。得到的express_cpm是矩陣，而不是數(shù)據(jù)框搅方。
save(express_cpm,filter_count,file="filterGeneCountCpm.Rdata")
說明
1比吭，Y 叔的 clusterProfiler 包（v4.4.4）中有個函數(shù) bitr() ，可以把Ensembl_ID 轉 Symbol姨涡。原來的數(shù)據(jù)有58735行衩藤，生成id2symbol這個變量，只有35309行涛漂，也就是只有這么多基因找到了symbol名字赏表。
2，用到cbind和merge這樣的函數(shù)怖喻，可以對數(shù)據(jù)進行合并操作底哗。cbind必須具有相同的行數(shù)，就可以把列合并了锚沸。R中的merge函數(shù)類似于Excel中的Vlookup跋选，可以實現(xiàn)對兩個數(shù)據(jù)表進行匹配和拼接的功能。by.x,by.y:用于制定進行合并的兩個數(shù)據(jù)集的共同列哗蜈，不必須行數(shù)相同前标。比如id2symbol的ENSEMBL行只有35309行，而rawcount的GeneID有58735行距潘。all.y意思是y矩陣（即rawcount）的行應該全在輸出文件里炼列。因為y的行多于x的行，那么輸出文件里就會有一些行里面有空的數(shù)據(jù)NA音比。
3俭尖，!a %in% b表示a不在b中為ture。也就是列名不是"ENSEMBL","Ensmble","median"這幾個的都要保留洞翩。
4稽犁，rowSums(rawcount>0)獲取count數(shù)大于0的樣本個數(shù)。本例中有9個樣本骚亿，所以每行rowSums后結果必定是0-9之間的整數(shù)已亥。
5，cpm函數(shù)是edgeR中的一種基因定量方式来屠，count-per-millon虑椎。原始的表達量除以該樣本表達量的總和震鹉，在乘以一百萬就得到了CPM值 。
6捆姜，express_cpm<-log2(cpm(filter_count)+1)這個命令里面传趾，要用+1。因為log2(0)是無意義的娇未，如果count里面有某個數(shù)是0墨缘，那么就沒法計算了。所以這樣正合適log2(1)=0零抬，log2(2)=1, log2(4)=2....
二镊讼、WT和AI組兩組間比較
#展示W(wǎng)T和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) #檢查一下組別數(shù)量
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(), #去掉網(wǎng)格線
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函數(shù)建立一個新的數(shù)據(jù)框平夜，有兩列分別為expression和sample蝶棋。 expression是變量 exprSet1里的數(shù)據(jù)。對應sample列是把變量 exprSet1的列名標注出來忽妒，each變量為重復次數(shù)玩裙，即變量 exprSet1的行數(shù)。
2段直，ggplot畫圖,data參數(shù)指定繪圖的數(shù)據(jù)吃溅，aes參數(shù)指定圖的x軸，y軸是什么鸯檬，fill是填充的顏色怎么區(qū)分决侈。通常為fill=分組組名。
3喧务，theme_bw() 類似默認背景赖歌，調整為白色背景和淺灰色網(wǎng)格線，無邊框
4功茴，theme參數(shù)里面有比較多的選項
axis.text.x = element_text(angle = 90)) 橫坐標軸標簽的方向
theme參數(shù)挺復雜的庐冯，但是可以使得圖更美觀
panel.grid = element_blank() # 網(wǎng)格線不要。如果想要這個參數(shù)就刪掉坎穿。
panel.background = element_blank() #繪圖區(qū)背景不要展父。其實有了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)
#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)
說明
PCA函數(shù)是FactoMineR包里面的一個重要命令，用于分析娜饵。
fviz_pca_ind函數(shù)是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)
# 計算線性模型的參數(shù)
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函數(shù)是edgeR包中的用來構建一個含有組別標記的DGElist對象电湘，是一個可以包含多種內容和統(tǒng)計的列表隔节。用DGEList函數(shù)讀取count matrix數(shù)據(jù)，需要提供一個現(xiàn)成的matrix數(shù)據(jù)寂呛，而不是指望它能讀取單獨的文件怎诫。前面已經(jīng)有了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函數(shù)生成的對象DEG1有兩個部分：counts（數(shù)據(jù)矩陣），sample（那些樣本狗唉，分組情況）初烘，這個列表還可以添加更多內容，如下分俯。
calcNormFactors這個函數(shù)肾筐，是計算每個樣本的標準化參數(shù)，每個樣本都用不同的參數(shù)缸剪，以消除由于測序深度以及有效文庫大小不同等帶來的影響吗铐。這一步之后每個樣本norm.factors變得不同了。之前都是1杏节。
estimateGLMCommonDisp(x,design)：為所有基因都計算同樣的dispersion唬渗。
estimateGLMTrendedDisp(x,design)：根據(jù)不同基因的均值--方差之間的關系來計算dispersion典阵，并擬合一個trended model。
estimateGLMTagwiseDisp(x,design)：計算每個基因的dispersion镊逝，并利用經(jīng)驗性貝葉斯方法（empirical bayes）壓縮至Trend Didspersion壮啊。
上面三個參數(shù)依此計算后，會在DGE1這個對象后面添加好幾個參數(shù)撑蒜，都是離散相關的系數(shù)之類歹啼。
之后，glmFit函數(shù)是根據(jù)design進行擬合座菠，glmLRT函數(shù)利用likelihood ratio test（LRT）進行統(tǒng)計檢驗狸眼。glmFit 和 glmLRT 函數(shù)是配對使用的，用于 likelihood ratio test (似然比檢驗)辈灼。 contrast這個參數(shù)里面1和-1的前后位置無所謂份企。因為都是把WT作為對照組，AI作為實驗組巡莹。在分組設定design1里面WT在前面就是對照司志。
之后生成數(shù)據(jù)框DEG_edgeR1
最后進行篩選。
#火山圖
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函數(shù)為這個對象做圖。
