0.輸入數(shù)據(jù)
在生信星球公眾號(hào)聊天窗口回復(fù)exp777即可獲取涂召,里面的數(shù)據(jù)如下:
rm(list = ls())
load("step2output.Rdata")
ls()
## [1] "exp" "Group" "gse_number" "ids"
exp[1:4,1:4]
## GSM1366348 GSM1366349 GSM1366350 GSM1366351
## 1007_s_at 8.130087 7.671477 8.671782 7.868316
## 1053_at 8.932805 8.679543 8.625015 8.637085
## 117_at 9.383421 8.605809 9.462774 9.898573
## 121_at 7.916751 8.500635 8.258467 8.553656
dim(exp)
## [1] 54675 22
table(Group)
## Group
## control RA
## 9 13
head(ids)
## probe_id symbol
## 1 1053_at RFC2
## 2 117_at HSPA6
## 3 121_at PAX8
## 4 1255_g_at GUCA1A
## 5 1316_at THRA
## 6 1320_at PTPN21
我的小R包坠非,從寫著玩到用起來,日漸豐滿了果正。今天的更新炎码,1.5以上的版本可用哦,去我的github下載即可秋泳。
#devtools::install_github("xjsun1221/tinyarray",upgrade = F)
library(tinyarray)
1.漂漂亮亮箱線圖
這個(gè)是表達(dá)矩陣?yán)锏幕?探針或者細(xì)胞豐度矩陣?yán)锏募?xì)胞潦闲,畫分組箱線圖,需要提供的輸入數(shù)據(jù)是表達(dá)矩陣和表達(dá)分組信息的因子迫皱。因子需要和表達(dá)矩陣的列一一對(duì)應(yīng)歉闰。
draw_boxplot(exp[1:40,],Group,width = 1)
自帶顯著性檢驗(yàn)了,默認(rèn)kruskal.test,不想要那些不顯著的行新娜,可以把他們?nèi)サ簦?/p>
draw_boxplot(exp[1:10,],Group,drop = T)
2.芯片表達(dá)矩陣的行名轉(zhuǎn)換
從GEO下載下來的表達(dá)矩陣行名是探針名赵辕,找到探針注釋,可以轉(zhuǎn)換為基因名
n = trans_array(exp,ids)
n[1:4,1:4]
## GSM1366348 GSM1366349 GSM1366350 GSM1366351
## RFC2 8.932805 8.679543 8.625015 8.637085
## HSPA6 9.383421 8.605809 9.462774 9.898573
## PAX8 7.916751 8.500635 8.258467 8.553656
## GUCA1A 5.085221 2.414033 1.718570 4.311794
3.熱圖概龄,改了改配色和細(xì)節(jié)
scale_before 這個(gè)參數(shù)以后就不用了还惠。這里畫的是表達(dá)矩陣?yán)飿?biāo)準(zhǔn)差最大的2000個(gè)基因哦。
n = names(tail(sort(apply(exp, 1, sd)),2000))
draw_heatmap(exp[n,],
Group,legend = T,
annotation_legend = T)
這里的函數(shù)寫的初衷是為了簡化一些常見圖表私杜,省掉無數(shù)的復(fù)制粘貼蚕键,如果想要自定義一些我沒有設(shè)置的參數(shù),那就把代碼拿下來自己編寫調(diào)整衰粹,比如:
trans_array
## function (exp, ids, from = "probe_id", to = "symbol")
## {
## a = intersect(rownames(exp), ids[, from])
## message(paste0(length(a), " of ", nrow(exp), " rownames matched"))
## ids = ids[!duplicated(ids[, to]), ]
## exp = exp[rownames(exp) %in% ids[, from], ]
## ids = ids[ids[, from] %in% rownames(exp), ]
## exp = exp[ids[, from], ]
## rownames(exp) = ids[, to]
## message(paste0(nrow(exp), " rownames transformed after duplicate rows removed"))
## return(exp)
## }
## <bytecode: 0x561616ae6340>
## <environment: namespace:tinyarray>
函數(shù)的本質(zhì)就是把參數(shù)代入進(jìn)了這些代碼锣光,研究一下能搞定的~
下面的幾個(gè)函數(shù)是原來的版本也有的~?
4.PCA火山圖韋恩圖 通通可以擁有
draw_pca(exp,Group)
#內(nèi)置數(shù)據(jù)集deg
head(deg)
## logFC AveExpr t P.Value adj.P.Val B probe_id
## 1 5.780170 7.370282 82.94833 3.495205e-12 1.163798e-07 16.32898 8133876
## 2 -4.212683 9.106625 -68.40113 1.437468e-11 2.393169e-07 15.71739 7965335
## 3 5.633027 8.763220 57.61985 5.053466e-11 4.431880e-07 15.04752 7972259
## 4 -3.801663 9.726468 -57.21112 5.324059e-11 4.431880e-07 15.01709 7972217
## 5 3.263063 10.171635 50.51733 1.324638e-10 8.821294e-07 14.45166 8129573
## 6 -3.843247 9.667077 -45.87910 2.681063e-10 1.487856e-06 13.97123 8015806
## symbol change ENTREZID
## 1 CD36 up 948
## 2 DUSP6 down 1848
## 3 DCT up 1638
## 4 SPRY2 down 10253
## 5 MOXD1 up 26002
## 6 ETV4 down 2118
draw_volcano(deg,pkg = 4)
x = list(Deseq2=sample(1:100,70),edgeR = sample(1:100,70),limma = sample(1:100,70))
draw_venn(x,"test")
搞定~