拟南芥转录组完整流程2

发表于 | 分类于
| 字数统计: 5,805 | 阅读时长 ≈ 23

前言

这是根据“思考问题的熊“的课程《生信技能树转录组数据分析入门实战》进行整理的，这是一个标准的转录组分析流程，所用模式生物是拟南芥，整个课程信息如下所示：

课程题目：

生信技能树转录组数据分析入门实战

课程内容：

  1. 课程介绍
  2. Linux操作技巧
  3. 生物信息软件安装
  4. 转录组分析整体思路
  5. 转录组数据预处理
  6. 转录组数据序列比对
  7. 转录组数据表达定量
  8. 转录组数据差异分析

所用数据，如下所示：

mark

（二）inux操作技巧

Linux常用客户端

Xshell：用于登录远程服务器。
Putty：用于登录远程服务器。
WInSCP：用于下载服务器文件。

Linux命名

这个比较基础，略过，只讲一个alias
alias用于缩简一些命令，例如这次的练习数据我放在了/mnt/hgfs/biotest/rnaseq/这个目录下，每次都输入都太麻烦了，因此可以使用alias命名进行添加，如下所示：

1
2
cd #进行主目录
vim .bashrc # 进入bashrc文件

在.bashrc文件中添加如下字段：

1
alias wp='cd /mnt/hgfs/biotest/rnaseq/'

保存并退出，其中需要注意的是等号两边不要有空格，否则会出错。
此时直接输入wp就可以直接进行目录。

（三）生物信息软件安装

主要是利用conda进行安装，比较简单，详细的教程可以参照青山屋主的教程，这个是我看到的比较好的教程。

（四）转录组分析整体思路

作者回顾了一下以前的思路，如下所示：

大图如下所示：
mark

（五）转录组数据预处理

数据介绍

按照作者教程中的文件夹整理了所要分析的数据，如下所示：

现在解释一下这几个文件夹，其中：

  1. 00ref存放的是参考基因组文件（fasta格式）和注释信息（gtf/gff格式）文件；参考序列以及注释信息的下载地方通有Ensembl，JGI，此外，还有一些针对某些模式生物的专门数据库，例如araport是拟南芥的数据库，本文所用的数据就是拟南芥的测序数据（已经下载好，过程略）。

  2. 01_raw_data中的是原始数据（测序公司给你的数据），其中，sample1_R1.fastq.gz和sample1_R2.fastq.gz表示的是双端测序（PE），一个是左边的数据，一个是右边的数据。

  3. 02clean_data用于存放清洗了原始数据后的数据；

  4. 03align_data用于存放比对后的数据（例如BAM或SAM），以及分析的数据；

  5. others文件夹中的数据不清楚，这一部分的数据是从教程的主页上下载下来的，估计是分析后的数据，先不管。

解压缩参考序列文件

刚下载的参考序列文件与注释文件是以gz结尾的压缩文件，因此先进行解压，gunzip *gz的作用就是将压缩文件解压到当前文件夹，如下所示：


两个文件已经解压好了，一个是gtf文件，一个是fasta文件，可以用head命令查看一下序列：

此外，还能通过cat，less，more，head，tail这几个命令来查看序列，现在查看一下注释文件：

MD5值的计算

MD5主要用于检验数据传输后，前后两个数据是否完整，简单来说就是，测序文件很大，通过MD5我就能知道上传或下载后的数据跟我原来的数据是不是一样。检测过程如下所示：
生成4个原始数据的MD5值：

1
2
cd ../01raw_data 
md5sum *gz > md5.txt

上述两个命令是切换到原始数据的文件夹，通过md5sum命令是生成该文件夹下的所有gz格式文件的MD5值，并将这些MD5值给md5.txt文件，现在查看一下md5.txt文件less md5.txt：

检验MD5值

用md5sum命令来检验一下MD5的值：

如果文件完整，后面会显示OK。

测序文件的质控

质控软件是fastqc，如果没有安装可以利用ocnda进行安装，conda install fastqc，安装后，查看一下版本fastqc -v，如下所示：

第一种方法质控：常规方法

现在用fastqc进行质控，如下所示：
fastqc sample*gz，

运行过程与测序文件的大小以及计算机的性能有关，我们这个数据量不大，比较快，结束后，用tree来查看一下结果：

其中以html结尾的数据就是测序的结果，关于测序结果的解读，以前参考过青山屋主的《使用 FastQC 做质控》，可以看一下。

第二种方法质控：用for循环

此外，还可以用for循环来做，先把原来的zip和html文件删除掉rm -f *zip和rm -f *html，用于for i inls *gz;do fastqc $i;done来进行循环，同样能达到相同的目的。

第三种方法质控：利用管道命令

先把原来的zip和html文件删除掉rm -f *zip和rm -f *html，然后输入命名ls *gz |xargs -I [] echo 'nohup fastqc [] &'，这个命名是的运行过程就是列出raw_data文件夹下的所有测序文件（就是4个原始数据文件），然后导入到xargs命名中，利用xargs命令调用nohup命名。关于这个命名的解释：

  1. ls *gz是列出当前目录下的所有gz文件。
  2. |表示将上面列出的文件导入到下一步（xargs命名），这叫管道命令，很形象，就是将上面的结果像管道一下，输入到下一步。
  3. xargs命令是给其他命令传递参数的一个过滤器，也是组合多个命令的一个工具，使用-I指定一个替换字符串[]，这个字符串在xargs扩展时会被替换掉，当-I与xargs结合使用，每一个参数命令都会被执行一次，我们使用上述的这段命名时，它表示，用ls *gz列出的几个文件名字替换掉nuhup fastqc []的中括号。
  4. echo用于输出字符串；
  5. nohup是将操作置于后台，它要与&配合使用。
  6. 有关管道与xargs的用法，这个文章介绍得不错《xargs命令详解，xargs与管道的区别》。

再看一下流程：
ls *gz |xargs -I [] echo ‘fastqc []’ 就是将第一个[]里的内容替换掉第二个[]，如下所示：

将上述命令保存到fastqc.sh文件中，并运行fastqc.sh文件，如下所示：

结果如下：

第四种方法质控：multiqc

这种方法严格来说并不是做质控，而是把fastqc的质控结果给合并起来，直接在fastqc分析后的结果的目录中输入multiqc .就可以，如下所示：

打开multiqc_report文件，如下所示：

数据的过滤（去接头）

数据过程主要使用的是Trimmomatic软件，此软件的功能是用于去掉原始数据的接头序列，先看一下测序的质量文件中的Adaptor Concent，如下所示：

这张图衡量的是序列中两端adapter的情况，从图中我们可以看出这些测序文件的接头序列并没有去掉，因此我们要去掉这些软件，Trimmomatic就是实现这个功能的，它的一些介绍如下所示：

软件特点：

  1. 支持多线程，处理数据速度快；
  2. 主要用于去除Illumina平台接头；
  3. 根据碱基质量对fastq进行筛选；
  4. 支持SE和PE测序数据，支持gzip和bzip2压缩文件。

质量过滤的依据：

  1. ILLUMINACLIP：过滤reads中的Illumina接头；
  2. LEADING：从reads开头切除质量值低于阈值的碱基；
  3. TRAILING：从reads末尾切除质量值低于阈值的碱基；
  4. SLIDINGWINDOW：从reads的5’端开始，进行滑窗过滤，切掉碱基质量平均值低于阈值的滑窗；
  5. MINLEN：丢弃经过剪切后长度低于阈值的这条reads；
  6. TOPHRED33：将reads的碱基质量值体系转换为phred-33（一般都选择phred-33）；
  7. TOPHRED66:将reads的碱基质量值体系转换为phred-64。

软件使用
1. 接头序列的选择：

  • 软件默认提供了几个接头文件，在FastQC的结果中会提示是哪种接头，如下所示：

    文件显示的是TruSeq Adapter，因此要选择“TrueSeq3-SE.fa” and “TruSeq3-PE.fa”进行过滤（我们前面知道，本次案例中的是双端测序，因此用到的是TruSeq3-PE.fa）。

    • 但是如果显示的是Illumina Single End”/“Illumina Paired End”就要选择“TruSeq2-SE.fa”和“TruSeq2-PE.fa”，

2. 去接头参数的选择
选择true或false，随后会介绍。

trimmomatic使用案例

双端测序会输入2个文件，分别是正向测序序列和反向测序序列，过滤完数据后，会输出4个文件，双端测序都保留的是paired，如果其中一端序列过滤后被丢弃了，另一端序列被保留了，则为unpaired。

现在使用trimmomatic去掉adapter，如下所示：

1
2
3
4
5
6
7
8
9
trimmomatic PE -threads 4 \
sample1_R1.fastq.gz sample1_R2.fastq.gz \
../02clean_data/sample1_paired_clean_R1.fastq.gz \
../02clean_data/sample1_unpair_clean_R1.fastq.gz \
../02clean_data/sample1_paired_clean_R2.fastq.gz \
../02clean_data/sample1_unpair_clean_R2.fastq.gz \
ILLUMINACLIP:/home/biotest/ENTER/share/trimmomatic-0.36-5/adapters/TruSeq3-PE-2.fa:2:30:10:1:true \
LEADING:3 TRAILING:3 \
SLIDINGWINDOW:4:20 MINLEN:50 TOPHRED33

代码的解释：

  1. -threads:后面的数据表示是线程数，可不设置，程序会自行检测并选择核心数。
  2. ILLUMINACLIP:过滤 reads 中的 Illumina 测序接头和引物序列，并决定是否去除反向互补的 R1/R2 中的 R2。在此命令中输入的是/home/biotest/ENTER/share/trimmomatic-0.36-5/adapters/TruSeq3-PE-2.fa，它是接头序列文件，可以在trimmomatic的安装目录中找到，这一步要使用绝对路径；
  3. LEADING:3 TRAILING:3
  4. SLIDINGWINDOW:从 reads 的 5’ 端开始，进行滑窗质量过滤，切掉碱基质量平均值低于阈值的滑窗。
  5. LEADING: 从 reads 的开头切除质量值低于阈值的碱基。
  6. TRAILING: 从 reads 的末尾开始切除质量值低于阈值的碱基。
  7. MINLEN: 如果经过剪切后 reads 的长度低于阈值则丢弃这条 reads，被丢弃的 reads 数会被统计在 Trimmomatic 日志的 dropped reads 中。
  8. AVGQUAL: 如果 reads 的平均碱基质量值低于阈值则丢弃这条 reads。
  9. TOPHRED33: 将 reads 的碱基质量值体系转为 phred-33。
  10. TOPHRED64: 将 reads 的碱基质量值体系转为 phred-64。
    运行过程如下所示：

    这个过程根据序列的长度与计算机的性能有所不同，我在虚拟机中跑了大概5分钟，运行完后，如下所示：

这个结果显示，输入了250000个序列，其中196007个序列（78.40%）的序列留下，14.26%留在了一端，2.39%留在了另一端，另外，扔掉了4.94%的序列。如果将ILLUMINACLIP:/home/biotest/ENTER/share/trimmomatic-0.36-5/adapters/TruSeq3-PE-2.fa:2:30:10:1:true 中的true改为false，则结果如下所示：

这个结果显示，输入了250000个序列，其中90442个序列（36.02%）的序列留下，56.64%留在了一端，1.23%留在了另一端，另外，扔掉了6.11%的序列。
另外一个sample2样本也按照同样的操作（选true），去掉接头的所有文件如下所示：

（六）转录组数据序列比对

软件介绍

序列比对分为两种，一种是有参比对，另一种是无参比对。这个过程就是将公司测序得到的reads比对到参考基因组上，如果我们所研究的物种有参考基因组的话，就是将reads直接比对到参考基因组上就行，如果没有参考基因组，则先需要拼接出一个基因组，然后再进行比对。这两种方式的差别如下所示：


这张图里面涉及三种比对，分别是基因组，有参转录组，无参转录组，可以参考一下（图片出处：Conesa A, Madrigal P, Tarazona S, et al. A survey of best practices for RNA-seq data analysis. Genome Biology. 2016;17:13. doi:10.1186/s13059-016-0881-8.）

转录组测序的比对通常分为基因组比对和转录组比对两种，顾名思义，基因组比对就是把reads比对到完整的基因组序列上，而转录组比对则是把reads比对到所有已知的转录本序列上。无参转录组所用到的工具是trinity，有参转录本所用到的工具是STAR。

如果是基于基因组比对，所用的工具是STAR和Hisat2。如果是基于转录组进行比对，则可以使用RSEM（可以结合bowtie2和STAR进行比对和定量分析）。关于转录组所使用的一些工具2017年，Nature Com去年发一篇文章，讲的39个工具，可以看一下。

本案例所用工具是STAR。过程包括①建立索引；②进行比对（分为两种操作，分为简单版与复杂版）；③查看比对文件。

建立索引

这一步位于质控之后，比对之前。在比对的时候为什么要建立索引，这涉及到比对软件的一些算法。大致原理就是，如果直接将reads比对到参考基因组上，非常消耗计算资源，通常不这干。通常情况下就是，利用索引，将reads比对到参考基因组上，这样会比较节省计算资源，具体的原理我也不太清楚。

建立索引的步骤如下：

先建立一个目录，用于存放索引文物，即在rnaseq目录下，输入mkdir arab_STAR_genome，建立一个名叫arab_STAR_genome的目录。接着安装STAR软件，命令为conda install star，记得是小写。

现在切换到biotest/rnaseq目录下（这个目录就是00ref，arab_STAR_genome的上层目录，一定要在这个目录下，否则运行会出错），运行软件，代码如下：

1
2
3
STAR --runThreadN 6 --runMode genomeGenerate \ 
--genomeDir arab_STAR_genome \ 
--genomeFastaFiles 00ref/TAIR10_Chr.all.fasta \ --sjdbGTFfile 00ref/Araport11_GFF3_genes_transposons.201606.gtf \ --sjdbOverhang 149

运行过程中出现了错误，STAR的进程被杀死了，因为我用的是虚拟机，虚拟机的内存是1G，而STAR这个软件很耗内存，因此把虚拟机的内存上调到8G，再次运行，出现finished successfully即表示运行结束，如下所示：

进程还是被杀死，因此最后切换到笔记本的Ubuntu中进去，STAR能顺利运行，运行结束后，如下所示：

可以查看一下arab_STAR_genome文件夹，下图就是运算的结果 ：

比对

接着进行第二步，先进行简单版的比对，简单版的比对输入的参数比较少，以比对sample1为例 说明，如下所示：

1
2
3
4
5
STAR  --runThreadN 5 --genomeDir arab_STAR_genome \
--readFilesCommand zcat \
--readFilesIn 02clean_data/sample1_paired_clean_R1.fastq.gz \
02clean_data/sample1_paired_clean_R2.fastq.gz \
--outFileNamePrefix 03align_out/sample1_

在03align_out文件 中可以看到这些文件 ：

SJ文件包括了一些可变剪接的信息，还有3个Log文件，它们包括了一些处理过程的信，sample1_Log.final.out含有了一些统计信息，如下所示 ：

sample1_align.out.sam文件则是比对的结果，例如是否有错配，哪些reads比对上了。接着用复杂版进行比对 ，
这次用sample2进行复杂比对，如下所示：

1
2
3
4
5
6
7
8
STAR  --runThreadN 5 --genomeDir arab_STAR_genome \
--readFilesCommand zcat \
--readFilesIn 02clean_data/sample2_paired_clean_R1.fastq.gz \
02clean_data/sample2_paired_clean_R2.fastq.gz \
--outFileNamePrefix 03align_out/sample2_ \
--outSAMtype BAM SortedByCoordinate \
--outBAMsortingThreadN 5 \
--quantMode TranscriptomeSAM GeneCounts

出现successfully表示比对结束，如下所示：

比对后的结果如下所示：

与sample1的简单比对相比，sample2的复杂比对结果中，少了一个sam文件，而多了2个bam文件，其中一个是经过坐标转换的，形成的转录本文件（sample2_Aligned.toTranscriptome.out.bam），另外一个是排过序的bam文件（sample2_Aligned.sortedByCoord.out.bam）。此外，还有一个文件，叫sample2_ReadsPerGene.out.tab，这个文件所含的信息是，一个基因上含有多少个reads。可以查看一下这个文件，如下所示：

第一列是基因名称，后面3列是基因的reads数。

现在查看一下sample2排序后的数据，所用命令为samtools view 文件名，如下所示：

1
samtools view sample2_Aligned.sortedByCoord.out.bam | head

转录组数据表达定量：

处理原始比对文件
picard/samtools
将sam格式转换为bam格式
对bam文件进行排序
去除比对得分比较低的序列
如果需要（当质量不好时）可以去除重复reads

定量表达的流程：

STAR+RSEM：输出结果选择转录本定量或者基因定量；定量单位包括feature count,FPKM，TPM；操作相对复杂；先比对再定量

STAR+HTSeq：输出结果为原始read count；结果可用于差异表达分析;操作相对 简单；先比对再定量

Kallisto（free-alignment）：速度快省内存；基于转录本定量；不产生Bam文件，不方便基它后续分析（不进行比对，直接定量）。

新建一个文件夹，arab_RSEM

1
rsem-prepare-reference --gtf 00ref/Araport11_GFF3_genes_transposons.201606.gtf \ 00ref/TAIR10_Chr.all.fasta \ arab_RSEM/arab_rsem

开始 运行，如下所示：

生成了7个文件，如下所示：

其中比较重要的是arab_rsem.transcripts.fa，这个文件提供了转录本的信息;接着，再新建一个目录 ，命名为04rsem_out，然后输入以下命令：

1
rsem-calculate-expression --paired-end --no-bam-output --alignments -p 5 -q 03align_out/sample2_Aligned.toTranscriptome.out.bam arab_RSEM/arab_rsem 04rsem_out/sample2_rsem

运行结束，一共生成 了3个文件
，如下所示：

在这三个文件中，sample2_rsem.stat含有统计信息，这些信息可以用于进一步的分析，例如可视化。而sample2_rsem.genes.results则是基因的结果，可以查看一下，如下所示：


其中，第1列是基因id，第2列是一个基因对应的转录本.

接着看一下基于转录本的信息，即文件sample2_rsem.isoforms.results的信息，如下所示：

第1列就变成了转录本的id，第2列就成了基因的id。
以上是STAR+RSEM的定量流程。

用kallisto进行定量。

如果没有安装kallisto，则用conda进行安装，安装后，新建一个目录，arab_kallisto，进入该目录，随后建立 索引 ，如下所示：

1
2
3
mkdir arab_kallisto
cd arab_kallisto
kallisto index -i arab_kallisto ../arab_RSEM/arab_rsem.transcripts.fa

开始运行，运行过程如下所示：

整个过程可能需要几分钟，最终会生成一个叫arab_kallisto的文件，这个文件就是后面进行定量的索引文件，如下所示：

现在再新建一个目录，用于存放定量的结果，如下所示：

1
2
mkdir 05kallisto_out
kallisto quant -i arab_kallisto/arab_kallisto -o 05kallisto_out/sample2 02clean_data/sample2_paired_clean_R1.fastq.gz 02clean_data/sample2_paired_clean_R2.fastq.gz

这个过程很快，结果如下所示：

然后进入05kallisto_out目录，里面还有一个称为sample2的目录，进入，如下所示：

sample2目录里有3个文件，其中tsv文件是可读的，可以查看一下，如下所示：

数据差异分析：

STAR+featureCounts是STAR+HTSeq的升级版，在这一部分中，利用STAR的结果，通过featureCounts来计算。

表达结果转换为表达矩阵。

为了方便进行差异分析，我们需要根据前文提到的内容，用STAR的复杂分析，将所有的样本都转化为相应的BAM格式，原始数据有sample1和sample2，2个样本又有2个生物学重复，因此是4个数据，分别为sample1，sample1_2，sample2和sample2_2，转换后的数据如下所示：

安装subread工具，conda install subread
随后进行分析，计算结果非常快，如下所示：

1
featureCounts -p -a ../00ref/Araport11_GFF3_genes_transposons.201606.gtf -o out_counts.txt -T 6 -t exon -g gene_id sample*Aligned.sortedByCoord.out.bam

03align_out中会出现2个文件，分别为out_count.txt与out_counts.txt.summary。可以查看一下out_counts.txt文件，如下所示：

这个文件就是定量的结果。

将定量的结果转换为表达矩阵。

根据前文的案例，将所有的样本都转化为表达矩阵，如下所示：

1
rsem-generate-data-matrix *_rsem.genes.results > output.matrix

第1列是基因名称，第2到第5列是样本的整合，如下所示：

有的基因表达量是0,随后可以将基因表达量为0的基因去掉，并且将表头进行更改，

1
awk 'BEGIN{print "geneid\ta1\ta2\tb2\tb2\n"}{if($2+$3+$4+$5>0)print $0}' output.matrix > deseq2_input.txt

运行后，查看一下原来的行数
wc -l output.matrix与wc -l deseq2_input.txt，如下所示：

上图中，去掉基因表达量为0的基因前，矩阵是37337行，去掉后为17971行，后续用deseq2_input.txt进行分析 。现在 建立一个06deseq_out的文件夹，进行差异分析，并将deseq2_input.txt复制过去：

1
2
mkdir 06deseq_out
mv deseq2_input.txt ../06deseq_out/

接着用R语言进行分析，代码如下所示：

加载各种所需要的包：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
getwd()
# 获取当前目录
setwd("D:\\netdisk\\biotest\\rnaseq")
# 将工作目录定位到测序结果的目录
source("http://bioconductor.org/biocLite.R")
# 调用biocLite
biocLite("DESeq2")
biocLite("edgeR")
# 安装分析数据所用的包
install.packages("ggrepel")
library(DESeq2)
library(edgeR)
library(ggrepel)
library(ggplot2)

导入数据：

1
2
3
4
5
6
7
8
9
10
11
input_data <- read.table(
  "D:\\netdisk\\biotest\\rnaseq\\06deseq_out\\deseq2_input.txt",
  header = TRUE, 
  row.names=1)
# 导入数据，其中添加row.names=1是将第1列设为矩阵的名称
input_data <- round(input_data,digits = 0)
# reads数都是整数，因此要保留整数部分，否则后面会出错
input_data <- as.matrix(input_data)
condition <- factor(c(rep("ctl",2),rep("exp",2)))
coldata <- data.frame(row.names = colnames(input_data),condition)

构建差异分析矩阵：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
dds <- DESeqDataSetFromMatrix(countData = input_data,            
                              colData = coldata,
                              design = ~condition)
# dds是分析矩阵
dds <- DESeq(dds) # 进行差异分析
res <- results(dds,alpha = 0.05) # 提取分析结果
summary(res)
# 查看分析结果
names(resdata)[1] <- "Gene"
res <- res[order(res$padj),]
resdata <- merge(as.data.frame(res),
                 as.data.frame(counts(dds,normalized = TRUE)),
                 by = "row.names",sort = FALSE)
write.table(resdata,file="diffexpr_results.txt",
            sep="\t",
            quote = F,
            row.names = F)
plotMA(res)

M是纵轴，表示log2 fold change，A是和横轴，代表标准后的cont的平均值，红点表示padj值小于0.1，落在绘图区域外的值用小三角表示。

p-value设定为NA的情况：

  1. 在同一行中，所有样品的count值都为0，这时baseMean为0，log2 fold change、p-value以及adjusted p-value都可设为NA
  2. 如果一行中包含了一个落在可信区域之外的极值，p-value 以及 adjusted p-value可设定为NA 。极值的确定可根据Cook距离，也可自行设定。
  3. 如果一行因为含有过低的mean normalized count而被自动独立过滤（automatic independent filtering）过滤掉了，这时可把adjusted p-value设定为NA

构建maplot函数：

1
2
3
4
5
6
maplot <- function(res,thresh = 0.05, labelsig = TRUE, ...){
  with(res, plot(baseMean, log2FoldChange, pch=20, cex=.5, log="x",...))
  with(subset(res,padj<thresh),points(baseMean, log2FoldChange, col="red",pch=20, cex = 1.5))
}
maplot(resdata,main="MA plot")

构建的maplot函数的绘图：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
resdata$significant <- as.factor(resdata$padj<0.05 & abs(resdata$log2FoldChange) >1)
ggplot(data=resdata, aes(x=log2FoldChange, y = -log10(padj),color = significant)) +
  geom_point()+
  ylim(0, 8)+
  scale_color_manual(values = c("black","red"))+
  geom_hline(yintercept = -log10(0.05),lty=4, lwd=0.6,alpha=0.8)+
  geom_vline(xintercept = c(1,-1),lty=4,lwd=0.6,alpha=0.8)+
  theme_bw()+
  theme(panel.border = element_blank(),
        panel.grid.major = element_blank(),
        axis.line = element_line(colour = "black"))+
  labs(title="Volcanoplot_biotest(by RVDSD)", x="log2(fold change)",
       y="-log10(padj")+
  theme(plot.title = element_text(hjust = 0.5))+
  geom_text_repel(data=subset(resdata, -log10(padj) > 6),
                  aes(label = Gene),col="black",
                  alpha = 0.8)

ggplot2绘制的火山图：