R實戰(zhàn) | NGS數(shù)據(jù)時間序列分析(maSigPro)
[圖片上傳失敗...(image-5355e7-1653975368931)]
[圖片上傳失敗...(image-533c51-1653975368931)]
一個答疑教程竖伯。
[圖片上傳失敗...(image-95c0e7-1653975368931)]
22
示例數(shù)據(jù)
#BiocManager::install('maSigPro')
library(maSigPro)
# 載入示例數(shù)據(jù)
data(data.abiotic)
data.abiotic[1:5,1:5]
data(edesign.abiotic)
head(edesign.abiotic)
> data.abiotic[1:5,1:5]
Control_3H_1 Control_3H_2 Control_3H_3 Control_9H_1 Control_9H_2
STMDF90 0.13735714 -0.3653065 -0.15329448 0.44754535 0.287476796
STMCJ24 NA NA NA NA NA
STMJH42 0.07864449 0.1002328 -0.17365488 -0.25279484 0.184855409
STMDE66 0.22976991 0.4740975 0.46930716 0.37101059 -0.004992029
STMIX74 0.14407618 -0.4801864 -0.07847999 0.05692331 0.013045420
> head(edesign.abiotic)
Time Replicate Control Cold Heat Salt
Control_3H_1 3 1 1 0 0 0
Control_3H_2 3 1 1 0 0 0
Control_3H_3 3 1 1 0 0 0
Control_9H_1 9 2 1 0 0 0
Control_9H_2 9 2 1 0 0 0
Control_9H_3 9 2 1 0 0 0
注意：data.abiotic是已經(jīng)標準化過的基因表達矩陣。
建立回歸模型
生成回歸矩陣(makeDesignMatrix)
design <- make.design.matrix(edesign.abiotic, degree = 2)
design$groups.vector
示例數(shù)據(jù)有三個時間的，故考慮二次回歸模型(degree = 2)姑尺。
> design$groups.vector
[1] "ColdvsControl" "HeatvsControl" "SaltvsControl" "Control"
[5] "ColdvsControl" "HeatvsControl" "SaltvsControl" "Control"
[9] "ColdvsControl" "HeatvsControl" "SaltvsControl"
尋找重要基因(p.vector)
F檢驗確定回歸方程的顯著性慎恒，采用BH的校正方式和悦，校正多重假設(shè)檢驗的p值赋焕。
校正后的p值小于p.vector的參數(shù)Q的基因就作為候選基因，進行下一步的分析强挫。通過fit$SELEC可以得到候選基因的表達量信息岔霸。
fit <- p.vector(data.abiotic, # 標準化的表達矩陣
design, # 實驗設(shè)計的矩陣 make.design.matrix 生成
Q = 0.05, # 顯著性水平
MT.adjust = "BH",
min.obs = 20
# 最低表達樣本數(shù) 不應(yīng)小于(degree+1)xGroups+1
)
fit$i # 顯著性基因的數(shù)量
fit$SELEC # 顯著性基因表達矩陣
尋找顯著性差異(T.fit)
上述的回歸方程是基于所有的自變量的組合構(gòu)建的，接下來就是通過逐步回歸法確定最佳的自變量組合俯渤。
tstep <- T.fit(fit, # p.vector結(jié)果
step.method = "backward",
alfa = 0.05) # 在逐步回歸中用于變量選擇的顯著性水平
在挑選最佳的自變量組合時呆细，通過每種自變量組合對應(yīng)的回歸模型的擬合優(yōu)度值R-squared來進行判斷，R-squared取值范圍為0到1八匠，數(shù)值越大絮爷，越接近1，回歸模型的效果越好坑夯。
獲取顯著性基因列表(get.siggenes)
sigs <- get.siggenes(tstep, # T.fit結(jié)果
rsq = 0.6, # 逐步回歸中的R-squared截至值
vars = "groups")
# vars參數(shù)有3種
# all 每個基因直接給出一個最佳的回歸模型
# groups 只給出不同實驗條件下相比control組中的差異基因
# each 會給出時間點和實驗條件的所有組合對應(yīng)差異基因列表
names(sigs)
names(sigs$sig.genes$ColdvsControl)
sigs$sig.genes$ColdvsControl$sig.profiles # 查看cold vs control的結(jié)果
結(jié)果可視化
韋恩圖(suma2Venn)
suma2Venn(sigs$summary[, c(2:4)]) # 左圖
suma2Venn(sigs$summary[, c(1:4)]) # 右圖
# 這個韋恩圖面積大小與數(shù)量不成比例 較普通
[圖片上傳失敗...(image-53f699-1653975368931)]
see.genes()
see.genes(sigs$sig.genes$ColdvsControl, # 差異基因表達矩陣
show.fit = T, # 是否顯示回歸擬合線(虛線)
dis =design$dis, # 回歸設(shè)計矩陣
cluster.method="hclust" , # 聚類方法
cluster.data = 1,
k = 9) # 聚類數(shù)目
[圖片上傳失敗...(image-e14690-1653975368931)]
這一步生成兩個圖柜蜈，如圖可分別查看缆娃。注意調(diào)整圖片顯示區(qū)域大小崇渗，以免報錯纵寝。
[圖片上傳失敗...(image-d192c2-1653975368931)]
[圖片上傳失敗...(image-6a30f8-1653975368931)]
PlotGroups()
選擇某一特定genes的表達進行可視化。
# 選取STMDE66基因
STMDE66 <- data.abiotic[rownames(data.abiotic)=="STMDE66", ]
PlotGroups (STMDE66,
edesign = edesign.abiotic)
[圖片上傳失敗...(image-ba61b6-1653975368931)]
# 添加回歸擬合線
PlotGroups (STMDE66,
edesign = edesign.abiotic,
show.fit = T,
dis = design$dis,
groups.vector = design$groups.vector)
[圖片上傳失敗...(image-cb99f-1653975368931)]
示例數(shù)據(jù)和代碼領(lǐng)取
詳見：https://mp.weixin.qq.com/s/MFdPrXrD_gErao-Jtsu6Fg
參考
Bioconductor - maSigPro(https://bioconductor.org/packages/release/bioc/html/maSigPro.html)
往期內(nèi)容
- (免費教程+代碼領(lǐng)取)|跟著Cell學(xué)作圖系列合集
- Q&A | 如何在論文中畫出漂亮的插圖必搞？
- Front Immunol 復(fù)現(xiàn) | 4. 使用estimate包評估腫瘤純度
- R繪圖 | 氣泡散點圖+擬合曲線
- 跟著 Cell 學(xué)作圖 | 桑葚圖(ggalluvial)
- R繪圖 | 對比條形圖+連線
- R繪圖 | 一幅小提琴圖的美化之旅
- R實戰(zhàn) | 給聚類加個圈圈(ggunchull)
- R繪圖 | 描述性統(tǒng)計常用圖(散點圖+柱狀圖+餅圖)
[圖片上傳失敗...(image-d6ce5c-1653975368931)]# R實戰(zhàn) | NGS數(shù)據(jù)時間序列分析(maSigPro)
[圖片上傳失敗...(image-50f56f-1653975368931)]
[圖片上傳失敗...(image-3758c5-1653975368931)]
一個答疑教程敬察。
[圖片上傳失敗...(image-1297de-1653975368931)]
22
示例數(shù)據(jù)
#BiocManager::install('maSigPro')
library(maSigPro)
# 載入示例數(shù)據(jù)
data(data.abiotic)
data.abiotic[1:5,1:5]
data(edesign.abiotic)
head(edesign.abiotic)
> data.abiotic[1:5,1:5]
Control_3H_1 Control_3H_2 Control_3H_3 Control_9H_1 Control_9H_2
STMDF90 0.13735714 -0.3653065 -0.15329448 0.44754535 0.287476796
STMCJ24 NA NA NA NA NA
STMJH42 0.07864449 0.1002328 -0.17365488 -0.25279484 0.184855409
STMDE66 0.22976991 0.4740975 0.46930716 0.37101059 -0.004992029
STMIX74 0.14407618 -0.4801864 -0.07847999 0.05692331 0.013045420
> head(edesign.abiotic)
Time Replicate Control Cold Heat Salt
Control_3H_1 3 1 1 0 0 0
Control_3H_2 3 1 1 0 0 0
Control_3H_3 3 1 1 0 0 0
Control_9H_1 9 2 1 0 0 0
Control_9H_2 9 2 1 0 0 0
Control_9H_3 9 2 1 0 0 0
注意：data.abiotic是已經(jīng)標準化過的基因表達矩陣鳞芙。
建立回歸模型
生成回歸矩陣(makeDesignMatrix)
design <- make.design.matrix(edesign.abiotic, degree = 2)
design$groups.vector
示例數(shù)據(jù)有三個時間的，故考慮二次回歸模型(degree = 2)贷盲。
> design$groups.vector
[1] "ColdvsControl" "HeatvsControl" "SaltvsControl" "Control"
[5] "ColdvsControl" "HeatvsControl" "SaltvsControl" "Control"
[9] "ColdvsControl" "HeatvsControl" "SaltvsControl"
尋找重要基因(p.vector)
F檢驗確定回歸方程的顯著性镊尺，采用BH的校正方式朦佩，校正多重假設(shè)檢驗的p值。
校正后的p值小于p.vector的參數(shù)Q的基因就作為候選基因庐氮，進行下一步的分析语稠。通過fit$SELEC可以得到候選基因的表達量信息。
fit <- p.vector(data.abiotic, # 標準化的表達矩陣
design, # 實驗設(shè)計的矩陣 make.design.matrix 生成
Q = 0.05, # 顯著性水平
MT.adjust = "BH",
min.obs = 20
# 最低表達樣本數(shù) 不應(yīng)小于(degree+1)xGroups+1
)
fit$i # 顯著性基因的數(shù)量
fit$SELEC # 顯著性基因表達矩陣
尋找顯著性差異(T.fit)
上述的回歸方程是基于所有的自變量的組合構(gòu)建的弄砍，接下來就是通過逐步回歸法確定最佳的自變量組合仙畦。
tstep <- T.fit(fit, # p.vector結(jié)果
step.method = "backward",
alfa = 0.05) # 在逐步回歸中用于變量選擇的顯著性水平
在挑選最佳的自變量組合時，通過每種自變量組合對應(yīng)的回歸模型的擬合優(yōu)度值R-squared來進行判斷音婶，R-squared取值范圍為0到1慨畸，數(shù)值越大，越接近1衣式，回歸模型的效果越好寸士。
獲取顯著性基因列表(get.siggenes)
sigs <- get.siggenes(tstep, # T.fit結(jié)果
rsq = 0.6, # 逐步回歸中的R-squared截至值
vars = "groups")
# vars參數(shù)有3種
# all 每個基因直接給出一個最佳的回歸模型
# groups 只給出不同實驗條件下相比control組中的差異基因
# each 會給出時間點和實驗條件的所有組合對應(yīng)差異基因列表
names(sigs)
names(sigs$sig.genes$ColdvsControl)
sigs$sig.genes$ColdvsControl$sig.profiles # 查看cold vs control的結(jié)果
結(jié)果可視化
韋恩圖(suma2Venn)
suma2Venn(sigs$summary[, c(2:4)]) # 左圖
suma2Venn(sigs$summary[, c(1:4)]) # 右圖
# 這個韋恩圖面積大小與數(shù)量不成比例 較普通
[圖片上傳失敗...(image-70965a-1653975368931)]
see.genes()
see.genes(sigs$sig.genes$ColdvsControl, # 差異基因表達矩陣
show.fit = T, # 是否顯示回歸擬合線(虛線)
dis =design$dis, # 回歸設(shè)計矩陣
cluster.method="hclust" , # 聚類方法
cluster.data = 1,
k = 9) # 聚類數(shù)目
[圖片上傳失敗...(image-a2527d-1653975368931)]
這一步生成兩個圖，如圖可分別查看碴卧。注意調(diào)整圖片顯示區(qū)域大小弱卡，以免報錯。
[圖片上傳失敗...(image-a09102-1653975368931)]
[圖片上傳失敗...(image-a4c051-1653975368931)]
PlotGroups()
選擇某一特定genes的表達進行可視化住册。
# 選取STMDE66基因
STMDE66 <- data.abiotic[rownames(data.abiotic)=="STMDE66", ]
PlotGroups (STMDE66,
edesign = edesign.abiotic)
[圖片上傳失敗...(image-1045d2-1653975368931)]
# 添加回歸擬合線
PlotGroups (STMDE66,
edesign = edesign.abiotic,
show.fit = T,
dis = design$dis,
groups.vector = design$groups.vector)
[圖片上傳失敗...(image-2d7f07-1653975368931)]
示例數(shù)據(jù)和代碼領(lǐng)取
點贊婶博、在看 本文，分享至朋友圈集贊20個并保留30分鐘荧飞，截圖發(fā)至微信mzbj0002領(lǐng)取凡人。
木舟筆記2022年度VIP可免費領(lǐng)取。
木舟筆記2022年度VIP企劃
權(quán)益：
-
2022年度木舟筆記所有推文示例數(shù)據(jù)及代碼(在VIP群里實時更新)垢箕。
<img src="https://mzbj-1312199408.cos.ap-shanghai.myqcloud.com/image-20220529221704900.png" alt="image-20220529221704900" style="zoom:67%;" />
木舟筆記科研交流群划栓。
半價購買
跟著Cell學(xué)作圖系列合集(免費教程+代碼領(lǐng)取)|跟著Cell學(xué)作圖系列合集兑巾。
收費：
**99￥/人**条获。可添加微信：`mzbj0002` 轉(zhuǎn)賬蒋歌，或直接在文末打賞帅掘。
<img src="https://mzbj-1312199408.cos.ap-shanghai.myqcloud.com/image-20220529221602758.png" alt="image-20220529221602758" style="zoom:50%;" />
參考
Bioconductor - maSigPro(https://bioconductor.org/packages/release/bioc/html/maSigPro.html)
往期內(nèi)容
- (免費教程+代碼領(lǐng)取)|跟著Cell學(xué)作圖系列合集
- Q&A | 如何在論文中畫出漂亮的插圖？
- Front Immunol 復(fù)現(xiàn) | 4. 使用estimate包評估腫瘤純度
- R繪圖 | 氣泡散點圖+擬合曲線
- 跟著 Cell 學(xué)作圖 | 桑葚圖(ggalluvial)
- R繪圖 | 對比條形圖+連線
- R繪圖 | 一幅小提琴圖的美化之旅
- R實戰(zhàn) | 給聚類加個圈圈(ggunchull)
- R繪圖 | 描述性統(tǒng)計常用圖(散點圖+柱狀圖+餅圖)
[圖片上傳失敗...(image-fc3c55-1653975368931)]
