现在我们再解读一下第三张图，如果你对视频感兴趣，还是可以继续留邮箱，我们在圣诞节（明天） 统一发邮件给大家全部的视频云盘链接和配套代码哈！
本章节我们的视频审查员（刘博-二货潜）
本文的目录如下：
(1)MAplot：图a和图b
(2)差异基因韦恩图:UpSetR/VennDiagram
(3)两样本 log 2 FC 相关性散点图
万事开头前先看包的说明书： DESeq2 manual ，里面是讲的很详细很全面的，所以下面与文章不相关的图就不展示出来了。
rm(list = ls())F )'../figure-01-check-gene-expression/all.counts.id.txt' , header = T )7 :16 ]# 第一列为基因名，将其赋值给行名, 做韦恩图需要 1 ]# 查看前四行和前四列信息 1 :4 , 1 :4 ]library (stringr)# 要擅用 str_split 切割字符串，表示按照下划线 "_" 对列名进行切割，取第一列；即样本名 # 三组，每个样品一组，即 PhoKO、SppsKO、WT '_' , simplify = T )[, 1 ] ###################################################################### ################### Firstly for DEseq2 ##################### ###################################################################### # 一步运行 # 这里我们主要是得到差异基因中间部分就不细说 if (T ){# 赋值表达矩阵和分组信息为一个新的变量，方便以后直接调用 library (DESeq2)) # 加载包 # 构建一个表达矩阵 "qc_dispersions.png" , 1000 , 1000 , pointsize = 20 )"Dispersion plot" )#write.csv(exprMatrix_rlog,'exprMatrix.rlog.csv' ) # normalizedCounts2 <- counts(dds, normalized=T) # it's the same for the tpm value # we also can try cpm or rpkm from edgeR pacage #write.csv(exprMatrix_rpm,'exprMatrix.rpm.csv' ) "DEseq_RAWvsNORM.png" , height = 800 , width = 800 )0.7 )if (n.sample > 40 ) par(cex = 0.5 )1.2 )2 , 2 ))"expression value" , las = 2 )"expression value" , las = 2 )library (RColorBrewer)8 , "Dark2" )[1 :length(unique(group_list))])# Sample distance heatmap #install.packages("gplots",repos = "http://cran.us.") library (gplots)"qc-heatmap-samples.png" , w = 1000 , h = 1000 , pointsize = 20 )F , trace = "none" ,100 , "black" , "white" ),10 , 10 ), main="Sample Distance Matrix" )"group_list" , "SppsKO" , "WT" ))"group_list" ,"PhoKO" ,"WT" ))'deg_output.Rdata' )简化一下，如果不需要中间的信息，我们只需要差异基因的话，那么只需要运行以下代码：
if (T ){# 赋值表达矩阵和分组信息为一个新的变量，方便以后直接调用 library (DESeq2)) # 加载包 # 构建一个表达矩阵 # 下面我们得到 Spps 突变后的差异基因 "group_list" , "SppsKO" , "WT" )) # 注意这里是前面比后面，别把顺序搞反了，到时候一不注意结果就是反的。所以拿差异倍数对着原始 reads 比较一下。 # 按照 padj 值进行排序 # 将数据转变为 data.frame 数据框 # 去除包含 NA 值的行 # 下面我们得到 Pho 突变后的差异基因 "group_list" ,"PhoKO" ,"WT" ))# 将关键结果以 Rdata 形式保存到本地，下次如有需要就不需要重新用 DESeq2 计算了 'deg_output.Rdata' )好了，上面我们就得到了差异基因。
什么是 MAplot ?DESeq2 包说明书中的一段话：
In DESeq2, the function plotMA shows the log2 fold changes attributable to a given variable over the mean of normalized counts for all the samples in the DESeqDataSet . Points will be colored red if the adjusted p value is less than 0.1. Points which fall out of the window are plotted as open triangles pointing either up or down.
也就是说表示的是变化倍数 log2 (Fold change) 与所有样本标准化后的 counts 数的平均值之间的关系。那么怎么画呢 ？ 如果看过 DESeq2 说明书就知道这是 MA-plot 部分。 由于我们这里是将三组都放在一个 dds 中，所以我们得分别挑出 Pho 和 Spps 进行可视化。
resultsNames(dds)1 ] "Intercept" "group_list_SppsKO_vs_PhoKO" "group_list_WT_vs_PhoKO" lfcShrink 收缩 FC 三种方法如下（ 这里直接放原文）：
normal is the the original DESeq2 shrinkage estimator, an adaptive Normal distribution as prior. This is currently the default, although the default will likely change to apeglm in the October 2018 release given apeglm’s superior performance.
apeglm is the adaptive t prior shrinkage estimator from the apeglm package (Zhu, Ibrahim, and Love 2018).
ashr is the adaptive shrinkage estimator from the ashr package (Stephens 2016). Here DESeq2 uses the ashr option to fit a mixture of Normal distributions to form the prior, with method="shrinkage".
绘制 Spps
# apeglm 方法需要安装 apeglm 包 # BiocManager::install("apeglm") # ashr 方法同样需要安装额外依赖的包 # BiocManager::install("ashr") "group_list_SppsKO_vs_PhoKO" , type = "apeglm" )"group_list_SppsKO_vs_PhoKO" , type = "normal" )"group_list_SppsKO_vs_PhoKO" , type = "ashr" )1 ,3 ), mar=c(4 ,4 ,2 ,1 ))1 ,1e5 ); ylim <- c(-3 ,3 )"apeglm" )"normal" )"ashr" )绘制 Pho
Pho_resLFC <- lfcShrink(dds, coef = "group_list_WT_vs_PhoKO" , type = "apeglm" )"group_list_WT_vs_PhoKO" , type = "normal" )"group_list_WT_vs_PhoKO" , type = "ashr" )1 ,3 ), mar=c(4 ,4 ,2 ,1 ))1 ,1e5 ); ylim <- c(-3 ,3 )"apeglm" )"normal" )"ashr" )
好了，我们选取 normal （ 开心就好，你选哪个 ），来绘制图 a 和 b
par(mfrow=c(1 ,2 ), mar=c(4 ,4 ,2 ,1 ))"Spps_normal" )"Pho_normal" )
差异基因韦恩图:UpSetR/VennDiagram
我们下面将用两种方式来展示交集
第一种：我们使用 R 包 UpSetR 来绘制差异基因之间的韦恩图（ 多组时候，这种更加一目了然 ）
library (UpSetR)'deg_output.Rdata' )0.05 ,'stable' ,0 , 'up' , 'down' ))0.05 , 'stable' ,0 , 'up' , 'down' ))'up' ,])'down' ,])'up' ,])'down' ,])in % SppsKO_up),in % SppsKO_down),in % PhoKO_up),in % PhoKO_down))
第二种：我们使用 VennDiagram 来展示，也是就是文中那种图 # 这里直接 copy 琪同学的 # 链接: https://mp.weixin.qq.com/s/Pg0mjz7mD73atMnHz7jv1A # 导入本地字体，不然 `Arial` 字体会警告 library ("extrafont" )'deg_output.Rdata' )0.05 ,'stable' ,0 , 'up' , 'down' ))0.05 , 'stable' ,0 , 'up' , 'down' ))'up' ,])'down' ,])'up' ,])'down' ,])library (VennDiagram)"Up in phoKO" = PhoKO_up,"Down in phoKO" = PhoKO_down,"Up in SppsKO" = SppsKO_up,"Down in SppsKO" = SppsKO_down"Venn.png" , # 保存到当前目录 # stroke 颜色 形式 粗细 "#000000" , lwd = 3 , lty = 1 ,# 填充 透明度 # 颜色可以参考前一篇，使用 takecolor 自己取色 "#D3E7F0" , "#9FBEDB" , "#A0D191" , "#6DAE8A" ),0.50 ,# 里外字体大小形式 1.5 ,"Arial" ,"bold" ,1.5 ,0.2 ,"Arial" ,# 图像整体位置 分辨率 0.2 ,300 )
与文章趋势基本一致。可以看出 Spps 和 Pho 共同调控许多基因，说明这两基因有一定的关系。
load(file = 'deg_output.Rdata' )library (ggpubr)'SppsKO' , 'PhoKO' ,"reg.line" , # Add regressin line "blue" , fill = "lightgray" ), # Customize reg. line TRUE # Add confidence interval # Add correlation coefficient "pearson" , 10 , label.y = 10 ) # 相关系数显示位置
从上图可以看出，两者的相关系数高达0.53 ，这在两个基因间是算具有很强的相关的相关性了，更加佐证了上图的韦恩图的结果。
好了，此部分就到这里了。
你可能会需要： 广州专场（全年无休）GEO数据挖掘课，带你飞（1.11-1.12） 和 生信入门课全国巡讲2019收官--长沙站
