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)
![图片](http://image109.360doc.com/DownloadImg/2021/12/0910/235577351_2_2021120910450152)
自带显著性检验了,默认kruskal.test,不想要那些不显著的行,可以把他们去掉:
draw_boxplot(exp[1:10,],Group,drop = T)
![图片](http://image109.360doc.com/DownloadImg/2021/12/0910/235577351_3_20211209104501880)
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)
![图片](http://image109.360doc.com/DownloadImg/2021/12/0910/235577351_4_20211209104502569)
这里的函数写的初衷是为了简化一些常见图表,省掉无数的复制粘贴,如果想要自定义一些我没有设置的参数,那就把代码拿下来自己编写调整,比如:
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)
![图片](http://image109.360doc.com/DownloadImg/2021/12/0910/235577351_5_20211209104503287)
#内置数据集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)
![图片](http://image109.360doc.com/DownloadImg/2021/12/0910/235577351_6_20211209104504130)
x = list(Deseq2=sample(1:100,70),edgeR = sample(1:100,70),limma = sample(1:100,70))
draw_venn(x,'test')