R语言绘图总汇 |覆盖教程所有图形 | 一

小杜的生信筆記

小杜的生信筆記，知乎、公众号同名。一名森信小白同学！！

前言

加入知乎已经快1年了（PS：还差4个月），自己也分享了一些教程。也获得了我们童靴的支持和点赞。（PS：自己仍然是个搬运工，哈哈哈哈哈）。现在回想，虽然分享了一些教程，但我们的同学找的时候仍是个比较麻烦的过程，那么为什么不可以做个总结呢？因此，这个教程就是总结一下前面教程的所有的绘图或分析的教程。最终目的，就是：方便你我他！

OK！

这就是这篇博客的初衷，会持续更新！

篇幅也会比较长，但是会有目录。 你可以直接跳转到你想要的内容，基本每个教程多会有实际数据 。 老粉都知道如何获取。如果不知道的童鞋，你可以关注我的公众号，每个教程会有相应的关键词，你回复即可哦！

--- Du

R语言中如何更改R包安装路径

R中包的安装位置默认是在C盘，经常出现C盘不够用，那我们就需要吧软件安装在其他盘中，包括R语言中R。

第一步：

## 我们操作在Rstuido中
## 新建一个初始启动文件
file.edit('~/.Rprofile')

第二步：

## 更改R包默认安装路径
.libPaths("D:/软件安装/R_libaray")   ## 输入你的路径
myPaths <- c("D:/软件安装/R_libaray")
.libPaths(myPaths)

现在就设置到你的R包安装路径，保存.Rprofile文件后，重启R即可。

当前你的文件位置

> .libPaths()
[1] "D:/软件安装/R_libaray"             
[2] "C:/Program Files/R/R-4.1.2/library"

1 导入所需包

#if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") #BiocManager::install("DESeq2") library(DESeq2) library(dplyr)

2 导入Counts数据矩阵

## 导入Counts数据矩阵 countdata <- read.table("Countdata.txt", header = T,row.names = 1) countdata[1:10,1:6] > countdata[1:10,1:6] Control_01 Control_02 Control_03 Treat_01 Treat_02 Treat_03 gene01 11 11 11 10 11 10 gene02 14 14 14 13 14 13 gene03 14 14 14 13 14 13 gene04 7 7 7 6 7 6 gene05 11 11 11 11 11 10 gene06 14 14 14 13 14 12 gene07 13 13 13 12 13 12 gene08 11 11 10 10 11 9 gene09 11 11 11 10 11 9 gene10 10 12 12 12 12 10

3 数据过滤

## 过滤在所有重复样本中小于1的基因 countdata = countdata[rowMeans(countdata) > 1,]

4 数据样本信息

样本注释信息自己在本地后制作后进行导入

coldata <- read.csv("countdata.csv", header = T, row.names = 1) head(coldata) ############# > head(coldata) condition Control_01 CK Control_02 CK Control_03 CK Treat_01 Treat Treat_02 Treat Treat_03 Treat

检查数据Counts文件与coldata数据是否匹配

all(rownames(coldata) %in% colnames(countdata)) all(rownames(coldata) == colnames(countdata))

当返回 TRUE 时，表明两个数匹配。

> all(rownames(coldata) %in% colnames(countdata)) # [1] TRUE > all(rownames(coldata) == colnames(countdata)) [1] TRUE

5 差异分析

## 制作差异矩阵 dds <- DESeqDataSetFromMatrix(countData = countdata,colData = coldata,design = ~ condition) dim(dds) ## 过滤 dds <- dds[rowSums(counts(dds)) > 1,] nrow(dds) ## 差异比较 dep <- DESeq(dds) res <- results(dep) diff = res diff <- na.omit(diff) ## 去除NA dim(diff) write.csv(diff,"all_diff.csv") # 导出所有的差异文件

DESeq输出的差异文件

> head(diff) log2 fold change (MLE): condition Treat vs CK Wald test p-value: condition Treat vs CK DataFrame with 6 rows and 6 columns baseMean log2FoldChange lfcSE stat pvalue padj <numeric> <numeric> <numeric> <numeric> <numeric> <numeric> gene01 10.6594 -0.01071206 0.788745 -0.0135811 0.989164 0.999747 gene02 13.6626 0.00959682 0.713528 0.0134498 0.989269 0.999747 gene03 13.6626 0.00959682 0.713528 0.0134498 0.989269 0.999747 gene05 10.8224 0.03339267 0.783215 0.0426354 0.965992 0.999747 gene06 13.4731 -0.03001613 0.716954 -0.0418662 0.966605 0.999747 gene07 12.6616 0.00389414 0.735478 0.0052947 0.995775 0.999747

筛选差异基因，使用P值和LogFC进行筛选，这里我们就不多做解释，详情可以看我们前面的推文。

# 设置筛选标准 foldChange = 1 padj = 0.05 diffsig <- diff[(diff$pvalue < padj & abs(diff$log2FoldChange) > foldChange),] dim(diffsig) write.csv(diffsig, "All_diffsig.csv")

注：上调和下调的差异基因的筛选，我们在这里不在讲述，请查看我们前面的推文即可

6 差异基因热图

6.1 标准化Count值

我们后面的热图使用Count值进行分析，以及后续的分析也可以使用count值进行分析。我们的Count值是基因的原始表达量，因此需要进行标准化出来，我们这里使用 vst 函数进行标准化处理。

## 标准化(标准化Counts值) vsd <- vst(dds, blind = FALSE) normalizeExp <- assay(vsd)

绘制差异基因热图

### 差异表达热图 df <- read.csv("All_diffsig.csv", header = T) head(df) df02 <- as.character(df$X) ## 标准化(标准化Counts值) vsd <- vst(dds, blind = FALSE) normalizeExp <- assay(vsd) ## 差异基因的Count diff_expr <- normalizeExp[df02,] head(diff_expr) ## 差异热图 library(pheatmap) annotation_col <- data.frame(Group = factor(c(rep("CK",3), rep("Treat",3)))) rownames(annotation_col) <- colnames(diff_expr) pdf("差异基因热图.pdf", height = 8, width = 12) p <- pheatmap(diff_expr, annotation_col = annotation_col, color = colorRampPalette(c("green","black","red"))(50), cluster_cols = F, show_rownames = F, show_colnames = F, scale = "row", ## none, row, column fontsize = 12, fontsize_row = 12, fontsize_col = 6, border = FALSE) print(p) dev.off()

limma包差异分析

1. 数据准备

差异分析是两两数据集间的比较，一个是对 照组样本(CK) ，一个是 处理组样本(Treat) 。

CK01 CK02 CK03 T1 T2 T3 gene0001 56 46 91 36 19 28 gene0002 0 0 5 0 31 8 gene0003 4 0 6 2 2 0 gene0004 257 254 235 241 162 183 gene0005 30 17 22 27 9 75 gene0006 503 565 490 369 426 480 gene0007 201 82 62 296 206 189 gene0008 56 6 12 108 62 0 gene0009 6 10 9 13 8 14 gene0010 19 0 0 27 0 0 gene0011 0 0 0 0 0 0

注：数据样本根据自己实验数据决定

2. 数据分析

2.1 导入所需的包

## 导入R包 library(limma) library(dplyr)

2.2 加载数据

df <- read.table("Counts_data.txt", header = T, sep = "\t", row.names = 1, check.names = F) head(df)

2.3 样本信息注释

list <- c(rep("Treat", 66), rep("CK",148)) %>% factor(., levels = c("CK", "Treat"), ordered = F) ##-------------- > head(list) CK Treat 1 0 1 2 0 1 3 0 1 4 0 1 5 0 1 6 0 1 ............. list <- model.matrix(~factor(list)+0) #把group设置成一个model matrix colnames(list) <- c("CK", "Treat") df.fit <- lmFit(df, list) ## 数据与list进行匹配

2.4 差异分析

df.matrix <- makeContrasts(tumor - normal, levels = list) fit <- contrasts.fit(df.fit, df.matrix) fit <- eBayes(fit) tempOutput <- topTable(fit,n = Inf, adjust = "fdr") > head(tempOutput) logFC AveExpr t P.Value adj.P.Val B gene11796 -1.2620803 5.551402 -8.392278 5.886136e-15 7.083964e-11 23.25729 gene7364 -0.9825962 7.471963 -7.598215 8.629186e-13 5.192613e-09 18.51600 gene11454 -1.3204341 6.238318 -7.368543 3.472935e-12 1.301549e-08 17.19354 gene11689 -1.0769861 4.967290 -7.331950 4.325880e-12 1.301549e-08 16.98503 gene11227 -1.2035217 5.925234 -7.263016 6.531456e-12 1.572122e-08 16.59390 gene3259 -2.3695741 9.467290 -6.962632 3.829470e-11 5.760959e-08 14.91576

到这部分，差异分析就结束了。对！！没错就是这么简简单单...................

后面的部分就是提取差异结果，和绘制火山图。

2.5 导出差异结果

## 导出所有的差异结果 nrDEG = na.omit(tempOutput) ## 去掉数据中有NA的行或列 diffsig <- nrDEG write.csv(diffsig, "all.limmaOut.csv")

2.6 筛选差异基因

差异基因基因的筛选，我们一般是使用P值和LogFC筛选，常用的筛选标准 P<0.05 , |LogFC| > 1 ，这是最常规的筛选标准，如果你的数据差异较大，也可以更改P值和LogFC的大小。

## 我们使用|logFC| > 0.5，padj < 0.05（矫正后P值） foldChange = 0.5 padj = 0.05 ## 筛选出所有差异基因的结果 All_diffSig <- diffsig[(diffsig$adj.P.Val < padj & (diffsig$logFC>foldChange | diffsig$logFC < (-foldChange))),] #--------------------- > dim(All_diffSig) [1] 0 6 ## 我们发现竟然没有差异基因，这是应该我这边的数据是随机的结果，如果你的数据有这样的问题，你需要在仔细检查一下哦。 ## 我们为了下面的操作正常进行，我们选用的P值（未矫正）进行筛选。 All_diffSig <- diffsig[(diffsig$P.Value < padj & (diffsig$logFC>foldChange | diffsig$logFC < (-foldChange))),] write.csv(All_diffSig, "all.diffsig.csv") ##输出差异基因数据集 #----------------- > dim(All_diffSig) [1] 335 6 ## 共有335个差异基因

7）筛选上调和下调的基因

diffup <- All_diffSig[(All_diffSig$P.Value < padj & (All_diffSig$logFC > foldChange)),] write.csv(diffup, "diffup.csv") diffdown <- All_diffSig[(All_diffSig$P.Value < padj & (All_diffSig < -foldChange)),] write.csv(diffdown, "diffdown.csv")

到这部分，我们差异分析就全部结束了。已经拿到差异文件，及上调和下调的基因文件。

3. 绘制火山图

在常规的差异分析之后，我们需要进行差异数据的可视化， 火山图 和 差异基因热图 是最常用的可视化图形。

## 导入R包 library(ggplot2) library(ggrepel) ## 绘制火山图 ## 进行分类别 logFC <- diffsig$logFC deg.padj <- diffsig$P.Value data <- data.frame(logFC = logFC, padj = deg.padj) data$group[(data$padj > 0.05 | data$padj == "NA") | (data$logFC < foldChange) & data$logFC > -foldChange] <- "Not" data$group[(data$padj <= 0.05 & data$logFC > 1)] <- "Up" data$group[(data$padj <= 0.05 & data$logFC < -1)] <- "Down" x_lim <- max(logFC,-logFC) # 开始绘图 pdf('volcano.pdf',width = 7,height = 6.5) ## 输出文件 label = subset(diffsig,P.Value <0.05 & abs(logFC) > 0.5) label1 = rownames(label) colnames(diffsig)[1] = 'log2FC' Significant=ifelse((diffsig$P.Value < 0.05 & abs(diffsig$log2FC)> 0.5), ifelse(diffsig$log2FC > 0.5,"Up","Down"), "Not") ggplot(diffsig, aes(log2FC, -log10(P.Value)))+ geom_point(aes(col=Significant))+ scale_color_manual(values=c("#0072B5","grey","#BC3C28"))+ labs(title = " ")+ geom_vline(xintercept=c(-0.5,0.5), colour="black", linetype="dashed")+ geom_hline(yintercept = -log10(0.05),colour="black", linetype="dashed")+ theme(plot.title = element_text(size = 16, hjust = 0.5, face = "bold"))+ labs(x="log2(FoldChange)",y="-log10(Pvalue)")+ theme(axis.text=element_text(size=13),axis.title=element_text(size=13))+ str(diffsig, max.level = c(-1, 1))+theme_bw() dev.off()

4. 绘制差异基因表达热图

注:本章节，我们是只使用count值进行热图绘制，后续的分子中，我们会对其进行优化

## 导入R包 library(pheatmap) ## 提取差异基因的表达量 DEG_id <- read.csv("all.diffsig.csv", header = T) # 读取差异基因的文件 head(DEG_id) ## 匹配差异基因的表达量 DEG_id <- unique(DEG_id$X) DEG_exp <- df[DEG_id,] hmexp <- na.omit(DEG_exp) ## 样本注释信息 annotation_col <- data.frame(Group = factor(c(rep("Treat", 66), rep("CK",148)))) rownames(annotation_col) <- colnames(hmexp) ## 绘制热图 pdf(file = "heatmap02.pdf", height = 8, width = 12) pheatmap(hmexp, annotation_col = annotation_col, color = colorRampPalette(c("green","black","red"))(50), cluster_cols = F, show_rownames = F, show_colnames = F, scale = "row", ## none, row, column fontsize = 12, fontsize_row = 12, fontsize_col = 6, border = FALSE) print(p) dev.off()

edgeR包差异分析

简介

edgeR 包也是用于做差异分析的包，也算是比较常用的方法。对于我自己来做分析， edgeR 包主要是用来做没有重复的数据差异分析。如果你的数据没有重复，那么建议使用edgeR包来做差异，但是在实验设计是没有设置重复，那么你得想清楚哦！（仅个人观点）

在 edgeR 给出的描述如下下，你可以详细阅读edgeR帮助文档。

1 数据准备

我们使用前面 DESeq2 差异分析中的数据（注：我们有重复），我们的数据仅是练习数据，没有任何参考价值。

Control_01 Control_02 Control_03 Treat_01 Treat_02 Treat_03 gene01 11 11 11 10 11 10 gene02 14 14 14 13 14 13 gene03 14 14 14 13 14 13 gene04 7 7 7 6 7 6 gene05 11 11 11 11 11 10 gene06 14 14 14 13 14 12

2 导入相关包

#install.packages("fdrtool") #install.packages("edgeR") library(edgeR) library("fdrtool")

3 导入数据

rawdata <- read.table("Countdata.txt", header = T, row.names = 1) head(rawdata)

制作样本文件的注释信息

group <- factor(c(rep("CK",3),rep("Treat",3)))

4 数据处理

过滤与标准化

y <- DGEList(counts = rawdata, genes = rownames(rawdata), group = group) keep <- rowSums(cpm(y)>1) >= 1 y <- y[keep,,keep.lib.sizes=FALSE] ##TMM 标准化 y <- calcNormFactors(y) y$samples

5 数据分析

###推测离散度，若样本是人，设置bcv = 0.4，模式生物设置0.1(此处是根据有相关教程进行设置，也可以你根据你的结果自行设置，最终得到你的理想值) bcv <- 0.2 et <- exactTest(y, dispersion=bcv^2) topTags(et) summary(de <- decideTestsDGE(et)) ###图形展示检验结果 png('0h_vs_2h_MAplot.png') detags <- rownames(y)[as.logical(de)]; plotSmear(et, de.tags=detags) abline(h=c(-4, 4), col="blue"); dev.off()

6 矫正P值

#矫正P值 DE <- et$table res <- DE head(DE) res$FDRP <- p.adjust(res$PValue,method = "fdr",n=length(res$PValue)) #res$logP <- -log10(res$FDRP) head(res)

7 筛选出差异基因

## 筛选出差异基因 ## 筛选标准 FDR < 0.05，|logFC| > 1 diffsig <- res[(res$FDRP < 0.05 & abs(res$logFC) > 1),] dim(diffsig) ###---------- > dim(diffsig) [1] 751 5

共得到751个差异基因

筛选上调和下调的基因，并进行标识，方法与前面的一致。

## 新增一列，标显著性 res$Group <- "Not" ###### res$Group[which((res$FDRP < 0.05) & (res$logFC > 1))] = "Up" res$Group[which((res$FDRP < 0.05) & (res$logFC < -1))] = "Down" #### 查看DE数目 table(res$Group) #----------- > table(res$Group) Down Not Up 517 14862 234

热图绘制，必须掌握的技能

pheatmap包绘制热图教程

原文：示例数据获取网址 ： R 数据可视化 01 | 聚类热图

1 数据准备

第一列，gene_id，后面是重复，3个重复，我的共4个处理数据，12组

2导入数据

#文件所在位置 setwd("I:FPKM") library(pheatmap) #BiocManager::install("heatmaps") #如没有按照“pheatmap包请安装” dataset <- read.table('FPKM_聚类图数据.txt',header = TRUE, row.names = 1)

3 构建矩阵

# 构建样本分类数据 sample_calss=c(rep('CS_Leaf',3), rep('CK_Lead',3), rep('CS_Root',3), rep('CK_Root',3)) level = c(1:12) annotation_c <- data.frame(sample_calss, level) rownames(annotation_c) <- colnames(exp_ds) gene_class=c(rep('High',30), rep('Low',30)) sample_type=c(rep('CS',20), rep('CK',20), rep('Immunology',20)) annotation_r <- data.frame(gene_class, gene_type) rownames(annotation_r) <- rownames(exp_ds)

画图

pheatmap(exp_ds, #表达数据 cluster_rows = T,#行聚类 cluster_cols = T,#列聚类 annotation_col =annotation_c, #样本分类数据 annotation_row = annotation_r, annotation_legend=TRUE, # 显示样本分类 show_rownames = T,# 显示行名 scale = "row", #对行标准化 color =colorRampPalette(c("#8854d0", "#ffffff","#fa8231"))(100), # 热图基准颜色 )

结果图形

热图中添加指定的基因

如想获得图例代码和数据，请到公众号“小杜的 生信筆記 ”中回复“20220117”，链接有效期2022.17-2022.2.17（特别申明：该数据仅是演示数据，无任何意义）

作为一个读者、用户者推荐：
基迪奥是第一个我认识并接触过，以及经常使用他们家云平台绘图的生信公司。毕业后也差一点加入他们。我个人的感受是使用他们家的云平台很方便，以及他们的绘图这块还是很强的。看到这篇文章，立即收藏起来。确实是一个很好的教程！！
今天又是来水一期教程，后期将使用他们的码进一步学习。
--- Du

图例： （都是大文章哦）

开始教程

前面的数据导入和绘图都是与我前面的热图的教程是一样（好看的热图绘制教程|重现) ），万变不离其宗。

后面就是重点

已经有这样的味道了！！！

3D主成分分析（PCA) | 图形绘制 | 教程（代码重新）

主成分分析（Principal Component Analysis，PCA)，一种线性降维的方法，它的目标是通过某种线性投影，将高维的数据映射到低维的空间中，并期望在所投影的维度上数据的信息量最大（方差最大），以此使用较少的数据维度，同时保留住较多的原数据点的特性。

不同的分析思路和分析对象，使用PCA结果解释的意义是不同的，主要还是看自己的数据以及自己对其的理解。

二维PCA

3D PCA图形

原文：High-resolution spatiotemporal transcriptomemapping of tomato fruit development and ripening

如果你感兴趣可以去了解一下，非常牛X的，优秀的，nice的一篇文章。数据量非常的大，做的也是非常的细致。

3D PCA 代码

# 导入包 library(rgl) Sys.setenv(LANGUAGE = "en") #显示英文报错信息 options(stringsAsFactors = FALSE) #禁止chr转成factor

导入表达量数据

# 读取表达量数据 df <- read.csv("PCA_inputdata.csv", header = T, row.names = 1) df[1:5,1:5] > df[1:5,1:5] sample001 sample002 sample003 sample004 sample005 AIM2 1.4 0.1 1.1 1.3 -1.2 ANXA11 -0.1 0.8 0.3 0.0 -0.1 APLN 0.3 0.0 0.5 -1.4 0.9 APOA1 0.0 -1.2 0.2 -0.7 0.2 APOA2 -0.1 -0.8 0.4 -0.5 -0. > dim(df) [1] 311 297

PCA分析

## 一个函数搞定 pca <- prcomp(t(df))

绘图

# 01.导入分类数据 lgg_sample <- read.table("group1.txt", header = T)$x ## class(lgg_sample) ## 查看数据类型 gbm_sample <- read.table("group2.txt", header = T)$x # 准备颜色 lgg_color <- rep("yellow", length(lgg_sample)) names(lgg_color) <- lgg_sample gbm_color <- rep("red", length(gbm_sample)) names(gbm_color) <- gbm_sample groups <- c(lgg_color, gbm_color)

3d PCA图

plot3d(pca$x[,1:3], # 取前三个主成分 xlab="Comp.1", ylab="Comp.2", zlab="Comp.3", col=groups, # 按groups填充颜色 type="s", # 画球，'p' for points, 's' for spheres, 'l' for lines, 'h' for line segments size=1, #球的大小 lwd=2, box=T) rgl.snapshot("PCA01.png")

到此步，你可以拖动你的鼠标来确定图形的方向和大小，很有意思。这是一个动态的图形。

但是比较遗憾，此图形只能输出 .png 格式的图形。如你其他的绘制方法，你也可以一起来分享一下哦！

## 五颜六色的图形，你也喜欢 plot3d(pca$x[,1:3], col=c("red","gray0"),# "blue","cyan","darkblue","green","darkgreen","lightpink"), size = 10, xlab="PC1",ylab="PC2",zlab="PC3")

结束语：

如果你看到这里，那么你是最棒的！我在前面的教程中也提及，我们最终还是得回归到绘图的本质。

那什么是“绘图的本质呢”，个人观点：本质就是你开始绘制图形的初心，我们一切都是为了用来解释说明数据的可用性，用图形的形式直观的表达出数据。数据是不会直接告诉你它所想要表达的意思，那只有看 分析者 使用图形来阐述。我们要理解这个图形的最基础的东西，如它的原理之类的，我们可以无需具体理解这个图形绘制的每一个步骤，因为你不是开发软件或开发数据包的。

但是，我们还是需要理解它运作的原理是什么？这部分，也许是可能被忽略的。我们记住： 万变不离其宗 ，越是最基础的东西，越容易被忽略，但也是最重要的。基本，所有 形 变，都是在 质 的基础上。

我们在追求高质量图形的同时，也需要追求高质量的意义！（不知你是否理解）

桑基图（SanKey） | 教程（代码重新）

## 设置路劲 setwd("D:\\小杜的生信筆記\\桑基图") ## 导入R包 #BiocManager::install("ggalluvial") library("ggplot2") library("ggalluvial")

导入数据：

df <- read.table("SangKey_input.txt", header = T, row.names = 1, sep = '\t') head(df)

数据格式:

> head(df) gene lncRNA miRNA 1 gene4 lncRNA1 miRNA2 2 gene4 lncRNA3 miRNA3 3 gene3 lncRNA5 miRNA3 4 gene2 lncRNA5 miRNA2 5 gene8 lncRNA4 miRNA1 6 gene8 lncRNA2 miRNA3

开始绘图:

#定义颜色 mycol <- rep(c("#223D6C","#D20A13","#FFD121","#088247","#11AA4D","#58CDD9", "#7A142C","#5D90BA","#029149","#431A3D","#91612D","#6E568C", "#E0367A","#D8D155","#64495D","#7CC767"),2)

格式转换:

#格式转换 UCB_lodes <- to_lodes_form(df[,1:ncol(df)], axes = 1:ncol(df), id = "Cohort")

> head(UCB_lodes) Cohort x stratum 1 1 gene gene4 2 2 gene gene4 3 3 gene gene3 4 4 gene gene2 5 5 gene gene8 6 6 gene gene8

--

> tail(UCB_lodes) Cohort x stratum 295 95 miRNA miRNA4 296 96 miRNA miRNA3 297 97 miRNA miRNA3 298 98 miRNA miRNA2 299 99 miRNA miRNA3 300 100 miRNA miRNA2

画图：

pdf("Sankey.pdf",width = 6, height = 6) ggplot(UCB_lodes, aes(x = x, stratum = stratum, alluvium = Cohort, fill = stratum, label = stratum)) + scale_x_discrete(expand = c(0, 0)) + geom_flow(width = 1/8) + #线跟方块间空隙的宽窄 geom_stratum(alpha = .9,width = 1/10) + #方块的透明度、宽度 geom_text(stat = "stratum", size = 3,color="black") + #文字大小、颜色 #不喜欢默认的配色方案，用前面自己写的配色方案 scale_fill_manual(values = mycol) + xlab("") + ylab("") + theme_bw() + #去除背景色 theme(panel.grid =element_blank()) + #去除网格线 theme(panel.border = element_blank()) + #去除外层边框 theme(axis.line = element_blank(),axis.ticks = element_blank(),axis.text = element_blank()) + #去掉坐标轴 ggtitle("")+ guides(fill = FALSE) dev.off()

R语言 | 堆积图绘制教程 | 纯代码分享

PS：如果你需要本教程的练习代码和数据，可以在 公众号 （点击）回复“20220121”即可获得。

绘图

## date: 2022.0121 ## author:小杜的生信筆記 # 堆积图 rm(list = ls()) library("ggplot2")

导入数据

setwd("D:\\小杜的生信筆記\\堆积图") df <- read.table("input_data.txt", header = T) head(df) > head(df) number Type Treatment 1 142 mRNAs Treat01 2 246 miRNAs Treat01 3 464 mRNAs Treat02 4 235 miRNAs Treat02 5 363 mRNAs Treat03 6 536 miRNAs Treat03

作图

基础图形

ggplot(df, aes(y= Treatment, x = number, fill = Type))+ ## 使用ggplot2语法 geom_bar(stat = "identity", position = "stack", color = "black",width = 0.8) ## 添加柱子

对图形美化

这个颜色的配置是我比较喜欢的颜色，如果你喜欢其他颜色，也可以自己进行配置

palette = Set1至Set**,Greens,`Accent,自己根据其选择吧。

ggplot(df, aes(y= Treatment, x = number, fill = Type))+ ## 使用ggplot2语法 geom_bar(stat = "identity", position = "stack", color = "black",width = 0.8)+ ## 添加柱子 theme_bw()+ scale_fill_brewer(palette = "Accent") ## 添加柱子颜色，个人非常喜欢这个配置的颜色

去掉方框

ggplot(df, aes(y= Treatment, x = number, fill = Type))+ ## 使用ggplot2语法 geom_bar(stat = "identity", position = "stack", color = "black",width = 0.8)+ ## 添加柱子 theme_bw()+ scale_fill_brewer(palette = "Accent")+ ## 添加柱子颜色，个人非常喜欢这个配置的颜色 theme_classic()+ theme(legend.title=element_blank())

调整柱子距离X轴和Y轴的距离

## 添加下面的吗即可 scale_x_continuous(expand = c(0.05,0))+ ## 限制x轴坐标的距离y轴的距离 theme(text = element_text(size=15), ## y轴字体大小 legend.position = c(0.8,0.4))+ #theme(legend.position='none') 不要主题##； theme(legend.position = c(0.8,0.4)) 调节位置

更改X轴和y轴的字体，颜色，大小

theme(axis.text.x = element_text(face = "bold", color = "black", size = 12), axis.text = element_text(colour = "black", size = 15) face = : 字体 color = : 颜色 size = : 字体大小

更改X轴和y轴的坐标

ylab("")+xlab("Number of **")+ ## 添加x轴坐标

更改或添加标题

labs(title = "input your title")+ ## 添加标题 theme(plot.title = element_text(hjust = 0.5, size = 12)) ##标题的调节

2022年8月22日

R语言绘制截断图小例子 | 直播笔记
截断图学习笔记

学习网址： Set Axis Break for ggplot2 (r-project.org)

导入相关的R包

## 导入R包 library(ggplot2) library(ggbreak) library(patchwork) # 下载你所需要的R包 # install.packages("") # BiocManager::install("") ###

导入数据

setwd("D:\\小杜的生信筆記\\截断图绘制") df <- read.table("inut.data.txt", header = T) head(df) 本教程使用的随机数据 d <- data.frame(x = 1:20, y = c(rnorm(5) + 4, rnorm(5) + 20, rnorm(5) + 5, rnorm(5) + 22) head(d) > head(d) x y 1 1 3.423056 2 2 5.148453 3 3 4.596218 4 4 4.641941 5 5 4.902195 6 6 20.323737

基础图形的绘制

p1 <- ggplot(d, aes(y, x)) + geom_col(orientation="y")

图形的调整，截断

## 在图形中随意位置加入你想要加入的字符，或标记 d2 <- data.frame(x = c(2, 18), y = c(7, 26), label = c("hello", "world")) p2 <- p1 + scale_x_break(c(7,17))+ ## 截断x轴 # scale_y_break(c("输入你的坐标")) ## 截断y轴 geom_text(aes(y, x, label = label), data = d2, hjust = 1, colour = 'firebrick')+ labs(y = "Time", x = "Expression level")+ theme_classic() ## 设置主题 # scale_fill_brewer(palette = "Accent") ## 设置颜色 ### 第二次截断 p2 + scale_x_break(c(18,21)) ### 第三次截断，第四次截断.......................... ## 建议：截断图，最好只截断一次

注意：建议：截断图，最好只截断一次

03. 竖着的柱状图-截断例子

ggplot(d , aes(x,y))+ geom_col()+ #基础图形绘制 scale_y_break(c(7,12), scales = 1.5)+ scale_y_break(c(18,21),scales = 2)+ # scale_y_reverse() labs(fill = "Time(h)",x = "", y = "Expression level")+ theme(text = element_text(size=12))+ ## 更改横纵坐标轴中的字体颜色、大小 theme(axis.text.x = element_text(color = "black",size = 10), axis.text.y = element_text(color = "black",size = 10))

04. 截断图中级篇

# 01. 准备 ***** # 02. 图形绘制 # 03. 图形美化 # 04. 图中字体大小更改美化

数据准备

d <- data.frame( x = 1:20, y = c(rnorm(5) + 4, rnorm(5) + 20, rnorm(5) + 5, rnorm(5) + 22), group = c(rep("A", 10), rep("B", 10)), face=c(rep("C", 5), rep("D", 5), rep("E", 5), rep("F", 5)) ## 查看数据 head(d) dim(d) > head(d) x y group face 1 1 3.809619 A C 2 2 1.729871 A C 3 3 3.031275 A C 4 4 5.819765 A C 5 5 3.860474 A C 6 6 20.886914 A D > dim(d) [1] 20 4

基础图形的绘制

p <- ggplot(d, aes(x = x, y = y))+ geom_col(orientation = "x")+ ##参数：orientation = "":t图形需要从哪一个参数进行转变 scale_y_reverse()+ facet_wrap(group~., scales="free_y", strip.position="right", nrow=2)+ coord_flip() ## 顺时针旋转90

截断x轴中的图形

p2 <- p+scale_y_break(c(7,17),scales = "free")

给图形添加颜色

p2 + aes(fill =group)+ theme(legend.position = "bottom")+ theme_classic()+ scale_fill_brewer(palette = "Accent")

require(ggplot2) library(ggbreak) set.seed(2019-01-19) d <- data.frame( x = 1:20, y = c(rnorm(5) + 4, rnorm(5) + 20, rnorm(5) + 5, rnorm(5) + 22), group = c(rep("A", 10), rep("B", 10)) p <- ggplot(d, aes(x=x, y=y)) + scale_y_reverse() + scale_x_reverse() + geom_col(aes(fill=group)) + scale_fill_manual(values=c("#00AED7", "#009E73")) + facet_wrap( group~., scales="free_y", strip.position="right", nrow=2 coord_flip() scale_y_break(c(7, 10), scales=0.5, ticklabels=c(10, 11.5, 13)) + scale_y_break(c(13, 17), scales=0.5, ticklabels=c(17, 18, 19)) + scale_y_break(c(19,21), scales=1, ticklabels=c(21, 22, 23))

加权基因共表达网络分析 (WGCNA）

一、前言

1.1 概念

加权基因共表达网络分析 (WGCNA, Weighted correlation network analysis)是用来描述不同样品之间基因关联模式的系统生物学方法，可以用来鉴定高度 协同变化 的基因集, 并根据基因集的内连性和基因集与表型之间的关联鉴定候补生物标记基因或治疗靶点。

该分析方法旨在寻找协同表达的基因模块(module)，并探索基因网络与关注的表型之间的关联关系，以及网络中的核心基因。

适用于复杂的数据模式，推荐5组(或者15个样品)以上的数据。一般可应用的研究方向有：不同器官或组织类型发育调控、同一组织不同发育调控、非生物胁迫不同时间点应答、病原菌侵染后不同时间点应答

2.2 原理

从方法上来讲，WGCNA分为 表达量聚类分析和表型关联 两部分，主要包括基因之间相关系数计算、基因模块的确定、共表达网络、模块与性状关联四个步骤。

第一步计算任意两个基因之间的相关系数（Person Coefficient）。为了衡量两个基因是否具有相似表达模式，一般需要设置阈值来筛选，高于阈值的则认为是相似的。但是这样如果将阈值设为0.8，那么很难说明0.8和0.79两个是有显著差别的。因此， WGCNA分析时采用相关系数加权值，即对基因相关系数取N次幂 ，使得网络中的基因之间的连接服从 无尺度网络分布(scale-freenetworks) ，这种算法更具生物学意义。

第二步通过基因之间的相关系数构建分层聚类树， 聚类树的不同分支代表不同的基因模块 ，不同颜色代表不同的模块。基于基因的加权相关系数，将基因按照表达模式进行分类，将模式相似的基因归为一个模块。这样就可以 将几万个基因通过基因表达模式被分成了几十个模块 ，是一个提取归纳信息的过程。

共表达网络 ：定义为加权基因网络。点代表基因，边代表基因表达相关性。加权是指对 相关性值进行冥次运算 (冥次的值也就是 软阈值 (power, pickSoftThreshold这个函数所做的就是确定合适的power))。无向网络的边属性计算方式为 abs(cor(genex, geney)) ^ power；有向网络的边属性计算方式为 (1+cor(genex, geney)/2) ^ power; sign hybrid的边属性计算方式为cor(genex, geney)^power if cor>0 else 0

Module(模块) ：高度內连的基因集。在无向网络中，模块内是高度相关的基因。在有向网络中，模块内是高度 正相关 的基因。把基因聚类成模块后，可以对每个模块进行三个层次的分析：1. 功能富集分析查看其功能特征是否与研究目的相符；2. 模块与性状进行关联分析，找出与关注性状相关度最高的模块；3. 模块与样本进行关联分析，找到样品特异高表达的模块。

Eigengene （eigen +‎ gene） ：基因和样本构成的矩阵， https:// en.wiktionary.org/wiki/ eigengene

Connectivity (连接度) ：类似于网络中 "度" (degree)的概念。每个基因的连接度是与其相连的基因的边属性之和。

Module eigengene E : 给定模型的第一主成分，代表 整个模型的基因表达谱 。这个是个很巧妙的梳理，我们之前讲过 PCA分析的降维作用，之前主要是拿来做可视化，现在用到这个地方，很好的用一个向量代替了一个矩阵，方便后期计算。(降维除了PCA，还可以看看 tSNE )

Intramodular connectivity : 给定基因与给定模型内其他基因的关联度，判断基因所属关系。

Module membership : 给定基因表达谱与给定模型的eigengene的相关性。

Hub gene: 关键基因 (连接度最多或连接多个模块的基因)。

Adjacency matrix (邻接矩阵) ：基因和基因之间的加权相关性值构成的矩阵。

TOM (Topological overlap matrix) ：把邻接矩阵转换为拓扑重叠矩阵，以降低噪音和假相关，获得的新距离矩阵，这个信息可拿来构建网络或绘制TOM图。

以上解释来自 ： WGCNA分析，简单全面的最新教程； https://www. jianshu.com/p/e9cc3f434 41d

分析案例

Integrative analysis of metabolome and transcriptome reveals the mechanism of color formation in pepper fruit (Capsicum annuum L.), 2020

Genome-Wide Identification and Characterization of Long Non-Coding RNAs in Peanut, 2019

--

NCBI搜索结果： WGCNA - Search Results - PubMed (nih.gov)

WGCNA分析，如果你的数据样本较大，建议在服务器中进行运行！！

二、WGCNA 代码教程

2.1 加载所需的R包

加载WGCNAR包

#install.packages("WGCNA") #BiocManager::install('WGCNA') library(WGCNA)

进行基础设置，（此步骤可以省略）

options(stringsAsFactors = FALSE) enableWGCNAThreads() ## 打开多线程

2.2 加载数据矩阵

WGCNA.fpkm = read.table("WGCNA.inputdata.txt",header=T, comment.char = "", check.names=F)

数据矩阵格式

> WGCNA.fpkm[1:10,1:10] sample MT_0_1 MT_0_2 MT_0_3 MT_3_1 MT_3_2 MT_3_3 MT_6_1 MT_6_2 MT_6_3 1 gene10 0 0.000000 0.000000 0.000000 0.000000 0.000000 0.225467 0.000000 0.116728 2 gene12 0 0.000000 0.000000 0.000000 0.140796 0.000000 0.051667 0.076227 0.109181 3 gene13 0 0.000000 0.000000 0.000000 0.000000 0.000000 0.090672 0.000000 0.000000 4 gene14 0 0.000000 0.247270 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 5 gene16 0 0.000000 0.213587 0.051469 0.433966 0.786561 0.000000 0.246451 0.972179 6 gene19 0 0.113876 0.060553 0.000000 0.331847 0.000000 0.099846 0.224528 0.191205 7 gene20 0 0.000000 0.000000 0.183355 0.531899 0.383033 0.000000 2.460091 2.441143 8 gene21 0 0.470667 0.000000 1.016470 0.583894 0.000000 1.074058 2.020643 1.007870 9 gene24 0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 10 gene25 0 0.054502 0.026962 0.026076 0.025938 0.000000 0.028589 0.023683 0.000000

查看数据名称

names(WGCNA.fpkm) [1] "sample" "MT_0_1" "MT_0_2" "MT_0_3" "MT_3_1" "MT_3_2" "MT_3_3" "MT_6_1" "MT_6_2" "MT_6_3" "MT_12_1" "MT_12_2" [13] "MT_12_3" "HG_0_1" "HG_0_2" "HG_0_3" "HG_3_1" "HG_3_2" "HG_3_3" "HG_6_1" "HG_6_2" "HG_6_3" "HG_12_1" "HG_12_2" [25] "HG_12_3"

对数据进行提取

names(WGCNA.fpkm) datExpr0 = as.data.frame(t(WGCNA.fpkm[,-1])) names(datExpr0) = WGCNA.fpkm$sample; rownames(datExpr0) = names(WGCNA.fpkm[,-1]) check missing value gsg = goodSamplesGenes(datExpr0, verbose = 3) gsg$allOK if (!gsg$allOK) # Optionally, print the gene and sample names that were removed: if (sum(!gsg$goodGenes)>0) printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ", "))) if (sum(!gsg$goodSamples)>0) printFlush(paste("Removing samples:", paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", "))) # Remove the offending genes and samples from the data: datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes] }

filter

meanFPKM= 0.5 ####--过滤标准，可以修改 n=nrow(datExpr0) datExpr0[n+1,]=apply(datExpr0[c(1:nrow(datExpr0)),],2,mean) datExpr0=datExpr0[1:n,datExpr0[n+1,] > meanFPKM] # for meanFpkm in row n+1 and it must be above what you set--select meanFpkm>opt$meanFpkm(by rp) filtered_fpkm=t(datExpr0) filtered_fpkm=data.frame(rownames(filtered_fpkm),filtered_fpkm) names(filtered_fpkm)[1]="sample" head(filtered_fpkm) write.table(filtered_fpkm, file="mRNA.symbol.uniq.filter.txt", row.names=F, col.names=T,quote=FALSE,sep="\t")

2.3 对input data进行聚类

#############Sample cluster########### sampleTree = hclust(dist(datExpr0), method = "average") pdf(file = "1.sampleClustering.pdf", width = 15, height = 8) par(cex = 0.6) par(mar = c(0,6,6,0)) plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 2, cex.axis = 1.5, cex.main = 2) ### Plot a line to show the cut #abline(h = 180, col = "red")##剪切高度不确定，故无红线 dev.off()

如果数据中有个别数据与其他数据有较大的差异性，为了让后续分析更精准，可以将其进行适当的删除即可。如上图所示。

2.4 加载性状数据

allTraits = read.table("type.txt",row.names=1,header=T,comment.char = "",check.names=F) dim(allTraits) names(allTraits) head(allTraits) > head(allTraits) Cold_3h Cold_6h Cold_12 MT_0_1 0 0 0 MT_0_2 0 0 0 MT_0_3 0 0 0 MT_3_1 3 0 0 MT_3_2 3 0 0 MT_3_3 3 0 0

表达量与性状数据进行匹配

fpkmSamples = rownames(datExpr0) traitSamples =rownames(allTraits) traitRows = match(fpkmSamples, traitSamples) datTraits = allTraits[traitRows,] rownames(datTraits) collectGarbage() pdf(file="2.Sample_dendrogram_and_trait_heatmap.pdf",width=20,height=12) plotDendroAndColors(sampleTree2, traitColors, groupLabels = names(datTraits), main = "Sample dendrogram and trait heatmap",cex.colorLabels = 1.5, cex.dendroLabels = 1, cex.rowText = 2) dev.off()

2.5 选择最优的sftpower值

# 选择最优的sftpower值 # softPower = sft$powerEstimate softPower = 9

经验power (无满足条件的power时选用)

# 无向网络在power小于15或有向网络power小于30内，没有一个power值可以使 # 无标度网络图谱结构R^2达到0.8，平均连接度较高如在100以上，可能是由于 # 部分样品与其他样品差别太大。这可能由批次效应、样品异质性或实验条件对 # 表达影响太大等造成。可以通过绘制样品聚类查看分组信息和有无异常样品。 # 如果这确实是由有意义的生物变化引起的，也可以使用下面的经验power值。 if (is.na(power)){ power = ifelse(nSamples<20, ifelse(type == "unsigned", 9, 18), ifelse(nSamples<30, ifelse(type == "unsigned", 8, 16), ifelse(nSamples<40, ifelse(type == "unsigned", 7, 14), ifelse(type == "unsigned", 6, 12)) 2.6 网络构建 adjacency = adjacency(datExpr0, power = softPower) TOM = TOMsimilarity(adjacency); dissTOM = 1-TOM # Call the hierarchical clustering function geneTree = hclust(as.dist(dissTOM), method = "average"); # Plot the resulting clustering tree (dendrogram) #sizeGrWindow(12,9) pdf(file="4_Gene clustering on TOM-based dissimilarity.pdf",width=24,height=18) plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity", labels = FALSE, hang = 0.04) dev.off()

层级聚类树展示各个模块

# 选择模块的数量，可选大小 minModuleSize = 30 # Module identification using dynamic tree cut: dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, deepSplit = 2, pamRespectsDendro = FALSE, minClusterSize = minModuleSize); table(dynamicMods) # Convert numeric lables into colors dynamicColors = labels2colors(dynamicMods) table(dynamicColors) # Plot the dendrogram and colors underneath #sizeGrWindow(8,6) pdf(file="5_Dynamic Tree Cut.pdf",width=8,height=6) plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut", dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05, main = "Gene dendrogram and module colors") dev.off()

模块间进行矫正合并

计算模块合并的高度

# Calculate eigengenes MEList = moduleEigengenes(datExpr0, colors = dynamicColors) MEs = MEList$eigengenes # Calculate dissimilarity of module eigengenes MEDiss = 1-cor(MEs); # Cluster module eigengenes METree = hclust(as.dist(MEDiss), method = "average") # Plot the result #sizeGrWindow(7, 6) pdf(file="6_Clustering of module eigengenes.pdf",width=7,height=6) plot(METree, main = "Clustering of module eigengenes", xlab = "", sub = "") MEDissThres = 0.4######剪切高度可修改 # Plot the cut line into the dendrogram abline(h=MEDissThres, col = "red") dev.off()

重新绘制合并后的模块层次聚类图

# Call an automatic merging function merge = mergeCloseModules(datExpr0, dynamicColors, cutHeight = MEDissThres, verbose = 3) # The merged module colors mergedColors = merge$colors # Eigengenes of the new merged modules: mergedMEs = merge$newMEs table(mergedColors) #sizeGrWindow(12, 9) pdf(file="7_merged dynamic.pdf", width = 9, height = 6) plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors), c("Dynamic Tree Cut", "Merged dynamic"), dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05) dev.off()

2.7 各个模块与数量性状间的相关性计算

## 模块与数量性状间的相关性 moduleColors = mergedColors # Construct numerical labels corresponding to the colors colorOrder = c("grey", standardColors(50)) moduleLabels = match(moduleColors, colorOrder)-1 MEs = mergedMEs nGenes = ncol(datExpr0) nSamples = nrow(datExpr0) moduleTraitCor = cor(MEs, datTraits, use = "p") moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples) pdf(file="8_Module-trait relationships.pdf",width=10,height=10) textMatrix = paste(signif(moduleTraitCor, 2), "\n(", signif(moduleTraitPvalue, 1), ")", sep = "") dim(textMatrix) = dim(moduleTraitCor) par(mar = c(6, 8.5, 3, 3)) labeledHeatmap(Matrix = moduleTraitCor, xLabels = names(datTraits), yLabels = names(MEs), ySymbols = names(MEs), colorLabels = FALSE, colors = greenWhiteRed(50), textMatrix = textMatrix, setStdMargins = FALSE, cex.text = 0.5, zlim = c(-1,1), main = paste("Module-trait relationships")) dev.off()

2.8 可视化基因网络 (TOM plot )

数据处理

nGenes = ncol(datExpr0) nSamples = nrow(datExpr0) nSelect = 400 # For reproducibility, we set the random seed set.seed(10) select = sample(nGenes, size = nSelect) selectTOM = dissTOM[select, select] selectTree = hclust(as.dist(selectTOM), method = "average") selectColors = moduleColors[select] plotDiss = selectTOM^7 diag(plotDiss) = NA library("gplots") pdf(file="Network heatmap plot_selected genes.pdf",width=9, height=9) mycol = colorpanel(250,'red','orange','lemonchiffon') TOMplot(plotDiss, selectTree, selectColors, col=mycol ,main = "Network heatmap plot, selected genes") dev.off()

2.9 绘制模块之间相关性图

pdf(file="Eigengene dendrogram and Eigengene adjacency heatmap.pdf", width=5, height=7.5) par(cex = 0.9) plotEigengeneNetworks(MEs, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle= 90) dev.off()

pdf(file="Eigengene dendrogram_2.pdf",width=6, height=6) par(cex = 1.0) plotEigengeneNetworks(MEs, "Eigengene dendrogram", marDendro = c(0,4,2,0), plotHeatmaps = FALSE) dev.off()

pdf(file="Eigengene adjacency heatmap_2.pdf",width=6, height=6) par(cex = 1.0) plotEigengeneNetworks(MEs, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2), plotDendrograms = FALSE, xLabelsAngle = 90) dev.off()

参考：

WGCNA分析，简单全面的最新教程

一文看懂WGCNA 分析(2019更新版)

minimap2序列比对(pfar包可视化)

前言

在前面的某一天，我无意间看到关于使用pfar可视化比对后的教程，感觉图还是比较好看的，因此，自己使用番茄基因组进行测试。理想是丰富的，但是现实是残酷的。两个番茄基因组的比较后最终的到22G的比对文件。emm。我看到后就感觉是没有戏的，但还是依旧进行后续的操作，直到今天，依旧是没有最终的结果。后面我又试问自己，为什么要使用基因组的序列进行比较呢？原计划是想看两个基因组间的差异，但............

自己根据前面博主们的教程，重复代码。其次，也找相关的资料进行学习，比如minimap2的比对软件，‘pafr’R包的使用。这些软件在前面我都是不知道的，比对软件，我试用较多的是samtools软件，minimap2与samtools比对结果是一样一样的。对于个人来说，目前学习生物信息学还是比较方便的，有很多的博主都会分享很多好的教程，非常的方便。但，也是比较困难的。

为什么这样说呢，因为我们大多数人都是想把码直接粘贴复制过来，这样确实非常的方便。但，反问自己一句，粘贴复制以后你还记得吗？不会记得多少吧。我一直建议我们初学小白，要自己动手打码，一方面是让自己记码熟悉码，另一方面是为后期自己遇到报错问题做铺垫。对于不是学生信的小白，只有不停的练习，你才有进步哦！！骚年加油哦！！

--du

1 比对序列准备

##基因组大小 > du -sh * 759M Sl.fa 958M Spenn.fasta

各个基因组的序列数量：

> grep ">" Sl.fa >SL4.0ch00 >SL4.0ch01 >SL4.0ch02 >SL4.0ch03 >SL4.0ch04 >SL4.0ch05 >SL4.0ch06 >SL4.0ch07 >SL4.0ch08 >SL4.0ch09 >SL4.0ch10 >SL4.0ch11 >SL4.0ch12 >Spenn-ch00 >Spenn-ch01 >Spenn-ch02 >Spenn-ch03 >Spenn-ch04 >Spenn-ch05 >Spenn-ch06 >Spenn-ch07 >Spenn-ch08 >Spenn-ch09 >Spenn-ch10 >Spenn-ch11 >Spenn-ch12

2. 使用minimap2进行比对

Minimap2是李恒大牛在2018年开发的针对于三代测序数据进行比对的工具。

优点： minimap2的优势是速度快。

缺点： 耗费内存。

github_LH： https:// github.com/lh3/minimap

## 安装 conda search minimap2 conda install -y minimap2

minimap2 比对

minimap2 -ax asm5 Sl.fna Spenn.fna > SL_SPenn.paf

3. minimap2介绍

网址：nanopore测序技术专题（十八）：minimap2比对；

https:// zhuanlan.zhihu.com/p/92 701077?from voters page=true

序列比对

二代测序的核心步奏就是将测序得到的数据重新比对到基因组上，这个基因组可以是通过测序数据拼接得到的，也可以是近源参考序列。这个过程叫做短序列比对，也叫做“reads mapping”。由于测序数据本身就是从基因组上测到的，那么很显然还能定位回去。比对之后，相当于将测序数据与基因组数据组合在一起，那么整个数据分析所有的内容都打包到这个比对结果中了，后续所有的分析都将从这个文件中获取，例如可以知道每条序列的比对情况，也可以得到基因组上每个位点的细节信息，例如每个位点被测序了多少次，所有碱基是相同还是不同。

minimap2简介

二代测序时代利用bwa软件完成这个过程，而相比于illumina测序，nanopore的测序读长更长，测序错误更多，因此必须采取新的比对策略。这个时候，我们的大神级人物，bwa软件的作者李恒，紧跟技术发展趋势，及时开发除了适用于三代测序数据（pacbio，nanopore）的比对工具minimap。现在已更新到minimap2版本，minimap2与bwa必读策略不同，要适应长读长，高测序错误的数据。如果对比对原理算法感兴趣，可以查阅minimap的文献，对于很大一部分用户，了解如何使用工具即可。

minimap2安装

没错，利用minimap2将nanopore测序数据比对到参考序列上就是整个nanopore数据分析的核心，因为序列拼接当中要用到minimap2的比对，如果查看一些拼接软件的源代码，就会发现很多软件都要调用minimap（或者minimap2）比对。而对于变异检测，也是先利用将测序数据（fastq格式）与参考序列（fasta格式）进行比对，得到比对结果（bam格式）。然后在利用一些变异检测软件得到潜在变异结果（vcf格式）。软件的安装也非常容易，可以直接使用bioconda安装，也可以执行变异，推荐大家尝试下执行编译安装。

curl -L https://github.com/lh3/minimap2/releases/download/v2.17/minimap2-2.17_x64-linux.tar.bz2 tar -jxvf minimap2-2.17_x64-linux.tar.bz2 cd minimap2-2.17_x64-linux/ make

软件使用案例

minimap2有多种比对功能，处理支持nanopre数据之外，还支持pacbio数据。比对模式可以是reads与reads之前比对，reads与基因组比对，基因组与基因组比对，以及短序列与基因组的比对，不同的比对有不同的作用，千万不能设置错误了。下面我们拿具体案例演示一下。 minimap2最常用的功能就是将测序数据比对到基因组上，这个过程与bwa比对类似，需要先建立索引，然后比对，最终得到sam格式的比对结果，如果熟悉bwa比对，那么这个操作就非常容易。 minimap2的比对输入文件为测序的reads，一般为fastq或者fasta格式，参考基因组，一般为fasta格式。minimap2可以输出paf格式以及sam格式，默认为paf格式。第一步：建立索引

minimap2 -d co92.min co92.fna

第二步：比对

minimap2 -ax map-ont co92.min ../4.filter/clean.filtlong.fq.gz >s1037.sam

常用选项参数

主要分成五大类，索引（Indexing），回帖（Mapping），比对(Alignment),输入/输出（Input/Output）,预设值（Preset）。 -x ：非常中要的一个选项，软件预测的一些值，针对不同的数据选择不同的值 map-pb/map-ont: pb或者ont数据与参考序列比对； ava-pb/ava-ont: 寻找pd数据或者ont数据之间的overlap关系； asm5/asm10/asm20: 拼接结果与参考序列进行比对，适合~0.1/1/5% 序列分歧度； splice: 长reads的切割比对 sr: 短reads比对 -d :创建索引文件名 -a ：指定输出格式为sa格式，默认为PAF -Q ：sam文件中不输出碱基质量 -R ：reads Group信息，与bwa比对中的-R一致 -t：线程数，默认为3

自己本人数据敲击结果

•结果文件非常的大，有22G，可想而知这个需要本地的电脑性能多大，因此不建议使用本地分析。

4 R语言可视化

比对结果，使用R包pafr进行可视化

注意： 我们这里由于使用两个参考基因组进行比对，数据太大，最终还是不能进行可视化，还是按教程中所说的，使用pafr包中的数据进行进行可视化。

温馨提示 ：自己本地电脑不要轻易尝试大数据量的分析（尤其是机械硬盘的童鞋，你的硬盘会直接卡死的，或直接报废）。

OK！

————————————

导入相关的R包

## date: 2022.3.19 ## author:小杜的生信筆記 # PAFR可视化多序列比对结果 rm(list =ls()) ## set pwd setwd("E:\\小杜的生信筆記\\minimap2序列比对") BiocManager::install("pafr") library(pafr) library(ggplot2) library(tidyverse) library(patchwork)

导入比对结果

df <- read_paf("D:\\Software Install\\R-4.1\\library\\pafr\\extdata\\fungi.paf") head(df) ## 查看结果 df %>% as.data.frame() %>% head()

>df %>%as.data.frame() %>%head() qname qlen qstart qend strand tname tlen tstart tend nmatch alen mapq NM ms AS nn tp cm 1 Q_chr1 627342040169594369152 - T_chr2 589359416366961988044349909353117 60 3208335984337271 0 P 33702 2 Q_chr1 627342056848416086379 - T_chr5 4620715 174028 567211379711414343 6034632302264333848 0 P 36435 3 Q_chr1 627342033491003725888 - T_chr2 589359422803152677081370043402876 6032833295596326118 0 P 35266 4 Q_chr1 627342049653585217172 - T_chr5 462071510989861354853250581256570 60 5989231517235520 0 P 23778 5 Q_chr1 627342052193275466197 - T_chr5 4620715 8506861095265237454253570 6016116212804213220 0 P 22842 6 Q_chr1 627342021397622372339 - T_chr2 589359435001943714260210968235303 6024335165723180683 0 P 20290 s1 s2 dv 134049528250.0023 236863729250.0048 335950634650.0032 4240442 560.0039 523097712540.0039 620475832450.0067 1 11073M1I168M12D497M2D1652M6D712M1I48M3I29M7D6967M1I3094M11D121M1D1007M2I24M10D2221M1D1064M2I451M1I594M13D1001M2I38M1I253M10D529M9D2865M13D1087M9D2524M1I1933M7I39M3I84M1I23M3D513M8D148M6I10M1I190M1D2101M1I30M1I97M1I387M15D104M1I58M2D678M2D6M2D629M2D4389M2D2243M2D25M4I4572M18I31M1I59M4D52M2D163M8I479M3D27M3D18M1I445M5I1938M1D30M14I624M2I1642M11D30M2I5719M934I4841M2D11M5I168M4D442M35I194M1D150M2I28M1I244M1D785M1D1227M6D177M1I1280M2D342M1I23M2D100M12D209M18D72M24D1743M3D501M3D2179M18D224M9D210M2D747M4D2218M9I1779M13D272M1I1521M15I1783M1D515M5D604M5D7992M7D1117M7D3938M4D2262M229I10746M3D2071M3D510M6D90M1D366M1I1263M2I4658M2D6941M1D713M1I426M10I484M3D3626M27I1280M5I336M6D2409M11D2047M1D298M8D2194M1D3853M1I3302M1I3578M4D338M1I595M1I652M1I35M10I354M14D410M2I3417M4I283M1I2378M1I712M1I2451M1D352M2D428M14D45M9D694M1I65M1D1706M8D140M1D89M1I535M1I15M11D13M1D380M6D2986M21I3988M1I183M10D3918M1I3235M14I609M1I301M4D30M1D141M2I1398M4I163M18I841M1I607M2I2255M7D2372M10D269M7D749M1I2218M10D237M1D83M2I233M9I229M36D62M1D81M26I418M9I52M12D914M5D945M15D2479M6D406M14D58M36D150M1D152M8D2609M5D187M1I1886M11D1185M7D3942M3D315M1D226M10I1166M5I1121M2I341M15D1035M1I2081M2I339M108I1498M1D260M15D2013M6I671M3D32M1I148M1I16819M12D745M2I1788M4D3222M3I515M1D1191M15D537M4D82M2I219M2D1074M5D3992M1I61M6D27M1D549M81I290M1I2127M19D4617M10D3599M1I64M1D10266M2D639M5D438M10D434M1I3911M3I319M5D409M3D1223M4I483M1I1287M4I2160M12D4650M24D1320M1I5508M2D1223M9D1235M8I1207M7D1563M2I296M1I1268M1I2819M5D1798M1I34M1I663M13D526M7D1022M9D952M3D689M1I4543M14D588M2D24M4D4743M8D4029M10D616M5D2298M1I5986M1I2167M1D7042M 2 9600M29D33M20D1261M14I261M7D58M1D218M8D367M6D1769M6I518M2D82M7D246M7D159M25D842M3D365M6D1556M5D521M1I171M3D1802M18I987M11D8885M7D321M13D6M2D54M17I2031M8D628M5D29M1D681M1D3545M36I2624M6D144M8D36M3D881M1I1572M2D1274M15I221M6I2166M2D7M8D2033M81I653M3I2141M2D36M11D501M24D2208M15D8922M8I504M1I264M1I28M9I1982M1I966M10D2687M1D4176M17D1417M1D94M12D819M4D387M12I166M3D1231M1I717M12D757M6I84M6I1328M3I6715M1I1487M63I4234M2D470M1D607M1D64M1I336M2D790M1I332M9I541M3I1723M10D503M3D879M1D32M4D161M3D8M2I59M2I146M1I55M2D126M2D2150M16I1178M5D1309M3D195M1I53M1I136M9D234M46I137M12742I21M1I565M1D831M8D79M7I327M6D98M27D14M4D23M7I601M1I188M1D60M39I17M3D10M11D134M11D42M2D354M6D126M1I347M10D202M1I1504M6D3151M7I1188M1I410M10I4045M20I545M2D594M10D759M12D896M2D154M1I1202M15I139M5I652M36I1159M18I5417M8I5744M2I835M1D254M2D193M7D291M22D324M14I295M1I59M7D1794M8D37M7D7001M3I1870M1D13M1I14M1I1568M1D2170M1D3627M5D1030M3I176M2D203M18D211M4D136M1D25M2I249M12I4M1I400M3I423M1I528M7D182M1I889M1I481M1I3278M12D371M3I380M10D485M1I566M8D449M1I1910M2D434M3D33M3D577M2D661M7D193M1I19M4D135M37D187M3D446M2D373M9D355M3D229M4D3288M11024D1552M21D737M9D2306M6D816M3I44M9D1516M1I181M2D77M1D1836M2D43M2D989M14D540M1I23M2D2240M14D5582M3D506M15D89M5D16M1D950M17D418M3I422M7D593M2D7M2D630M1D2364M4I1283M1D269M1I1573M12D20795M14D510M13D464M1D34M27I2962M11D31M1I994M4D131M3D1521M6I1656M1I83M1I1678M22D22M14D506M21I1099M5D2192M1I1513M3D118M1I2170M1I53M6I2053M1I3258M1D12M1I47M5D2172M7558I3727M2I807M4D2332M42D995M1D3142M3D2005M1I583M4D4216M1I68M1D3891M2I2523M4D5444M20D20000M1D438M40I664M1D2333M2D1538M57I2969M41I13098M3I227M1I36221M2I4871M1I39M5D114M1D151M59D53M42D5366M1I1920M2I8071M1I2215M656D750M 3 2740M4D1011M18D67M24D15836M12D3357M2D1861M11D306M14D494M3D1097M42I314M2D1901M1D616M7D12M2D819M5D585M3D3179M2D306M48D10029M1D7584M15D1824M9D336M60I856M1D395M4D322M12D35M36I2343M90I11M1I511M1I318M2I676M2I386M1D103M2D8M5I684M12D137M21D10878M10D4176M1D1665M1D292M3I448M7I809M1I112M1D58M4D5489M16D62M1D225M1I24M36D11M1D124M1I17M51D38M2I294M10I2793M1D297M19388D7M3D52M1I1999M1D84M1D43M1D1170M1D4501M18D922M1D449M12I561M14I945M3D1096M1D137M2I592M1I1384M54D1340M20D1377M3D11M1I7M2I260M11D235M11I3432M2I936M7D275M9I7577M1I2080M3D543M1I609M1I4226M1I187M4I38M1D3231M1D583M1D25M1I244M1D614M3I2M3I20M2D269M1I178M1D388M10I6M4D3M1I343M1D16M1D1345M1I4116M2D4M1D127M1I223M1D2869M7I1055M9I131M11D155M7I545M2D2993M41I220M8D2472M2I2243M1D218M3I219M1I31M1I153M1I1323M8I6522M1D9M4D2015M1I991M1D1672M11D511M7D822M14D1823M6D709M4I563M7D519M9D2741M2I64M13I57M2D1565M18D220M1D727M11D454M13D614M2D3036M29I8525M3I2154M1I139M2D216M1I1052M1I1731M1I492M1D265M3D2389M2I150M1D269M1D58M5D1056M18D2559M28I286M1I2931M1D65M4D1761M1D18M1D335M1D1824M3D167M3D2358M1D5M1D3516M10D1692M1I2566M1D654M8D663M1I311M1D21M1D1916M3I157M42D653M10D25M4D94M21D333M1I4024M7D7152M3I221M18I17M1I51M14D2955M2D300M1D111M1I6367M9D2987M2I379M2D5925M12D1640M1I2611M1I176M1I2526M3I259M2D111M1D3225M1I226M4D614M6D742M1I240M3D710M1D11399M1D1357M1I23M1I1147M1D991M2I89M2D102M12D5533M8D2001M3D885M2I742M1D473M4D3012M23D178M5D2435M18D1510M2I187M9D1737M5440I59M1I53M1D1670M13D31M8D1232M4D1446M7D3493M3I209M8I1276M5475D1140M24D215M17I396M6D4395M1D75M9D8M2D2220M1D31M1I12M8D140M6D236M3I2018M8I421M1D145M10D671M1I2483M1I2392M9I91M19D92M40D2438M1I3679M2I5574M21D4498M1D2079M1D224M1D34M8D43M36I3M36D8963M12I1838M14D1083M1D3075M1D46M1I3164M10D400M1I1169M2D616M1I4091M11D988M8D202M8D8M3D376M1I153M6D379M1I90M4D95M 4 2068M1I339M6I2565M1D568M2I647M6D77M3I2788M4D319M15I7388M1D1929M1I318M8D916M4I198M2D20M1I602M3I740M326D1268M27I6186M11D135M39D710M20I630M9D745M1D5390M1I1472M1I259M6D2970M12I1453M2D826M1I729M1D550M8D1667M1D80M1I49M7I108M10D652M2D3441M6D548M4I2253M6I689M10D336M143D5333M27D1644M2I428M4D6770M2D7381M1D1946M2D679M10D522M44I906M21D3941M2D230M1I68M4I88M1D3521M1I527M3D959M1I6349M5I346M1I7144M1I297M1I52M3D225M1I554M14I45M1D299M11D1438M8D63M1D3142M9D5128M20D356M10I2718M22I4192M7D3625M13I2890M6D44M8I482M6I209M2I508M1I65M1D150M32I209M4I810M2I1333M6D3308M4D355M11D656M9I2500M6D39M34D1276M3I256M1D134M90D1307M7D1055M16D356M10I147M2I312M1I11302M1I408M3D36M20D3180M2I471M1I21M2D3770M2I389M2I32M2D1433M27I1332M3D295M4D6M11D671M1D1591M9I387M3D1020M18I344M17D2045M2I269M28D1564M8D339M32I97M20D17M3D60M13D55M1I180M4I92M1I107M2I2412M3D352M2I228M2D404M1D16M2D99M2I125M6I1231M10I146M1I309M1I3393M1D384M10D22M4D18M1D15M2I30M5D5396M2I21M1I65M4D205M1I1992M32I348M1D3409M2D583M1I251M1I1398M4D151M2I148M2I2153M3I452M1D1126M1D1846M3D789M6I416M4D39M5D1205M9I932M24I292M97D74M2I6M6D2426M4I324M32D93M8I453M9I2490M2I1999M1D467M3D183M3D2492M7D22M4I759M1D583M1D52M60D58M2I33M28I173M1I4479M2I598M15D207M1D51M3I902M27I637M6I858M2I7M8I193M1I396M7I3795M52D18M16I285M1I93M10I59M3D194M3I42M1D1900M11I1921M3326D94M2I416M7D1784M3D2312M1I152M9D1523M3D1576M1D155M11D2173M2I128M12I594M9D3023M 5 1093M12I1176M1D244M11D2707M1I1986M20I156M1I103M3D124M42D8444M1I120M15D11M1D454M8I1250M30D1244M2I1873M1D118M1I315M3I312M3D4043M1I1137M1D1899M23D750M18I175M30I2410M1I19M3I98M6I886M6D13M1I16M20I1041M3D1116M10D394M2D2168M14D3252M6D535M10D3000M6D1920M1D162M1D305M8I520M1I129M2D41M8I3424M8D51M4I135M3I678M1I507M6D325M14D28M3I300M2I261M6D1815M1I2189M1D287M1D4813M2D6M5D140M2D134M20I606M6D359M1I145M4D1066M11D130M12D2717M14D947M26D1035M3D591M18I1188M2I1344M20I723M24D661M18I1927M15D1994M7D391M30D343M36D480M20D4146M8D143M1D978M9D306M9D1078M6D1047M50I23M13I169M6I8M2I56M1D2324M15D3901M1I2706M12D388M1I1348M1I739M19D1050M6I8M2I4968M3I454M4D7929M6I1492M9D91M1D194M1I285M7D7136M1I155M1D479M1D2302M15I1980M1D237M1D1119M8D2033M11D87M1I918M1D29M3D5241M9D150M4I605M12I1145M14I2231M1I599M1I426M1I228M5I675M1D2628M1D11M9D390M5D8108M4D5426M11D823M10D20M8D1293M1I132M36D246M4D119M64I142M1I1213M9D410M13I676M1D1088M11D733M32I819M1D453M13D198M1D5645M84I4364M1I290M5500D652M2I5156M1I33M109I2801M2I522M16D50M1D151M28I1409M1D178M13D253M1I1474M19D18M2I360M2I800M1I2600M18D295M20D582M5D1883M6I818M14D28M12I445M2D946M1D173M35D2062M1I14944M1I2562M1I9M3D11M1D243M1D81M1I1015M1I1483M1I362M6I1184M27D230M2I200M14D378M1D66M3D164M6D173M4D20M9D941M1I153M1D1195M2I2241M9I358M1I49M10I478M1I255M4I42M3D3700M8201I1803M9D897M1D228M1D26M263D293M1D74M 6 977M1I6451M1I8266M1D5380M16D428M2I1494M1D410M18D126M1D3705M2I1038M18D995M11I149M1D265M5D827M3D1694M1I1619M11D42M37D576M8D4645M5D4889M3I54M7D80M1D3193M1I50M4I2696M4I4814M1D529M5D1451M1D89M5D501M2I3341M3D901M1I73M1D662M1I3788M1I388M2D11M1I2907M1D5M1I1289M4D16M5D325M1D4178M13I261M6I2M9D1334M12D443M1D1144M3D1580M20D306M9I1468M9I524M1I924M6D7916M2I866M3D348M1I5625M1D2204M4I2395M2133D12906M3I599M10D2989M2D1178M7D1015M1D41M4I156M3D1594M10D521M1I58M1D543M11D4516M3D140M3D3135M6I203M1D698M1D5481M18D256M3D45M11D85M4D131M9D3347M1D1040M5D122M2D42M6D692M4D4449M1D2253M10D5864M16D821M6D1009M1D71M1I1434M8D100M1I468M1I39M14D1815M6D327M1I115M1D1340M117I2676M8D2571M30I228M4D46M6D167M1I214M1D1664M5I140M52I55M2D13M6I225M46D101M1D293M13410I225M9D1915M18D81M24D6073M9D3031M1I4510M1D2692M1I304M1D1772M24D4704M3D1162M2I716M1I1357M1I306M1I13M1D1140M4D36M8D806M12D413M1I198M1D195M7508I91M14D6912M 6NA

绘制共线性的点图

###绘制共线性的点图 pdf("共线性的点图.pdf", width = 8, height = 8) dotplot(df) dev.off()

覆盖度

函数plot_coverage提供了一个有用的方法来查看一个给定的基因组有多少被包含在基因组排列中。默认情况下，它将目标染色体中的每个序列显示为一个矩形框，阴影区域代表目标基因组中被包含在比对中的部分。

p1 <- plot_coverage(df) p2 <- plot_coverage(df,fill = 'qname') p3 <- plot_coverage(df, fill = 'qname')+ scale_fill_brewer(palette = 'Set1') pdf("覆盖度.pdf", width = 8, height = 6) p1 + p2 + theme(legend.position = 'none')+ p3 + theme(legend.position = 'none') dev.off()

以上pafr代码部分来自：

Pafr官方文档：

https:// cran.r-project.org/web/ packages/pafr/vignettes/Introduction to pafr.html

点图：

ggplot(df, aes(alen, dv)) + geom_point(alpha=0.6, colour="steelblue", size=2) + scale_x_continuous("Alignment length (kb)", label = function(x) x/ 1e3) + scale_y_continuous("Per base divergence") + theme_pubr()

我们可以看到，有许多短的排列组合，其中一些是非常分歧的。但是也有一些非常长的排列，所有这些排列都显示出高度的相似性。因为pafr对象实际上是一个数据框架，我们可以再次使用标准的R函数来检查或分析它。例如，我们可以计算以查询基因组中每个序列为特征的排列的平均分歧水平。改变点阵图中序列的顺序默认的图是相当稀疏的，每个对齐的片段显示为一条暗线，查询和目标基因组中的序列边界显示为虚线。dotplot的附加参数让我们可以修改该图。例如，我们可以为每个查询和目标序列添加标签（label seqs），并改变目标序列的出现顺序。因为这里产生的点阵图是一个ggplot对象，我们也可以使用theme bw()来改变图的主题。

dotplot(df, label_seqs = TRUE, order_by = 'qstart') + theme_bw()

参数：参数order_by有三个可能的值。size"，"qstart "或 "supported"。 Size" （默认值）只是将查询序列和目标序列从大到小排列。 qstart "保持查询序列按大小排序，但按目标序列与查询序列的匹配程度重新排列。

基因密度图绘制

前言

应各位童鞋的留言，分享一下基因密度图的代码和数据，请需要的童鞋自取。代码和数据会上传到网盘。

-- du

公众号回复：20220411获得数据和代码

代码：

### 密度图 ## date: 2022.04.08 setwd("E:\\小杜的生信筆記\\基因密度图") #安装软件 #source("http://bioconductor.org/biocLite.R") #biocLite("gtrellis") #BiocManager::install("gtrellis") library(gtrellis) library(RColorBrewer) library("circlize") library("ComplexHeatmap") ## 导入数据 bed1<-read.table("Genome_length.txt",head=T,sep='\t') bed2<-read.table("gene_length.txt",head=T,sep='\t') ## gene密度大小转换 gene_density = genomicDensity(bed2,window.size = 1e5) col_fun = colorRamp2(seq(0, max(gene_density[[4]]),length=11),rev(brewer.pal(11, "RdYlBu"))) cm = ColorMapping(col_fun = col_fun) #绘制坐标 lgd = color_mapping_legend(cm, plot = TRUE, title = "",color_bar="continuous") ## 调整位置 gtrellis_layout(bed1,byrow = F,ncol = 1,xpadding = c(0.1, 0),gap = unit(2, "mm"), border = FALSE,asist_ticks=FALSE,track_axis = FALSE,legend=lgd) ## 绘制gene密度图 add_heatmap_track(gene_density, gene_density[[4]], fill = col_fun,track=1) ## 添加染色体名称 add_track(track = 1, clip = FALSE, panel_fun = function(gr) { chr = get_cell_meta_data("name") if(chr == "Chr9") { grid.lines(get_cell_meta_data("xlim"), unit(c(0, 0), "npc"), default.units = "native") } grid.text(chr,x =0.02, y = 0.38, just = c("left", "bottom")) add_track(track = 1, clip = FALSE, panel_fun = function(gr) { chr = get_cell_meta_data("name") grid.text(chr,x =0.02, y = 0.38, just = c("left", "bottom"))

ggcor包作图 | 相关性热图 | mental分析图

前言

在今天这个图形出来的时候就已经有大神在分享绘制的教程。在那时自己也获得绘制图形的代码（大概是在两年前吧）。最近，自己也看到几个比较好教程，在重新整理一下。

知乎主页 ： https://www. zhihu.com/people/xiaodu .com

B站主页 ： https:// space.bilibili.com/4462 19753?spm id from=333.1007.0.0

公众号：小杜的生信筆記

ggcor作图总结，本次教程主要针对ggcor包的作图来学习，主要图形如下，也许你在很多文章中看到过，这些图形在2020年左右出来时，在科研圈还是有很大的反响，以及很多大神就开发了ggcor这个作图包，得到很多人的使用。

图例

安装ggcor

ggcor包是比较难安装的，直接使用install.package()是无法安装的，在R的镜像中是无法存在的，很多大神都是将上传到GitHub中。如果你运气好，那么一次就可以安装，但如果你没有一点运气，那就需要花费很长的时间。我自己安转就花费了很多的时间，如果你按照我一下的方法也无法安装，那你就继续找”度娘“吧。

## # install.packages("devtools") devtools::install_github("houyunhuang/ggcor") # install.packages("devtools") devtools::install_git("https://gitee.com/houyunhuang/ggcor.git") ## ggcor下载 https://gitee.com/try-quiet/ggcor.git ## 下载下来后进行安装 devtools::install("~/R/x86_64-pc-linux-gnu-library/4.1/ggcor_master/",force = TRUE) devtools::install_local("~/R/x86_64-pc-linux-gnu-library/4.1/ggcor_master",force = TRUE)

教程一

教程网址： https:// mp.weixin.qq.com/s/7mYz gZngj-YetAQ0jCLfkw

--输出图形效果

代码区域：

## 导入相关R包 library(corrplot) library(vegan) library(ggcor)

本教程使用的数据是来自vegan包

> head(varechem) N P K Ca Mg S Al Fe Mn Zn Mo Baresoil Humdepth pH 18 19.8 42.1 139.9 519.4 90.0 32.3 39.0 40.9 58.1 4.5 0.3 43.9 2.2 2.7 15 13.4 39.1 167.3 356.7 70.7 35.2 88.1 39.0 52.4 5.4 0.3 23.6 2.2 2.8 24 20.2 67.7 207.1 973.3 209.1 58.1 138.0 35.4 32.1 16.8 0.8 21.2 2.0 3.0 27 20.6 60.8 233.7 834.0 127.2 40.7 15.4 4.4 132.0 10.7 0.2 18.7 2.9 2.8 23 23.8 54.5 180.6 777.0 125.8 39.5 24.2 3.0 50.1 6.6 0.3 46.0 3.0 2.7 19 22.8 40.9 171.4 691.8 151.4 40.8 104.8 17.6 43.6 9.1 0.4 40.5 3.8 2.7 > head(varespec) Callvulg Empenigr Rhodtome Vaccmyrt Vaccviti Pinusylv Descflex Betupube Vacculig Diphcomp Dicrsp Dicrfusc Dicrpoly 18 0.55 11.13 0.00 0.00 17.80 0.07 0.00 0 1.60 2.07 0.00 1.62 0.00 15 0.67 0.17 0.00 0.35 12.13 0.12 0.00 0 0.00 0.00 0.33 10.92 0.02 24 0.10 1.55 0.00 0.00 13.47 0.25 0.00 0 0.00 0.00 23.43 0.00 1.68 27 0.00 15.13 2.42 5.92 15.97 0.00 3.70 0 1.12 0.00 0.00 3.63 0.00 23 0.00 12.68 0.00 0.00 23.73 0.03 0.00 0 0.00 0.00 0.00 3.42 0.02 19 0.00 8.92 0.00 2.42 10.28 0.12 0.02 0 0.00 0.00 0.00 0.32 0.02 > dim(varechem) [1] 24 14 > dim(varespec) [1] 24 44

绘图：

mantel <- mantel_test(varespec, varechem, spec.select = list(Spec01 = 1:7, #依次定义四种物种作为Mantel的分析对象 Spec02 = 8:18, Spec03 = 19:37, Spec04 = 38:44)) %>% mutate(rd = cut(r, breaks = c(-Inf, 0.2, 0.4, Inf), labels = c("< 0.2", "0.2 - 0.4", ">= 0.4")),#定义Mantel的R值范围标签，便于出图 pd = cut(p.value, breaks = c(-Inf, 0.01, 0.05, Inf), labels = c("< 0.01", "0.01 - 0.05", ">= 0.05")))#定义Mantel检验的p值范围标签，便于出图 quickcor(varechem, type = "upper") +#绘制理化数据热图 geom_square() +#定义成方块状 anno_link(aes(colour = pd, size = rd), data = mantel) +#定义连线 scale_size_manual(values = c(0.5, 1, 2))+ guides(size = guide_legend(title = "Mantel's r",#定义图例 order = 2), colour = guide_legend(title = "Mantel's p", order = 3), fill = guide_colorbar(title = "Pearson's r", order = 4))

教程二

来源网址： https:// gitee.com/try-quiet/ggc or

本教程除了绘制Mental图外，还是用ggcor绘制相关性图形。

Correlation plot

## 导入相关的包 library(ggplot2) library(ggcor) ## 方块相关性 quickcor(mtcars) + geom_square()

## 圆形相关性 quickcor(mtcars, type = "upper") + geom_circle2()

参数设置：（参数设置与ggcorplt一致） type = " " 1、full:前部显示 2、upper:显示上部分 3、lower：显示下部分

一半显示图形，一半显示数字（P值，及显著性）

quickcor(mtcars, cor.test = TRUE)+ ## cor.test geom_square(data = get_data(type = "lower", show.diag = FALSE))+ geom_mark(data = get_data(type = "upper", show.diag = FALSE), size =2.5)+ geom_abline(slope = -1, intercept = 12)

#定义图例，在最后添加 guides(fill = guide_legend(title = "P value", order = 2))

参数设置：

##隐藏图例 guodes(fill=guide_legend(title=NULL)) ##修改图例 scale_fill_discrete(breaks=c("*","**","***")) ##颠倒图例顺序 guides(fill = guide_legend(reverse=TRUE))

Mantel Plot

该码与开始的码基本一致。

data("varechem", package = "vegan") data("varespec", package = "vegan") mantel <- mantel_test(varespec, varechem, spec.select = list(Spec01 = 1:7, Spec02 = 8:18, Spec03 = 19:37, Spec04 = 38:44)) %>% mutate(rd = cut(r, breaks = c(-Inf, 0.2, 0.4, Inf), labels = c("< 0.2", "0.2 - 0.4", ">= 0.4")), pd = cut(p.value, breaks = c(-Inf, 0.01, 0.05, Inf), labels = c("< 0.01", "0.01 - 0.05", ">= 0.05"))) quickcor(varechem, type = "upper") + geom_square() + anno_link(aes(colour = pd, size = rd), data = mantel) + scale_size_manual(values = c(0.5, 1, 2)) + scale_colour_manual(values = c("#D95F02", "#1B9E77", "#A2A2A288")) + guides(size = guide_legend(title = "Mantel's r", override.aes = list(colour = "grey35"), order = 2), colour = guide_legend(title = "Mantel's p", override.aes = list(size = 3), order = 1), fill = guide_colorbar(title = "Pearson's r", order = 3))

Circular heatmap

rand_correlate(100, 8) %>% ## require ambient packages quickcor(circular = TRUE, cluster = TRUE, open = 45) + geom_colour(colour = "white", size = 0.125) + anno_row_tree() + anno_col_tree() + set_p_xaxis() + set_p_yaxis()rand_correlate(100, 8) %>% ## require ambient packages quickcor(circular = TRUE, cluster = TRUE, open = 45) + geom_colour(colour = "white", size = 0.125) + anno_row_tree() + anno_col_tree() + set_p_xaxis() + set_p_yaxis()

绘制柱状图

导入数据

设置路劲

setwd("E:\\小杜的生信筆記\\20220505柱状图")

导入数据

df <- read.csv("data.csv", header = T) head(df) > head(df) Group Expression Sex 1 mpg 3.90 F 2 mpg 3.90 F 3 mpg 3.85 F 4 mpg 3.08 F 5 mpg 3.15 F 6 mpg 2.76 F

注意：本次数据是随机数据

绘制柱状图

自定义主题(非必须)

## 自定义图标主题（不是必须项） top.mar = 0.2 right.mar = 0.2 botton.mar = 0.2 left.mar = 0.2 ## 自定义主题 ## 合并上面的参数，将其合并成mythemel mytheme1 <- theme(panel.background = element_blank(), axis.ticks.length=unit(1.6,"mm"), plot.margin=unit(x=c(top.mar,right.mar,botton.mar,left.mar), units="inches")) #自定义主题2； #隐藏纵轴，并对字体样式、坐标轴的粗细、颜色、刻度长度进行限定； mytheme2<-theme_classic()+ theme(text=element_text(family = "sans",colour ="gray30",size = 12), axis.line = element_line(size = 0.6,colour = "gray30"), axis.ticks = element_line(size = 0.6,colour = "gray30"), axis.ticks.length = unit(1.5,units = "mm"), plot.margin=unit(x=c(top.mar,right.mar,botton.mar,left.mar), units="inches"))

绘制柱状图（教程一）

这个教程的好处就是，无需单独计算样本点均值等信息

#无需提前计算均值和标准差，直接绘制柱状图； ggplot(df, aes(x = Group, y = Expression, color = Group, fill = Group))+ stat_summary(fun = mean, geom = "bar", alpha = 0.5, width = 0.5, show.legend = F)

参数解析:

stat_summary()函数将自动计算每个样本的平均数 alpha = : 透明度 width:柱子宽度 show.legend:图例的显示与否

添加误差线

p2 <- p1 + stat_summary(fun.data = mean_sdl, fun.args = list(mult = 1/2), geom = "errorbar", width = 0.3, show.legend = F)

添加散点

p3 <- p2 + geom_jitter(shape = 16, alpha = 0.9, size =4, width = 0.15, show.legend = F)

在这里参数是非常简单的，不过多做解释。

添加定义主题

## 主题1 p3+scale_y_continuous(expand = expansion(add = c(0.1,0.3)))+ mytheme1

## 主题2 p3+scale_y_continuous(expand = expansion(add = c(0.1,0.3)))+ mytheme2

说实话，他给得这个自定义主题，并不喜欢，字体太小，各方面还是需要自己根据实际情况进行调整

#添加空心点； p3 <- p2+ geom_jitter(shape=1,alpha=1,size=2,width = 0.12,show.legend=F) #添加白心点； p3 <- p2+ geom_jitter(shape=21,fill="white", alpha=0.9,size=2,width = 0.12,show.legend=F) #应用自定义主题2； p3+scale_y_continuous(expand = expansion(add = c(0,0.3)))+ mytheme2

添加蜂群点

#尝试使用beeswarm包添加散点； geom_beeswarm(data = df, aes(color = Group), priority = "descending", alpha = 0.7, dodge.width = 0, size = 3, cex = 1.5, show.legend = F)

重要参数： priority = c("ascending","descending", "density", "random", "none") 默认为"ascending" 上升、升序； "descending" 降序、下降；而"density"朝向上下两个方向； "random"和"none" 效果相似，为随机状； dodge.width可用于设置图形左右两侧空余区域，避免图形被绘图区域剪切。

总结：

该柱状图的主要还是推荐绘制柱状图+散点图类型的。蜂群散点图不是那么的推荐。

绘制柱状图（教程二）

该教程需要计算均值和标准差，这个也是我前期一直使用的方法之一。

library(dplyr) df2 <- df %>% group_by(Group) %>% summarise( sd = sd(Expression), ## 计算标准差 Avg_Exp = mean(Expression) ## 计算均值 # A tibble: 4 x 3 Group sd Avg_Exp <chr> <dbl> <dbl> 1 gear 1.31 4.17 2 mpg 0.535 3.60 3 vs 0.971 2.67 4 wt 0.978 3.22

绘制图形

ggplot(df2, aes(x =Group, y = Avg_Exp))+ geom_bar(mapping = aes(color = Group, fill = Group), stat = "identity", alpha = 0.7, width = 0.5, show.legend = F)

添加误差线

ggplot(df2, aes(x =Group, y = Avg_Exp))+ geom_bar(mapping = aes(color = Group, fill = Group), stat = "identity", alpha = 0.7, width = 0.5, show.legend = F)+ geom_errorbar(mapping = aes(color = Group, ymin=Avg_Exp-sd, ymax=Avg_Exp+sd), width = 0.2, show.legend = F)

ggplot(df, aes(Group, Expression,fill = Sex))+ geom_bar(stat = "identity", position = "dodge") geom_text(aes(label = Expression))

这个图的数据是不对的，在这里只是顺便绘制一下，，具体图形数据需要你提前准备。

数量的统计

ggplot(df, aes(x = Group, fill = Sex))+ geom_bar(position = "dodge")

绘制差异基因富集柱状图

获得本教程原数据和代码

公众回复关键词：20220828 注：如需要获得以前的教程的代码和数据，请自行到知乎，公众号教程中查找关键词即可。（个人时间和能力有限，没有太多精力来帮你查看，谢谢！！）

问题咨询

如有问题咨询，可到知乎中付费提问题！！你问我答！

富集柱状图绘制

## 导入相关的包 library(clusterProfiler) library(enrichplot) library(ggplot2) library(plyr) library(dplyr

该教程需要进行对自己的差异基因进行注释，注释方法如下：

# ---------------基因注释------------------------- ## 方法一：使用R包msigdbr里的注释 ## 方法二：使用本地gmt文件 # 方法三：自己定义注释文件 # 如果你研究的领域比较新或小众，那么现有的注释文件可能不能满足你的分析需求。这时往往会自己总结一个基因集： # 查文献：你对某篇文章发现的某些基因感兴趣，例如参与某个生物学过程的几十上百个基因； # 序列比对：如果你的目标物种现有的注释太少，可以跟模式生物做序列比对，用相似序列的同源基因来注释目标物种的基因。读取当前文件夹内的所有gmt文件

读取当前文件夹内的所有gmt文件

files = list.files(pattern="*.gmt") gmts <- setNames(lapply(files, read.gmt), gsub(".v7.2.entrez.gmt", "", files))

数据处理函数

merge_result2 <- function(enrichResultList, output = "compareClusterResult") { if ( !is(enrichResultList, "list")) { stop("input should be a name list...") if ( is.null(names(enrichResultList))) { stop("input should be a name list...") x <- lapply(enrichResultList, as.data.frame) names(x) <- names(enrichResultList) y <- ldply(x, "rbind") if (output == "compareClusterResult") { y <- plyr::rename(y, c(.id="Cluster")) y$Cluster = factor(y$Cluster, levels=names(enrichResultList)) return(new("compareClusterResult", compareClusterResult = y)) y <- plyr::rename(y, c(.id="Category")) if (output == "enrichResult") { return(new("enrichResult", result = y)) if (output == "gseaResult") { return(new("gseaResult", result = y)) keep_category <- function(em_ORA, n) { table_em <- as.numeric(table(em_ORA$Category)) start <- rep(0, length(table_em) - 1) for(i in seq_len(length(table_em) - 1)) { start[i] <- sum(table_em[seq_len(i)]) showCategorys <- sapply(table_em, function(x) min(n, x)) start <- c(0, start) + 1 end <- start + showCategorys - 1 keep <- NULL for(i in seq_len(length(start))) { keep <- c(keep, c(start[i] : end[i])) return(keep) enrich_filter <- function(em_result, showCategory) { keep <- keep_category(em_result, showCategory) em_result <- em_result[keep, ] if ("NES" %in% colnames(em_result)) em_result$Count <- em_result$core_enrichment %>% strsplit(split = "/") %>% vapply(length, FUN.VALUE = 1) return(em_result) ## 作图函数 em_plot <- function(em_1 = NULL, em_2 = NULL, showCategory = 2, fill = "p.adjust", hjust = 1) { fill <- match.arg(fill, c("Category", "p.adjust", "log10_p.adjust")) result1 <- enrich_filter(em_1, showCategory) if (is.null(em_2)) { result <- result1 } else { result2 <- enrich_filter(em_2, showCategory) result2$Count <- -result2$Count result <- rbind(result1, result2) result$Category <- gsub("\n.*", "", result$Category) result$log10_p.adjust <- log10(result$p.adjust) data_plot <- result[, c("ID", "Category", "p.adjust", "log10_p.adjust", "Count")] data_plot2 <- data_plot data_plot2$ID <- factor(data_plot2$ID, levels = unique(data_plot2$ID)) data_plot2 <- plyr::rename(data_plot2, c("Count" = "gene_number")) h_just <- ifelse(data_plot2$gene_number < 0, -hjust, hjust) ggplot(data_plot2, aes_string(x = "gene_number", y = "ID", fill = fill)) + geom_col() + geom_text(aes_(x =~ gene_number + h_just, label =~ abs(gene_number)), color="black") + scale_x_continuous(label = abs, expand = expansion(mult = c(.01, .01))) + #两侧留空 theme_classic() + ylab("") + theme(axis.title.x = element_text(size = 15)) + facet_grid(Category ~ ., scales="free", space="free")

导入差异数据矩阵

##导出差异基因数据矩阵 df.fc <- read.csv("input.csv") colnames(df.fc)[1] <- "SYMBOL" ## 给第一列命名名称 df.fc[1:3,]

转换数据

#把gene symbol转换为ENTREZ ID df.id <- bitr(df.fc$SYMBOL, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = "org.Hs.eg.db") df.id[1:5,]

数据合并与排序

#让基因名、ENTREZID、foldchange对应起来 df.fc.id <- merge(df.fc, df.id, by="SYMBOL", all=F) df.fc.id[1:5,1:5] #按照foldchange排序 df.fc.id.sorted <- df.fc.id[order(df.fc.id$logFC, decreasing = T),] #获得ENTREZID、foldchange列表，做为GSEA的输入 id.fc <- df.fc.id.sorted$logFC names(id.fc) <- df.fc.id.sorted$ENTREZID

GSEA富集分析

# 自定义函数 gsea_func <- function(x, genelist, readable = FALSE) { GSEA_result <- GSEA(genelist, TERM2GENE = x, eps = 0) if (readable & nrow(GSEA_result) > 0) GSEA_result <- setReadable(GSEA_result, 'org.Hs.eg.db', #物种 'ENTREZID') #转回gene symbol return(GSEA_result) ## 富集分析 em <- setNames(lapply(gmts, gsea_func, id.fc, readable = TRUE), names(gmts)) %>% merge_result2(output = "gseaResult") em[1:5,1:3]

基础绘图

## 基础作图 em_plot(em, showCategory = 5, fill = "log10_p.adjust", hjust = 5)

基础绘图

em_plot(em, showCategory = 5, fill = "Category", hjust = 5) + scale_fill_brewer(palette="Set1")

有正负之分，可以把负的画到左侧，正的画到右侧

em_GSEA1 <- em_GSEA2 <- em res <- em@result em_GSEA1@result <- res[which(res$NES > 0), ] # 被激活的通路 em_GSEA2@result <- res[which(res$NES < 0), ] # 被抑制的通路 em_plot(em_GSEA1, em_GSEA2, showCategory = 2, fill = "log10_p.adjust", hjust = 10)

美化

em_plot(em_GSEA1, em_GSEA2, showCategory = 2, fill = "Category", hjust = 12) + scale_fill_brewer(palette="Set1")

绘制箱线图

直接使用ggplot()进行绘制

p5 <- ggplot(df,aes(x=Group,y=Expression))+ geom_boxplot(data=df,mapping=aes(color=Group,fill=Group), alpha=0.5, show.legend=F)

绘制小提琴图

ggplot(df, aes(x = Group, y = Expression, fill = Group))+ geom_violin(df, mapping = aes(color =Group, fill = Group), alpha = 0.4, trim = F, show.legend = F)

#自定义颜色； mycolor <- c("#0077c1","#00a99e","#6bc72b","#ff5a20","#ff1620","#752995") p1 + scale_colour_manual(values=alpha(mycolor,1))+ scale_fill_manual(values=alpha(mycolor,0.4))

叠加箱形图

ggplot(df, aes(x = Group, y = Expression, fill = Group))+ geom_violin(df, mapping = aes(color =Group, fill = Group), alpha = 0.4, trim = F, show.legend = F)+ ## 添加箱线图 geom_boxplot(df, mapping = aes(color = Group, fill = Group), width = 0.1, show.legend = F)

绘制山峦图

#添加分位线(中位线)； #线的样式可选solid、dashed、dotted、dotdash、longdash、twodash等。 #也可指定单一颜色，如vline_color="black"； ggplot(df, aes(x = Expression, y = Group, color = Group, fill = Group))+ geom_density_ridges(scale = 1.6, alpha = 0.3, rel_min_height=0.005, quantile_lines = TRUE, quantiles = 0.5, vline_size=0.3, vline_linetype="dashed")

注：出现这样的图形是由于数据的原因，不是所有的数据都适合绘制山峦图。

柱状图中添加差异标记

前言：

在一年以前，在做分析的时候，自己第一次使用R语言绘制柱状图，并模仿文献中的配色进行绘制。当时，绘制时候是把已经处理好的数据进行绘制。以及绘制图形的代码还是比较繁琐的。今天在公众号看到如何进行差异分析并添加显著水平标记？教程，是一个非常好的教程。自己在来学习一遍，希望自己可以理解得更透彻一点吧！

一、前期绘制的柱状图

这个图形自己还会比较喜欢的，不知道是否符合你的”口味“呢？

1.1 数据准备

这里的数据是已经进行了方差分析的数据，这是最开始时候绘制的，后面给他们做分析时，也做了从数据到图形的教程。

> head(df1) group value SD 1 CK 1.0000 0.073 2 1d 0.6879 0.073 3 3d 1.5000 0.100 4 6d 1.7399 0.233

绘图

## 加载相关的包 library(ggplot2) library(ggthemr) library(tidyverse) library(dplyr) library(ggpubr) library(devEMF) library(ggsignif) ## 导入数据 setwd("D:\\小杜的生信筆記\\20220505柱状图绘制") df1 <- read.csv("data_diffsing.csv", header = T) head(df1) names(df1) <- c("group","value","SD") compared <- list(c("CK","1d"),c("CK", "3d"),c("CK","6d")) ## 绘图 ggplot(df1, aes(x=group, y=value, fill = group))+ geom_bar(stat = "identity", position = "dodge",color = "black")+ scale_y_continuous(expand = c(0,0),limits = c(0,3.4))+ ##y从0开始 #更改柱子的颜色 scale_fill_manual(values = c("#0f0f0f","#5f5f5f", "#d7d7d7","#ffffff"))+ theme_bw()+theme(legend.position = "none")+ ## 添加误差线 geom_errorbar(aes(ymax = value+SD, ymin = value-SD), position = position_dodge(0.01),width = 0.2)+ ## 添加显著标记 geom_signif(annotations = c("*","**","***"), y_position = c(2.5,2.9,3.2),xmin = c(1,1,1), xmax = c(2,3,4))+ labs(y = "Relative expression leves", x = "", title = "gene 01")

教程二

这个教程来自"基迪奥公众号"，我一直认为基迪奥在做分析这块做的是比较好的。尤其是在做图这块，我是比较推荐大家可以关注他的。我不止一次说过，基迪奥做分析做得好。

我自己在公司里面做过分析，但是相比起来确实整体格局是不能先比较的。

2.1 导入数据

在这里数据所用的也是前面柱状图教程中的数据基于R 语言绘制柱状图教程混合整理.md 。在此，不再做阐述。

## 导入数据 dt <- read.csv("data.csv", header = T) head(dt) > head(dt) Group Expression Sex 1 mpg 3.90 F 2 mpg 3.90 F 3 mpg 3.85 F 4 mpg 3.08 F 5 mpg 3.15 F 6 mpg 2.76 F

2.2 方差齐性检验

## 方差齐性检验 nom <- bartlett.test(Expression~Group, data = dt) #如果是（多）因数，使用interaction()函数，leveneTest(y~interaction(var1,var2),data= data) #输出p值，当p>0.05时，方差是齐性的。 nom$p.value

2.3 单因素方差分析

oneway <- aov(Expression~Group, data = dt) anova(oneway) > anova(oneway) Analysis of Variance Table Response: Expression Df Sum Sq Mean Sq F value Pr(>F) Group 3 38.337 12.7789 13.138 1.664e-07 *** Residuals 124 120.611 0.9727 Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

2.4 组间多重比较

#组间多重比较； #TukeyHSD法（Tukey’s Honestly Significant Difference） com<-TukeyHSD(oneway,ordered =FALSE,conf.level = 0.95) ## 结果 Tukey multiple comparisons of means 95% family-wise confidence level Fit: aov(formula = Expression ~ Group, data = dt) $Group diff lwr upr p adj mpg-gear -0.5743750 -1.21647061 0.06772061 0.0968900 vs-gear -1.5005937 -2.14268936 -0.85849814 0.0000001 wt-gear -0.9536875 -1.59578311 -0.31159189 0.0010016 vs-mpg -0.9262187 -1.56831436 -0.28412314 0.0014881 wt-mpg -0.3793125 -1.02140811 0.26278311 0.4178404 wt-vs 0.5469062 -0.09518936 1.18900186 0.1239996

2.5 LSD法（Fisher’s Least Significant Difference）

#install.packages("agricolae") library("agricolae") #BiocManager::install("agricolae") out <- LSD.test(oneway,"Group",p.adj="bonferroni") #LSD法检验微小的差异，比较方便的是直接得出显著性字母标记，不需人工标记

> out $statistics MSerror Df Mean CV t.value MSD 0.9726721 124 3.413773 28.89007 2.681244 0.6610885 $parameters test p.ajusted name.t ntr alpha Fisher-LSD bonferroni Group 4 0.05 $means Expression std r LCL UCL Min Max Q25 Q50 Q75 gear 4.170937 1.3055696 32 3.825861 4.516014 2.300 7.300 3.30000 3.9500 5.050 mpg 3.596563 0.5346787 32 3.251486 3.941639 2.760 4.930 3.08000 3.6950 3.920 vs 2.670344 0.9710387 32 2.325267 3.015420 1.512 5.424 2.19750 2.4625 2.645 wt 3.217250 0.9784574 32 2.872173 3.562327 1.513 5.424 2.58125 3.3250 3.610 $comparison $groups Expression groups gear 4.170937 a mpg 3.596563 ab wt 3.217250 bc vs 2.670344 c attr(,"class") [1] "group"

2.6 绘制图形

library(ggplot2) library(ggpubr) ## 自定义图表主题，对图表做精细调整； top.mar=0.2 right.mar=0.2 bottom.mar=0.2 left.mar=0.2 #主题1； #隐藏坐标轴轴，并对字体样式、颜色、刻度长度等进行限定； mytheme1<-theme_classic()+ theme(text=element_text(family = "sans",colour ="gray30",size = 12), axis.line = element_blank(), axis.ticks = element_line(size = 0.6,colour = "gray30"), axis.ticks.length = unit(1.5,units = "mm"), plot.margin=unit(x=c(top.mar,right.mar,bottom.mar,left.mar), units="inches")) #主题2； #对字体样式、坐标轴的粗细、颜色、刻度长度等进行限定； mytheme2<-theme_classic()+ theme(text=element_text(family = "sans",colour ="gray30",size = 12), axis.line = element_line(size = 0.6,colour = "gray30"), axis.ticks = element_line(size = 0.6,colour = "gray30"), axis.ticks.length = unit(1.5,units = "mm"), plot.margin=unit(x=c(top.mar,right.mar,bottom.mar,left.mar), units="inches")) #自定义颜色； mycolor <- c("#0077c1","#00a99e","#6bc72b","#ff5a20","#ff1620","#752995")

ggboxplot(dt, x = "Group", y = "Expression", color = "Group", fill = "Group", alpha =0.2, add = "none")

散点图

## 散点图 ggboxplot(dt, x = "Group", y = "Expression", color = "Group", add = "jitter")+ scale_colour_manual(values = alpha(mycolor,1))

教程中的自定义颜色

ggboxplot(dt, x = "Group", y = "Expression", color = "Group", add = "jitter")+ scale_colour_manual(values = alpha(mycolor,1))+ mytheme1

2.7 添加显著性标识

#获取分组标签； labs <- unique(dt$Group) head(labs) > head(labs) [1] "mpg" "wt" "vs" "gear"

#使用ggpubr进行组间均值比较，默认wilcox.test法； #可通过method参数进行指定方法： compare_means(Expression~Group, dt)

# A tibble: 6 x 8 .y. group1 group2 p p.adj p.format p.signif method <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr> 1 Expression mpg wt 0.0316 0.063 0.0316 * Wilcoxon 2 Expression mpg vs 0.000000181 0.0000011 1.8e-07 **** Wilcoxon 3 Expression mpg gear 0.0784 0.078 0.0784 ns Wilcoxon 4 Expression wt vs 0.00388 0.016 0.0039 ** Wilcoxon 5 Expression wt gear 0.00653 0.02 0.0065 ** Wilcoxon 6 Expression vs gear 0.00000250 0.000013 2.5e-06 **** Wilcoxon

## 多组比较(one-way ANOVA test)：anova (parametric) 和 kruskal.test (non-parametric)； compare_means(Expression~Group, dt, method = "anova") #两组比较：t.test (parametric) and wilcox.test (non-parametric)； compare_means(Expression~Group, dt, method = "t.test") compare_means(Expression~Group, dt, method = "wilcox.test") ###------- # 指定参考组（对照组）； compare_means(Expression ~ Group, dt, ref.group ="mpg")

以上，差异分析就已经完成

compare_means(Expression~Group, dt, ref.group = "mpg", method = "t.test") #指定感兴趣得比较组； my_comparisons<- list( c("mpg", "wt"), c("mpg", "vs"),