現(xiàn)在我們再解讀一下第三張圖，如果你對視頻感興趣，還是可以繼續(xù)留郵箱，我們在圣誕節(jié)（明天） 統(tǒng)一發(fā)郵件給大家全部的視頻云盤鏈接和配套代碼哈！
本章節(jié)我們的視頻審查員（劉博-二貨潛）
本文的目錄如下：
(1)MAplot：圖a和圖b
(2)差異基因韋恩圖:UpSetR/VennDiagram
(3)兩樣本 log 2 FC 相關(guān)性散點圖
萬事開頭前先看包的說明書： DESeq2 manual ，里面是講的很詳細很全面的，所以下面與文章不相關(guān)的圖就不展示出來了。
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 ){# 賦值表達矩陣和分組信息為一個新的變量，方便以后直接調(diào)用 library (DESeq2)) # 加載包 # 構(gòu)建一個表達矩陣 "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 ){# 賦值表達矩陣和分組信息為一個新的變量，方便以后直接調(diào)用 library (DESeq2)) # 加載包 # 構(gòu)建一個表達矩陣 # 下面我們得到 Spps 突變后的差異基因 "group_list" , "SppsKO" , "WT" )) # 注意這里是前面比后面，別把順序搞反了，到時候一不注意結(jié)果就是反的。所以拿差異倍數(shù)對著原始 reads 比較一下。 # 按照 padj 值進行排序 # 將數(shù)據(jù)轉(zhuǎn)變?yōu)?nbsp;data.frame 數(shù)據(jù)框 # 去除包含 NA 值的行 # 下面我們得到 Pho 突變后的差異基因 "group_list" ,"PhoKO" ,"WT" ))# 將關(guān)鍵結(jié)果以 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.
也就是說表示的是變化倍數(shù) log2 (Fold change) 與所有樣本標準化后的 counts 數(shù)的平均值之間的關(guān)系。那么怎么畫呢 ？ 如果看過 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 共同調(diào)控許多基因，說明這兩基因有一定的關(guān)系。
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 ) # 相關(guān)系數(shù)顯示位置
從上圖可以看出，兩者的相關(guān)系數(shù)高達0.53 ，這在兩個基因間是算具有很強的相關(guān)的相關(guān)性了，更加佐證了上圖的韋恩圖的結(jié)果。
好了，此部分就到這里了。
你可能會需要： 廣州專場（全年無休）GEO數(shù)據(jù)挖掘課，帶你飛（1.11-1.12） 和 生信入門課全國巡講2019收官--長沙站
