生信代码:热图、箱线图、火山图、PCA、韦恩图通通可以拥有

↑↑↑   关注选刊说   ↑↑↑

您的SCI选刊助手!

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下载即可。

关于这个包之前写的介绍:

我写了一个R包,简化芯片的差异分析

我完善了那个R包,可以简化多组的差异分析啦

披着ggplot皮的pheatmap,深夜激动更新我的包

太好用了!批量生存分析加画图,一步到位,还支持最佳截点~

两个检验给ceRNA锦上添花

批量计算表达矩阵任意两个基因的相关性

#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")

学术交流群

小编,请备注: 单位+研究方向

免责声明:本文转载仅仅是出于传播信息的需要,并不意味着代表本公众号观点或证实其内容的真实性;如其他媒体、网站或个人从本公众号转载使用,须保留本公众号注明的“来源”,并自负版权等法律责任;作者如果不希望被转载,请与我们接洽。版权属于原作者

请大家点击右下角
科研菌分享给更多有需要的人
生物医学科研方法

手把手教你用 Excel 绘图 | 平均值对比图

2021-1-15 18:54:47

生物医学科研方法

散点图散点重叠过多时,怎么办?

2021-1-15 19:24:56