0.输入数据
在生信星球公众号聊天窗口回复exp777即可获取，里面的数据如下：
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.漂漂亮亮箱线图
这个是表达矩阵里的基因/探针或者细胞丰度矩阵里的细胞，画分组箱线图，需要提供的输入数据是表达矩阵和表达分组信息的因子。因子需要和表达矩阵的列一一对应。
draw_boxplot(exp[1:40,],Group,width = 1)
自带显著性检验了，默认kruskal.test，不想要那些不显著的行，可以把他们去掉：
draw_boxplot(exp[1:10,],Group,drop = T)
2.芯片表达矩阵的行名转换
从GEO下载下来的表达矩阵行名是探针名，找到探针注释，可以转换为基因名
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.热图，改了改配色和细节
scale_before 这个参数以后就不用了。这里画的是表达矩阵里标准差最大的2000个基因哦。
n = names(tail(sort(apply(exp, 1, sd)),2000))
draw_heatmap(exp[n,],
Group,legend = T,
annotation_legend = T)
这里的函数写的初衷是为了简化一些常见图表，省掉无数的复制粘贴，如果想要自定义一些我没有设置的参数，那就把代码拿下来自己编写调整，比如：
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>
函数的本质就是把参数代入进了这些代码，研究一下能搞定的~
下面的几个函数是原来的版本也有的~
4.PCA火山图韦恩图 通通可以拥有
draw_pca(exp,Group)
#内置数据集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")
搞定~
网友评论