R语言笔记

概览

R语言基础

R语言简介

R语言应用

谷歌学术2022年8月10日截图：

软件安装

如果电脑分成系统盘（C盘）和非系统盘（非C盘），先在非C盘的位置创建文件夹R，然后安装R和RStudio的时候选择安装在这个路径。

R

RStudio

rtools

rtools的作用可以简单理解为有时候有些R包的安装需要用rtools编译才能安装。

注意：安装时默认路径安装！！！

软件设置

快捷键

常见报错

R包管理

查找合适的包

R包管理

镜像问题

Bioconductor包

安装R包

R包迁移

文件操作

认识数据

数据清洗

工作目录

读写文本文件

读写EXcel

读写其他文件

数据结构

数据类型

数据结构

变量操作

数据索引

数据索引

判断类型

修改内容

数据操作

合并

排序

筛选

分组

计算

ggplot2基础

资料

绘图语法

  • 理解区分分类变量和连续变量；
  • +必须位于每行绘图代码的末尾；
  • 区分形状颜色和填充色；

image-20220826222040082

  • 不同的几何对象可以筛选数据；

image-20220826223036734

图层概念

主题修改

颜色修改

学习资源

坐标修改

分面操作

常见图形

学习资源

散点图

柱状图

堆积柱状图加上误差线

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
StackedBarErrorBar <- function(data, value, group, level) {
  StackedBarErrorBar.sub <- function(data, value, group, level) {
    res <- NULL

    col.name <- colnames(data)

    data %>%
      dplyr::rename(
        value = value,
        group = group
      ) -> data

    for (i in 1:length(level)) {
      data %>%
        dplyr::filter(group %in% level[i:length(level)]) %>%
        dplyr::group_by(group) %>%
        dplyr::mutate(mean.value = mean(value)) %>%
        dplyr::ungroup() %>%
        dplyr::select(group, mean.value) %>%
        dplyr::filter(!duplicated(group)) %>%
        dplyr::select(mean.value) -> res.1

      mean.value <- sum(res.1$mean.value)

      data %>%
        dplyr::filter(group == level[i]) %>%
        dplyr::mutate(mean.posi = mean.value) %>%
        dplyr::mutate(
          mean.value = mean(value),
          sd.posi = sd(value),
          n = n(),
          se = sd.posi / sqrt(n)
        ) %>%
        rbind(res) -> res
    }

    col.name <- c(col.name, "mean.posi", "mean.value", "sd.posi", "n", "se")

    res %>%
      magrittr::set_names(col.name) -> res

    return(res)
  }

  # contiction
  if ("grouped_df" %in% class(data)) {
    data %>%
      do(StackedBarErrorBar.sub(.,
        value = value,
        group = group,
        level = level
      )) -> res
  } else {
    StackedBarErrorBar.sub(
      data = data,
      value = value,
      group = group,
      level = level
    ) -> res
  }
  return(res)
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
df %>% 
  dplyr::filter(group == "干重") %>% 
  tidyr::pivot_longer(cols = 3:6) %>% 
  dplyr::select(-group)%>% 
  dplyr::group_by(treatment) %>% 
  StackedBarErrorBar(.,
                     value = "Value",
                     group = "name", 
                     level = c("叶","茎", "主根","须根")) %>% 
  dplyr::mutate(name = factor(name, level = c("叶","茎", "主根","须根")),
                label = "a") %>% 
  ggplot(aes(treatment, y = mean.value/n, fill = name)) +
  geom_bar(stat = "identity") +
  geom_errorbar(aes(ymin = mean.posi - se,
                    ymax = mean.posi + se),
                width = 0.2) +
  geom_text(aes(treatment,y = mean.posi + se + 2, label = label), check_overlap = TRUE) +
  geom_hline(yintercept = 62, color = NA) +
  scale_fill_aaas() +
  scale_y_continuous(expand = c(0,0)) +
  mytheme() +
  theme(legend.position = c(0.1, 0.85),
        legend.title = element_blank(),
        legend.background = element_blank())

直方图

箱线图

韦恩图

小提琴图

热图

基础统计

学习资料

线性模型

  • 求最小二乘均值并标注显著性

    1
    2
    3
    4
    5
    6
    7
    library(tidyverse)
    library(multcomp)
    library(emmeans)

    lm(Sepal.Length ~ Species, data = iris) %>% 
      emmeans(~Species) %>% 
      cld(alpha = 0.5, Letters = letters, adjust="tukey")

  • 什么是模型的p值：和空模型相比，是否显著。

  • 什么是模型的$R^2$：多少的因变量被模型解释。

  • glm模型不会返回p值和$R^2$，利用函数anova可以计算p值，函数nagelkerke可以计算伪$R^2$。

t检验

方差分析

卡方检验

wilcox检验

Fisher精确检验

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
# create a dataframe
df <- data.frame("cured" = c(60, 30), "noncured" = c(10, 25), row.names = c("treated", "nontreated"))
df
# output
           cured noncured
treated       60       10
nontreated    30       25

mosaicplot(df, color = TRUE)  

library(stats)
fisher.test(df)
# output

	Fisher Exact Test for Count Data

data:  df
p-value = 0.0002357
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
  1.983312 13.107997
sample estimates:
odds ratio 
  4.930093

# create a dataframe
df <- data.frame("no_disease" = c(40, 30, 25), "disease" = c(10, 20, 25), 
                    row.names = c("drugA", "drugB", "drugC"))
df
# output
      no_disease disease
drugA         40      10
drugB         30      20
drugC         25      25

mosaicplot(df, color = TRUE)  

library(stats)
fisher.test(df)
# output
	Fishers Exact Test for Count Data

data:  df
p-value = 0.005304
alternative hypothesis: two.sided

library(rstatix)
pairwise_fisher_test(as.matrix(df), p.adjust.method = "fdr")

# output
  group1 group2     n       p   p.adj p.adj.signif
* <chr>  <chr>  <dbl>   <dbl>   <dbl> <chr>       
1 drugA  drugB    100 0.0486  0.0729  ns          
2 drugA  drugC    100 0.00305 0.00915 **          
3 drugB  drugC    100 0.422   0.422   ns

贝叶斯统计

学习资料

基础分析

回归分析

相关性分析

机器学习

PCA

LDA

随机森林

参考帖子

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
#随机森林寻找 Biomarker 和环境驱动因子，以及气泡图 + 柱状图 + 热图组合可视化
otu<-read.csv("D://test/test_otu.csv",row.names = 1)
env<-read.csv("D://test/test_design.csv")
#链接：https://pan.baidu.com/s/10MnqI7AJEk3f_qVoPDtjyA 提取码：5nnb #已经更新
library(vegan)
library(randomForest)
library(vegan)
library(ggplot2)
library(reshape2)
#构建数据框
train_otu<-cbind(env$Group,t(otu[1:100,]))
train_otu<-as.data.frame(train_otu)
###### 构建模型
set.seed(123)
rf_train<- randomForest (train_otu$V1~., data = train_otu, ntree=1000,importance = TRUE) #ntree 可用来调整模型效果
rf_train
###### 挖掘预测模型中最重要的关键微生物 biomarker
df<-data.frame(importance(rf_train), check.names = FALSE)
###### heatmap of top 20 preditors
#筛选前 20 重要 OTU, 依据 Increase in MSE%
df<-df[order(df$MeanDecreaseAccuracy,decreasing = T),]
df<-df[1:20,]
#获取 top20otu 物种信息和丰度表
otu_top<-subset(otu,rownames(otu)%in%rownames(df))
pp<-otu_top[,1:24]
pp_new<-as.data.frame(cbind(rowSums(pp[,1:12])/12,rowSums(pp[,13:24])/12))
colnames(pp_new)<-c("AA","BB")
pp_new$AA<-as.numeric(as.character(pp_new$AA))
pp_new$BB<-as.numeric(as.character(pp_new$BB))
io<-pp_new%>% mutate(pp_new=row.names(.))%>%melt()
#组图 A Biomarker 在样品之间的丰度气泡图
ggplot(io,aes(variable,pp_new,fill=value))+
  geom_tile(color="grey2",fill="white",size=0.7)+
  geom_point(pch=21,size=8.2)+
  scale_fill_gradientn(colours =colorRampPalette(c('#2D6DB1', 'white', '#DC1623'))(111))+
  labs(x = NULL,y = NULL,fill="")+
  theme_minimal()+
 # theme(legend.position = "down",legend.text = element_text(size = 6))
  theme(axis.text=element_text(colour='black',size=9))+
  theme(axis.text.x =element_text(angle =90,hjust =0.9,vjust = 0.9,size=8))
##### 组图 b 柱状图展示 top 20 OTU 的重要性
dataframe<-cbind(df,rownames(df))
ggplot(dataframe, aes(x=MeanDecreaseAccuracy, y=rownames(df) )) + 
  geom_bar(stat = "identity", width = 0.7,position = "dodge",colour="black",fill="#2D6DB1") + 
  labs(x="Increase in MSE%", y="",size = 4)+
  theme_bw(base_line_size = 1.1,base_rect_size =1.1)+
  theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank())+
  theme(axis.text=element_text(colour='black'))+
  theme(legend.position = "right",legend.text = element_text(size = 6))
#组图 c 相关矩阵展示 top 20 OTU 和环境因子关联性
#合并理化因子和微生物
library(corrplot)
library(Hmisc)
library(reshape2)
library(RColorBrewer)
df3<-cbind(t(otu_top[,1:24]),env[,5:12])
#计算 Spearman 相关分析系数
spearman_df <- corr.test(df3, method = 'spearman', adjust = 'none')
spearman_r<- data.frame(spearman_df$r) 
spearman_r<-spearman_r[1:20,21:28]
spearman_p<- data.frame(spearman_df$p) 
spearman_p<-spearman_p[1:20,21:28]
#绘制相关分析图
#转化为 ggplot 画图需要的长列表
io1<-spearman_r%>% mutate(pp=row.names(.))%>%melt()
io2<-spearman_p%>% mutate(pp=row.names(.))%>%melt()
df<-cbind(io1,io2$value,io2$value)
#添加显著性，标记星号，实际需要自行添加 (根据相关分析显著性结果 Spearman、Pearson 的 p 值)，这里仅供画图，数据无意义。
colnames(df)<-c("pp","variable","value","sig","p")
df$sig[df$p<=0.05]<-"*"
df$sig[df$p<=0.01]<-"**"
df$sig[df$p<=0.001]<-"***"
df$sig[df$p>=0.05]<-""
#计算 Spearman、Pearson 的 r 值后，其范围为 - 1 到 1
ggplot(df,aes(variable,pp,fill=value))+
  geom_tile(color="grey2",fill="white",size=0.7)+
  geom_point(pch=22,size=7.2)+
  geom_text(aes(label=df$sig), color="gray2", size=4) +
  scale_fill_gradientn(colours =colorRampPalette(c('#2D6DB1', 'white', '#DC1623'))(111))+
  labs(x = NULL,y = NULL,fill="R value")+
  theme_minimal()+
  theme(axis.text=element_text(colour='black',size=9))+
  theme(axis.text.x =element_text(angle =90,hjust =0.9,vjust = 0.9,size=9))

组学分析

转录组

群体遗传学

GWAS

学习资料

微生物组

大数据挖掘

其他资料

PERMANOVA

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
library(vegan)

data(dune)
data(dune.env)

dune.dist <- vegdist(dune, method="bray")

# default test by terms

dune.div <- adonis2(dune ~ Management*A1, 
                    data = dune.env, 
                    permutations = 999, 
                    method="bray")

dune.div %>% 
  as.data.frame()

# overall tests

adonis2(dune ~ Management*A1, 
        data = dune.env, 
        permutations = 999, 
        method="bray",
        by = NULL)

结构方程模型

带线性混合效应模型的 SEM 构建、解释以及实际应用

图片

1
2
3
4
5
6
7
8
9
10
11
12
#加载包，导入数据
library(piecewiseSEM)#加载包
sampld <- read.table("clipboard",header=T)#读入数据
sample_psem <- psem( 
  lme4::lmer(MBN ~ Plantdiversity + N_addition +(1|site), data = sampld),
  lme4::lmer(PC1 ~ Plantdiversity + N_addition +(1|site), data = sampld), 
  lme4::lmer(Plantdiversity ~ N_addition+(1|site), data = sampld),
  lme4::lmer(N2O ~ N_addition+MBN+PC1+Plantdiversity + (1|site), data = sampld)
)#建模
summary(sample_psem)#提取建模结果

plot(sample_psem,digits = 2)

图片

R开发

shiny

golem开发shiny

主要笔记：

  • 创建项目

下面的代码会创建一个R package格式的文件夹，

1
golem::create_golem(path = "path/to/package")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
.
+-- DESCRIPTION
+-- dev
|   +-- 01_start.R
|   +-- 02_dev.R
|   +-- 03_deploy.R
|   \-- run_dev.R
+-- GOKEGG4TCMP.Rproj
+-- inst
|   +-- app
|   |   \-- www
|   |       \-- favicon.ico
|   \-- golem-config.yml
+-- man
|   \-- run_app.Rd
+-- NAMESPACE
\-- R
    +-- app_config.R
    +-- app_server.R
    +-- app_ui.R
    \-- run_app.R

主要的操作都在dev这个文件夹下完成。

  • 编辑DESCRIPTION文件

shinyWidgets

各种有意思的组件

shiny案例

appsilon

CSS优化shiny

服务器部署Rshiny注意事项

  • 应用文件夹应该放在/srv/shiny-server文件夹下才能生效，否则页面无法找到；
  • 可以把应用放在/opt/shiny-server/app/下，然后软连接到/srv/shiny-server.

其他知识点

关于距离

各种距离计算

点击访问

git知识点

添加服务器git服务器

1
2
3
git remote add sysfwq ssh://lixiang@xx.xx.xx.xx:2222/home/xx/gitppi/phdata.git
git push sysfwq
git push --set-upstream sysfwq master

snakemake流程

生物信息学
#生物信息学 #R
R语言笔记
https://lixiang117423.github.io/article/sharer/
作者
小蓝哥
发布于
2022年8月10日
许可协议