1 为什么要进行生信分析

在如今的生物的相关实验过程中，会产生海量的实验数据，通过使用计算机相关的软件操作帮助我们进行数据的解析，有助于我们寻找更关键的数据。

1.1 生信的标准分析流程是什么？

根据使用的需求不同，生信的分析步骤也会有些微差异，但较为通用的一般可以归纳为几个常用的基础步骤（尽管使用的软件可能互不相同，但仍未跳出该步骤的原理）：

对基因以及转录组而言：
- 实验设计
  明确实验的需要——我们到底需要进行什么分析？（如差异表达分析、基因组变异分析）。
  选择适当的样本和根据实验需要选择适当的技术（如RNA-seq，WGS等）。
  根据实际情况选择样本的数量与重复的次数。
- 数据生成
  经过了实验的进行，一般会产生大量的原始数据生成，包括但不限于：测序数据、蛋白质谱、以及其他数据类型。
  有些可能需要自己提取原料后交予公司检测，有些则根据要求不同，使用方法不同，具体的情况优具体决定。
- 数据预处理
  质量控制（qc）：一般而言如果是交由公司进行数据生成的情况下，公司可能会完成最基本的质控。
  数据过滤：过滤一些“噪音”区域以提高下游分析的准确率，如低表达基因与低覆盖区域。
- 数据比对与组装
  将于处理后的数据与参考基因组比对，或进行基因组的组装
  - 比对（Alignment）：将处理后得到的测序数据（RNA-seq）比对到参考基因组上（基因组比对）或转录组（转录组比对）上。工具：hisat2。
  - 组装（Assembly）：在没有参考基因组的情况下，将得到的测序数据组装为一个基因组。工具：trinity（转录组）、SPAdes（基因组）。
- 功能注释与富集分析
  我们进行了数据比对或基因组、转录组的组装后，还需要通过生信的数据库进行蛋白质和基因的功能注释，注释后可以对组转与比对的结果进行样本的功能和通路的富集分析，通过该分析可以明确其生物学意义。
  - 基因功能注释
    为基因注释生物学功能，常用工具：BLAST、Ensembl。
  - GO富集分析
    通过GO数据库分析差异差异基因的功能差异，常用工具：clusterprofiler。
  - KEGG通路分析
    分析基因在该生物通路中的作用，常用工具：KEGG。
- 可视化
  将数据以图表方式直观的展示，以便于解释结果与交流。
  - 火山图：展示差异表达基因的显著性与表达变化倍数（log2foldchange）。
  - 热图：展示基因样式的整体表达情况。
  - 主成分分析：用于检测样本间的总体差异。
    在R中的常用工具：ggplot2、pheatmap。

以上步骤上相对于基因以及转录组而言比较常用的实验过程，后续则是对结果进行一个总结汇报，撰写为文章，后续还有关于该结果的验证实验。生信分析贯穿了实验的整个过程，在每一个步骤都有对应的工具以及使用方法。希望后续使用到该流程的同门可以学习以上的工具的基础使用方法。

2 转录组分析

将你实验与处理后得到的RNA-seq数据进行比对到参考基因组与转录组，进行功能的注释后进行富集分析，最后可视化。以下为转录组分析的一些通用的流程与工具介绍。

2.1 什么是转录组分析？

转录组原理是研究和测量细胞、组织或生物体在特定条件下的全部RNA分子（即转录本）的表达水平和组成。转录组代表了基因组在特定时刻或条件下活跃转录的部分。转录组分析旨在揭示基因表达的动态变化，并探索不同生物状态下的基因调控机制。通过RNA逆转录为cDNA后，对cDNA进行测序，得到基因表达活跃度，直观意义上，同样长度的基因，被阅读到的次数越多，则意味着他们的表达更活跃，表达量更多。

2.2 关于RNA-seq数据

公司返还的RNA-seq测序数据一般包含以下几种基础文件是我们需要的：测序结果文件（*.fq.gz），与样本信息表。

关于样本信息表：一般最低需要包含以下一些列：样本名、测序数据名、是否为双端测序文件，双端测序文件名。
以下为一个样本的例子

这是sample.txt中的信息
ZJFGonad1	T20210624_00036	T20210624_00036_1.fq.gz
ZJFGonad1	T20210624_00036	T20210624_00036_2.fq.gz
ZJFGonad2	T20210624_00037	T20210624_00037_1.fq.gz
ZJFGonad2	T20210624_00037	T20210624_00037_2.fq.gz
ZJFGonad3	T20210624_00038	T20210624_00038_1.fq.gz
ZJFGonad3	T20210624_00038	T20210624_00038_2.fq.gz
ZJDDMGonad3	T20210624_00041	T20210624_00041_1.fq.gz
ZJDDMGonad3	T20210624_00041	T20210624_00041_2.fq.gz
ZJMGonad2	T20210624_00043	T20210624_00043_1.fq.gz
ZJMGonad2	T20210624_00043	T20210624_00043_2.fq.gz
ZJMGonad3	T20210624_00044	T20210624_00044_1.fq.gz
ZJMGonad3	T20210624_00044	T20210624_00044_2.fq.gz
SYFGonad1	T20210624_00045	T20210624_00045_1.fq.gz
SYFGonad1	T20210624_00045	T20210624_00045_2.fq.gz
SYFGonad2	T20210624_00046	T20210624_00046_1.fq.gz
SYFGonad2	T20210624_00046	T20210624_00046_2.fq.gz
SYFGonad3	T20210624_00047	T20210624_00047_1.fq.gz
SYFGonad3	T20210624_00047	T20210624_00047_2.fq.gz
SYMGonad1	T20210624_00048	T20210624_00048_1.fq.gz
SYMGonad1	T20210624_00048	T20210624_00048_2.fq.gz
SYMGonad2	T20210624_00049	T20210624_00049_1.fq.gz
SYMGonad2	T20210624_00049	T20210624_00049_2.fq.gz
SYMGonad3	T20210624_00050	T20210624_00050_1.fq.gz
SYMGonad3	T20210624_00050	T20210624_00050_2.fq.gz
ZJFBrain1	T20210624_00051	T20210624_00051_1.fq.gz
ZJFBrain1	T20210624_00051	T20210624_00051_2.fq.gz
ZJFBrain2	T20210624_00052	T20210624_00052_1.fq.gz
ZJFBrain2	T20210624_00052	T20210624_00052_2.fq.gz
ZJFBrain3	T20210624_00053	T20210624_00053_1.fq.gz
ZJFBrain3	T20210624_00053	T20210624_00053_2.fq.gz
ZJDDMBrain1	T20210624_00054	T20210624_00054_1.fq.gz
ZJDDMBrain1	T20210624_00054	T20210624_00054_2.fq.gz
ZJDDMBrain2	T20210624_00055	T20210624_00055_1.fq.gz
ZJDDMBrain2	T20210624_00055	T20210624_00055_2.fq.gz
ZJDDMBrain3	T20210624_00056	T20210624_00056_1.fq.gz
ZJDDMBrain3	T20210624_00056	T20210624_00056_2.fq.gz
ZJMBrain1	T20210624_00057	T20210624_00057_1.fq.gz
ZJMBrain1	T20210624_00057	T20210624_00057_2.fq.gz
ZJMBrain2	T20210624_00058	T20210624_00058_1.fq.gz
ZJMBrain2	T20210624_00058	T20210624_00058_2.fq.gz
ZJMBrain3	T20210624_00059	T20210624_00059_1.fq.gz
ZJMBrain3	T20210624_00059	T20210624_00059_2.fq.gz
SYFBrain1	T20210624_00060	T20210624_00060_1.fq.gz
SYFBrain1	T20210624_00060	T20210624_00060_2.fq.gz
SYFBrain2	T20210624_00061	T20210624_00061_1.fq.gz
SYFBrain2	T20210624_00061	T20210624_00061_2.fq.gz
SYFBrain3	T20210624_00062	T20210624_00062_1.fq.gz
SYFBrain3	T20210624_00062	T20210624_00062_2.fq.gz
SYMBrain1	T20210624_00063	T20210624_00063_1.fq.gz
SYMBrain1	T20210624_00063	T20210624_00063_2.fq.gz
SYMBrain2	T20210624_00064	T20210624_00064_1.fq.gz
SYMBrain2	T20210624_00064	T20210624_00064_2.fq.gz
SYMBrain3	T20210624_00065	T20210624_00065_1.fq.gz
SYMBrain3	T20210624_00065	T20210624_00065_2.fq.gz

表格所展示的信息就是公司测序返还的组名与样本名的对应关系。

3 linux

生物信息学的数据所占用的储存空间都较大，且需要运算的数据量十分庞大，需要多线程核心与高额内存配合，且数据处理的大部分软件在linux环境下才可以运行，因此进行生信处理需要使用服务器进行数据处理。

linux是一门以C语言为底编写的编程语言降低了语言的使用门槛，代码更为简洁，功能内容页更为丰富，以便于人们使用。

一些最基础的使用命令在系统的使用linux前，你只需要学会一些最基础的命令，让自己可以使用linux，而不需要使用linux进行编写程序，以便快速开展实验。以下为一些代码的简单介绍，具体需要什么代码还可以根据需要查询linux命令手册

  ｜—命令—           ｜-内容 -             ｜
  ｜cd <path>        ｜进入一个目录        ｜
  ｜ls <dir or file> ｜列出路径内的内容    ｜
  ｜sh .sh           ｜运行sh脚本          ｜

如果你需要查看一些命令的使用说明可以使用 command -h -help以查询该命令的详细说明

3.1 如何进入服务器

首先我们要如何进入服务器？我们不可能随时都在服务器机房中使用服务器主机进行操作，所以通常情况需要使用到使用者电脑的ssh命令，以远程连接服务器，在自己的电脑主机上使用linux命令操纵服务器。
什么是ssh？
ssh是一种远程连接服务的协议名字，ssh命令的含义即为，使用ssh协议进行远程连接，使用者不需要理解ssh协议到底是什么内容，我们只需要负责使用即可。

在MacOS上，由于Mac和linux同为ARM架构下的产物，因此在终端中自带ssh命令，可以直接通过ssh 用户名@服务器的ip地址进行访问
在Windows上，则有所不同，由于其为x86架构的产物，在windows的命令行，其运行的环境为封闭环境，不提供ssh命令，因此需要借助一些软件来使得Windows可以使用ssh命令远程连接服务器。
- 下载X shell7 ，该软件为付费使用软件，但具有免费版本，可以给学校和学生免费获取授权，点击进入官网下载，下载后，注册邮箱并使用邮箱激活后即可使用。
完成以上步骤后，我们就已经可以远程连接来操作服务器了，那我们要如何传输文件呢？此处使用filezilla，点击进入官网下载，并选择Download Filezilla Client版本进行下载。

在做好以上的软件准备后，就可以进入服务器了。在命令行中，用户名的开头表示的就是当前使用的用户账号，以及所处目录[用户名@服务器名所处当前目录]$

3.1.1 关于命令行模式的一些基本介绍

什么是命令行？什么是图形化界面？
在计算机的执行层面上，实际上背后运作的逻辑并没有大家日常所用的计算机界面一样那么简洁直观，大家所使用的民用计算机系统的界面是各种图形化以后展示的结果。
可以通俗的理解为大家用鼠标点击打开一个文件夹这种操作，背后的逻辑则是程序工程师们将计算机可以识别的指令直观化，图形化为了用鼠标左键双击两下即可以点开一个文件夹。也就是简单的理解为命令行的cd 一个目录 等同于大家使用windows或者mac的时候的鼠标双击打开一个文件夹。

区别在哪？
图形化界面就是大家用可视化的操作简单点击，系统通过预设好的一些命令行，将命令行发送到需要处理的方位，也就是说计算机操作的本质是对这命令行使用这个系统编写时设定好的命令。而面向大家的使用界面则是工程师们经过不懈努力编写成的可以用对图像的操作简单快速的调用这个操作背后编写好的命令发送给计算机，让计算机进行响应。
- 图形化听起来那么好，那么有什么缺陷吗？
  当然有！他的使用者没有那么大的权限，只能进行工程师们设定好的一些操作，而直接操作命令行则可以完成除了图形化操作以外可以完成的一些功能。举个例子吧，或许他并不是那么合适：你对着文件夹A点击了两下，那么你只能进入文件夹A对吗？但是在linux中，进入目录的只是一个基本指令cd，你可以通过cd进入除了A以外的别的文件夹。虽然他丧失了直观可视的图形化操作，但是给予了用户更大的使用方便性。
我们为什么要学习命令行？
因为服务器版的系统并不是简单给予民用的，面对的对象不是毫无计算机基础的用户，因此服务器版本的系统并不带有图形化界面（或许有些带有，但是我们实验室的不具有图形化界面），需要使用者通过基本命令来调用服务器上的资源。
关于文件夹与目录
实际上大家在系统中看到的文件夹图标样式的所有“文件夹”，本质上都是在命令行中定义为dir（即为目录）。
- 想象一下，假设你有一个文件夹A，又在文件夹A里面创建了一个文件夹B，在文件夹B中有一个txt文件c.txt文件。
  此时，这就已经有一个目录了，即在命令行中，理解为你在目录A中有一个目录B，在目录B中有一个文件c.txt。这很好理解吧？即“文件夹”等同“目录”
- 那么什么叫路径（PATH）呢？
  路径即为，这个目录，或者这个文件，在这个计算机中存放的地方（这个描述并不准确，但是在linux或者windows以及mac中可以简单这样认为）。再用一下上面的例子吧，假设你登陆这个服务器的用户名叫AA，服务器名字叫fuwuqi，那么你可以简单认为，这个文件c.txt的存放的位置是：/AA/A/B/c.txt
  这个是十分具有逻辑性对吧，打个比方吧，你叫AA,你在这个名字叫fuwuqi的房间中，有了一块属于你的空间，那么和你有关的东西应该存在在属于你的空间内，这很好理解吧？然后你在你的空间内，放了一个文件袋A，这个文件袋里又放了一个文件袋B，然后里面放了一页写了字的纸，你给他取名叫c。在计算机命令行模式中，就表现为/AA/A/B/c.txt，这个例子希望可以有助于你理解什么叫路径。1
- 如果你理解了路径就是在计算机中，用以给计算机识别某个文件所在的位置的一种表达方式，那么你选择还需要知道，在计算机中，关于路径有很多种表达方式。
  首先计算机系统中，有一个专门存放用户数据的文件夹，一般叫user,给你分配的空间（也就是是叫你的用户名的文件夹）就是建立在这个user目录下的一个次级目录，表现形式通常为/user/用户名/。
  但是计算机系统的运行还有系统自己独立的一些文件夹，他们通常是独立为一个其他的文件夹，所以其实你在你用户名下，计算机真正识别的路径其实是/服务器目录/user/用户名/ 这种从服务器目录出发的路径又称为绝对路径。
  既然有绝对，那就也会有相对，因为一般登陆的服务器，或者使用自己的计算机的时候都需要先创建一个用户对吧？你在使用这个用户名的时候就可以以～/表示/服务器目录/user/用户名/，也因为他是根据你登陆的用户名而有的路径表达方式，所以也叫相对路径，他可以帮助你在复杂的多级目录下，快速帮助你定位到你当前用户名的位置，以便后续路径的定位。但是你要注意文件需要保存在你使用的用户名的后续目录中，因为～/不能指向别的用户，以及用户的母文件夹以上的路径。
- 除了～/这一种相对路径以外，还有别的更简单的路径表达方式吗？有的，比如当你现在正处在目录A这一等级中，你可以使用./快速表示自己所在的位置，如果你要进入目录B则可以写为cd ./B，这极大的节约了你输入路径的时间，以及降低了你因为路径错误带来的麻烦。
  又或者，我如果想快速回到上一级目录呢？则使用../可以快速回到当前目录的上一级目录，当你在目录A中输入了cd ../ 则可以快速回到/用户名/这一级目录。
  那假如我在多级目录下，想快速回到用户的主目录（也就是/用户名/）有办法吗？我们可以直接输入cd后续不带任何路径，他会默认将我们带回用户的主目录下。

以后我们还能学到很多别的linux命令行，这教程主要并不是阐述linux的使用，关于linux的基础学习可以参考以下视频。

3.2 linux在服务器中需要用到转录组分析软件

由于前期大量的数据处理需要占用大量内存与运算资源，因此，使用linux服务器进行数据分析是目前前期处理的主要方式。

3.2.1 前期准备,软件下载

让我们开始学习如何使用linux服务器吧，以一个转录组分析为开始。
首先请让我们登陆服务器使用ssh命令登陆我们的服务器吧 ssh example@192.168.1.1
登陆成功后显示为这样
(base) [yonghuming@fuwuqi suozaimulu]$括号内即为我们所在的服务器环境，目前是基础环境，用户名则显示为用户名@服务器名，～表示当前所处的目录，这个时候的～表示是当前用户的家目录（也就是主目录）。
图示为登陆后界面的介绍。
我们需要先准备进行后续处理的软件

软件的下载
- 在linux中提供手动下载的方式，使用 wget http://xxx.命令可以将网上提供的软件包下载到你所在的目录，此处我们下载一个’conda’软件试试。 wget https://repo.anaconda.com/archive/Anaconda3-2024.06-1-Linux-x86_64.sh #这个是conda完整版的下载地址，大小约1G，相对占用较多空间，但更稳定
  wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh #此为miniconda下载地址，具有更小的体积以及基础功能，不如完整版稳定
  在下载完成后可以获得一个.sh文件，含义为可执行shell脚本，可以通过sh命令运行该脚本，直接安装conda
  sh Anaconda3-2024.06-1-Linux-x86_64.sh #sh命令使用方法为 sh *.sh或如果你使用minicondash https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
  之后一路按回车，直到他要求输入yes，即可完成安装。
- 使用环境管理工具（如conda）进行安装使用命令conda install xxx
- 我们为什么需要使用管理工具进行安装？如果是使用linux提供的下载方式，我们只能获得安装包（如.sh，.gz等文件），还需要我们手动进行安装，而在linux中，我们安装成功安装包后，还需要设定一个环境变量来使得你可以在不同的目录下输入这个软件的名字都可以让系统识别。
  以conda为例我们需要使用
```
export PATH=$PATH:/path/to/conda/bin  
#使用export命令，添加路径变量(path),使得linux可以识别conda这个字符是一个软件，自动调用这个路径里的软件。
#这个命令只对你当前这一次登陆服务器有效，当你下线后失效，为了使其永久生效，需要将export命令写入~/.bashrc中，~的含义为当前用户的主目录（又称家目录）
conda更新后，在conda被安装时会自动写入路径变量，但不代表其他软件手动安装时会自动写入，还是需要手动写入的。
```
```
   # >>> conda initialize >>>
   # !! Contents within this block are managed by 'conda init' !!
__conda_setup="$('/ME4012_Vol0002/Ahome_Dir/gonggong/bio_conda/anaconda3/bin/conda' 'shell.bash' 'hook' 2> /dev/null)"
if [ $? -eq 0 ]; then
eval "$__conda_setup"
else
if [ -f "/ME4012_Vol0002/Ahome_Dir/gonggong/bio_conda/anaconda3/etc/profile.d/conda.sh" ]; then
    . "/ME4012_Vol0002/Ahome_Dir/gonggong/bio_conda/anaconda3/etc/profile.d/conda.sh"
else
    export PATH="/ME4012_Vol0002/Ahome_Dir/gonggong/bio_conda/anaconda3/bin:$PATH"
fi
fi
```
为.bashrc文件中，conda自动写入的环境变量，请观察export行。
怎么样查看文件呢？
使用less命令less .bashrc可以在一个窗口查看文件，按q可以退出查看。
使用cat命令cat .bashrc可以将文件打印在当前窗口，当文件过大时会一直滚动窗口，此时按ctrl+c可以取消进程，ctrl+c（control+c）都是通用的取消指令。
如果只想看这个文件头部几行呢？使用head命令 head -n 1 .bashrc就会打印出.bashrc的第一行。末端则是使用tail命令，原理是一样的。
- 我们要如何在linux中编辑文件？
  使用vim命令，我们可以打开vim编辑器进行文件编辑，vim .bashrc即可打开图片编辑随后以英文输入法按i，可以进入编辑模式,按esc退出，返回到visual（查看）模式。
我们需要使用的软件有哪些？只对于这一次转录组分析而言，我们需要使用hisat2与samtools软件在使用conda安装前，需要为conda添加搜索软件的仓库（channels）输入命令conda config --add channels conda-forge依靠同理我们还需要添加r、bioconda两个库，如果anaconda网站很卡，那么可以使用清华镜像站：https://mirrors-i.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/ 的bioconda，同理也可以在清华镜像站中找到其他库。使用conda命令安装conda install hisat2、conda install samtools

3.3 进行参考基因组与转录组比对

我们拿到了转录组测序数据后，第一时间需要进行比对，在此我推荐使用hisat2软件进行基因组的比对。进行比对我们首先需要获得该物种的参考基因组、基因组注释文件以及样本测序数据。

如何获得参考基因组首先参考基因组需要什么文件呢？需要，FASTA文件（参考基因组文件）以及gtf或gff文件（基因组注释文件）可以在以下网站寻找自己实验物种的基因：
NCBI数据库:https://www.ncbi.nlm.nih.gov/
Ensembl数据库：https://useast.ensembl.org/index.html
UCSC数据库：https://genome.ucsc.edu/

通过以上的数据库寻找自己实验的物种的参考基因组。又或者可以选择自己组装的基因组。

3.3.1 hisat2的使用

hisat2是一个比对软件，可以将目的基因mapping到参考基因组上，得出比对的差异数据。

在使用hisat2前，需要先使用之前搜集到的参考基因组构建一个hisat2使用的比对文库。
使用hisat2-build程序构建。

hisat2-build [options]* <reference_in> <ht2_index_base>
以下是代码的解释
[options]   hisat2-build运行的设置参数
<reference_in>    参考基因组输入
<ht2_index_base>  hisat2比对文库路径

建库成功以后，可以使用hista2进行比对。
hisat2的基本使用方法：

hisat2 [options]* -x <ht2-idx> {-1 <m1> -2 <m2> | -U <r>} [-S <sam>]
以下是代码的解释
[option]        hisat2运行的一些额外参数，使得hisat2有其他拓展性的功能
-x <ht2-idx>    ht2-idx表示hisat2构建的参考比对库 -x表示参考库的路径
-1 / -2         对于双端测序的文件，有1.fq.gz 2.fq.gz两个测序文件，使用对应路径
-U              单端测序文件路径    
[-S <sam>]      输出为.sam格式文件

仅仅依靠基本的参数，hisat2也已经可以运行进行比对，但是仍然需要一些参数达到更强大的功能（比如让他调动更多的线程进行运算，具有更快的运算速度，又或者规定他的输出路径），以下提供一些别的参数可以在需要的情况下使用。

[option]
-q                    输入的文件为fq/fastq（默认）
-f                    输入的文件为fa/fasta
-I/-- minins <int>    最小读长，仅在无碎片比对中有效
-X/-- maxins <int>    最大读长，同上。
--summary-file <path> 输出比对总结到指定文件
--new-summary         将比对总结以一个对计算机识别更友好的风格输出到指定文件中
-p                    使用线程数（无法超过cpu线程数）
-help                 帮助文档
-version              查看hisat2的版本

完成以上步骤后，你将得到一些.sam文件，为比对结果总和，由于sam文件过于庞大，但sam格式的文件是人类可读文件，我们需要将其转换为bam的二进制文件，精简其信息。以便于后面的输出。

在此之前，我们如何得到这个软件的运行情况呢？可以通过hisat2 [options]* -x <ht2-idx> {-1 <m1> -2 <m2> | -U <r>} [-S <sam>] 1>file.log 2>&1将软件的运行日志输出为一个log文件，以下是对log文件的一些解释。

1 表示标准输出，即这个软件的标准运行结果会反馈到该日志中。
2 表示错误信号，当软件运行出现错误时，会讲错误日志输出。
2>&1 表示将错误日志输出到标准输出的日志中，并添加到日志末尾，而不是覆盖该日志。以便更容易理解软件的运行。

3.3.2 sam文件的转换

sam文件转换，我们一般使用samtools软件

samtools <command> [option]

将sam文件成功转换为bam文件以后，还需要建立bam文件的索引bai文件，他会自动被读取，但没有该索引的情况下无法正确导向到bam文件。
sra-tools软件同样具备该功能 sra-tools需要一些额外的参数完成这些工作，接下来将解释参数的含义

<command>
sort 将比对文件排序
index 构建比对索引

[option]
-o 输出文件路径

总的来说samtools的使用方式较为明了，将其变为bam文件后有利于后续的处理。

3.3.3 将gff3文件转换为gtf文件

有时候公司返回的测序数据中的注释文件为gff3文件而非gff文件，在生信处理中,由于后续软件需要使用到的注释文件格式为信息内容更为简单的.gtf文件,因此需要将GFF3文件转换为.gtf格式。

使用gffread软件。
gffread [-g <genomic_seqs_fasta> | <dir>] [-s <seq_info.fsize>]

可以通过-T参数将gffread工具变为,将gff变为gtf格式
-T ——main out put will be GTF instead of gff3 使用-o 参数指定输出路径 gffread input.gff3 -T -o output.gtf

gff3文件示例

## zjChr1   AUGUSTUS    gene    53688   61862   .   -   .   ID=zjg1;
## zjChr1   AUGUSTUS    mRNA    53688   61862   0.97    -   .   ID=zjg1.t1;Parent=zjg1;
## zjChr1   AUGUSTUS    stop_codon  53688   53690   .   -   0   ID=zjg1.t1.stop1;Parent=zjg1.t1;
## zjChr1   AUGUSTUS    CDS 53688   53714   1   -   0   ID=zjg1.t1.CDS1;Parent=zjg1.t1;
## zjChr1   AUGUSTUS    exon    53688   53714   .   -   .   ID=zjg1.t1.exon1;Parent=zjg1.t1;
## zjChr1   AUGUSTUS    intron  53715   53879   1   -   .   ID=zjg1.t1.intron1;Parent=zjg1.t1;
## zjChr1   AUGUSTUS    CDS 53880   53938   1   -   2   ID=zjg1.t1.CDS2;Parent=zjg1.t1;
## zjChr1   AUGUSTUS    exon    53880   53938   .   -   .   ID=zjg1.t1.exon2;Parent=zjg1.t1;
## zjChr1   AUGUSTUS    intron  53939   54026   1   -   .   ID=zjg1.t1.intron2;Parent=zjg1.t1;
## zjChr1   AUGUSTUS    CDS 54027   54143   1   -   2   ID=zjg1.t1.CDS3;Parent=zjg1.t1;

gtf文件示例

## zjChr1   AUGUSTUS    transcript  53688   61862   0.97    -   .   transcript_id "zjg1.t1"; gene_id "zjg1"
## zjChr1   AUGUSTUS    exon    53688   53714   .   -   .   transcript_id "zjg1.t1"; gene_id "zjg1";
## zjChr1   AUGUSTUS    exon    53880   53938   .   -   .   transcript_id "zjg1.t1"; gene_id "zjg1";
## zjChr1   AUGUSTUS    exon    54027   54143   .   -   .   transcript_id "zjg1.t1"; gene_id "zjg1";
## zjChr1   AUGUSTUS    exon    54815   54893   .   -   .   transcript_id "zjg1.t1"; gene_id "zjg1";
## zjChr1   AUGUSTUS    exon    54970   55038   .   -   .   transcript_id "zjg1.t1"; gene_id "zjg1";
## zjChr1   AUGUSTUS    exon    55282   55346   .   -   .   transcript_id "zjg1.t1"; gene_id "zjg1";
## zjChr1   AUGUSTUS    exon    55641   55698   .   -   .   transcript_id "zjg1.t1"; gene_id "zjg1";
## zjChr1   AUGUSTUS    exon    56065   56145   .   -   .   transcript_id "zjg1.t1"; gene_id "zjg1";
## zjChr1   AUGUSTUS    exon    56231   56267   .   -   .   transcript_id "zjg1.t1"; gene_id "zjg1";

通过该文件示例我们可以看到，gtf的文件统一将命名改为了ggene_id，统一了基因的id。

3.3.4 进行功能注释

将转录组mapping到基因组后，需要进行表达的定量分析，此处需要使用到除了linux以外的另外一种编程语言R。

3.4 什么是R

R是一种专注于数据处理与图形化的现代轻量级语言，R中具有丰富的关于生信处理的库（如bioconducter等），还有一个优质的管理插件BiocManager，因此对于处理后的生信数据分析一般使用R语言进行分析。

数据类型在各种不同的编程语言中，都对我们输入的字母、数字与符号有规定的类型，但是古早的开发型语言，对于数据的类型有严格的规定，而建立在这些开发型的较为底层的语言之上的应用型语言则对于数据的类型具有较为宽松的，不同于开发型语言的规定。
如之前linux的一些软件中所标注的类型举例： int 又称为整数型数据，即为只能输入整数数值，无法输入字母，字符，以及带有小数点的浮点数。
string 为字符串类型，什么是字符串呢？即为一串你输入的任何字符，字符串生成后实际就无法改动了，他们变为一个常量。
希望对于一些数据类型的讲解可以帮助你理解，数字，字符在计算机中到底是怎么定义的，但深入讲解则超出了本文想要表达的含义，大体你可以理解为，纯数字的一串，和带有字母的一串。
实际上对于计算机的理解而言，哪怕他们位数相等，但是他们这两个数据类型是不相等的，字符串类型即可以包含数字，也可以包含字母等，而int只能包含数字。因此实际上他们最大的上限位数不同，也就是这两种数据类型能储存的上限不同，在更为底层的如C语言中还有如<float>等其他数据类型。但我们不需要理解这个原理。
说了那么多，R语言中的数据类型是什么呢？
分别为：
numberic–数值型
integer –整数型
character –字符型
logical –逻辑型
complex –复杂型
raw –原始型
知道了解即可，实际上并不需要理解。这需要与后续的变量设定关联起来进行说明。
变量
由于在linux中，使用者只需要根据已经设定好的程序操作，输入即可，不需要涉及有关变量的操作，因此挪到R语言这一章中进行说明。
- 什么是变量？
  变量是一种用于存储数据的命名实体，数据可以是数值、字符串、对象等不同类型。你可以为变量赋值、读取它的值，或者在程序执行过程中改变它的值。
  变量是一个抽象的概念，实际应用中，变量被称为是一个具有名字的，有特定的规则属性的一个储存单元。
  而这个名字由变量的创造者命名。变量的创造并不是指创造出储存单元，而是从储存空间中划分出一个储存单元用于这个变量的数据，此时变量就被视为已经创造
- 变量具有什么规则属性？
  首先需要明确，变量由于需要储存数据，因此变量中的数据有数据类型规定数据的类型，而对于变量（即这个储存单元本身）有对象类型来规定变量的类型。
- 什么是对象类型。
  对象类型是用户自定义的一种复合数据类型，它封装了数据结构和用于操纵这些数据结构的过程和函数。因此本质上其也是数据类型，但是在R中，对象类型被简单化了。可以理解为这个变量的属性。
- 对象类型有哪些？
  vector –向量(R中默认的变量类型)
  list –列表
  matrix –矩阵
  array –数组
  factor –因子
  data.frame –数据框
  scalar –标量
- 它有什么含义吗？
  明确了对象类型的几种类型以后，你应该可以明白它的含义了——他规定了变量中的数据的表现形式。
  如是以列表的方式，列成表输出。
  还是以数组的方式，成为一个数字组输出。
  取决于你对这个变量的用法的需求。
  以下对变量的一些对象类型进行一些举例.
  在图中可以清晰直观的看到，matrix类型中，只可以有同一类型的数据类型，而在data.frame格式中却可以有多种数据类型，这便是对象类型的规则。

3.4.1 为什么要使用R?

R是一个开源的编程语言，它被创造出来的用途一开始就是用于进行统计分析以及数据的可视化处理。其代码由整个R语言社区一起维护管理，在此基础上，衍生出了一个R语言专门用于生物信息方面处理的社区Bioconduct，并为R语言编写了很多用于生信处理的软件，极大的降低了用户使用时的学习成本，规避了真正的编程问题（即不是创造程序，而是使用代码）。

4 以一个例子开始R的使用之旅

4.1 下载R

R语言在服务器中也可以进行部署，我们的后续步骤会需要使用到linux系统下的R语言环境，因此需要先下载R语言（如服务器中没有）
此处推荐使用conda直接安装R。 conda install r-base
安装成功后在服务器中输入命令R即可进入R语言环境
当我们不需要时使用q()即可退出R。
我们已经获得了一个由hisat2软件进行比对后产生sam文件，又经过samtools重新排序后生成的含有基因比对信息的.bam文件
sam文件包含的信息读者了解即可，现在已经有成熟的程序可以自行进行后续处理，我们只需要调整程序的参数。
获得一个以sam文件处理得到的bam文件后，我们需要进行对比数据的定量分析。

4.2 什么是定量分析，为什么要进行定量分析？

1.什么是定量分析：定量分析是一种通过收集、处理和解释数字数据，以分析和理解现象或问题的研究方法。
2.为什么要进行定量分析：在生信中我们只是得到了基因比对的差异，但是我们仍需要一个标准来解释到底什么基因发生了什么样的变化，是基因的表达上调了？还是下调了？都需要通过进行定量表达的分析后得出结论,下一步对基因表达量片段进行统计。

4.3 定量分析需要做什么(程序准备)？

首先我们需要准备R语言环境中需要的代码库（或者说软件包）

R #在linux中输入该字，进入R语言环境
install.package("BiocManager")    #insatll.package是R语言中安装程序包的命令，安装的BiocManager是生信软件管理工具，后续使用该管理工具进行生信相关程序包下载。
BiocManager::install("Rsubread")
BiocManager::install("DESeq2")
BiocManager::install("edgeR")    #使用BiocManager工具下载"Rsubread"、"limma"、"edgeR"三个软件包

这三个软件包具体有什么用？
请点击“Rsubread”,“edgeR”，DESeq2
经过该步骤后，我们已经准备好了进行定量统计的R语言软件包，可以开始进行定量统计了。
此处需要Rsubread包中的核心函数featureCounts，其主要用途为进行测序数据中计数比对到基因组或其他特征的读段，统计比对数据。

## 加载必要环境
library("limma")
library("Rsubread")
library("edgeR")
#使用featureCounts进行基因计数，并输出为一个数据框文件（fc_list），参数的解释（输入ban文件，注释使用文件 =输入gtf文件，使用GTF进行注释，是双端测序）
fc_list <- featureCounts(*.bam,annot.ext = *.gtf,isGTFAnnotationFile = TRUE, isPairedEnd = TRUE )

经过以上步骤我们可以生成一个数据框表格，但是只有该计数结果不足以作为最终的结果，因为缺少标准化，由于基因长度与测序深度不同，需要对定量结果进行统计学处理以得到标准结果,我们需要一个共同认可的标准化参数TPM/FPKM来解释基因表达的特异性。

什么是TPM/RPKM？
- 是两种常用于 RNA-Seq 数据分析中的基因表达量标准化方法。它们都旨在解决不同样本间测序深度和基因长度差异对表达量计算的影响。RPKM与TPM都是先将不同样本长度（b）中得到的样本表达量去除以样本长度（a），得到样本基因概率学上表达等级（c）。（此步骤去除样本长度带来的表达量差异）
- RPKM得到c后，将c除以所有样本总表达量（d），得到RPKM差异值
- TPM得到c后，将c除以c中表达量之和（e），得到TPM值
虽然RPKM在后续处理中使用的是样本总表达量，在统计学计算中样本总表达量应视为：$样本单位平均表达量（c）*基因长度（b）$，即我们假定样本取样是均匀的，越长的基因理应得到更多的reads比对结果。因此，对于RPKM的处理中，在最后一步重新加入了基因长度的影响。
但是对于我们需要的结果为：基因表达长度影响中，我们可以得到当整样本总表达数量之和去用c去除时，对于结果为：样本差异化表达参数该计算时，样本结果不受基因长度的影响，得到的数据展示的对比结果在结论上不受影响。

下一步为将fc_list进行标准化处理

fc_ref <- cbind(fc_list$annotation[,1],fc_list$counts,dg_rpkm,dg_tpm)
fc_ref <-setNames(data.frame(fc_list$stat,fc_list$counts,dg_rpkm,dg_tpm),c("gene_id","counts","rpkm","tpm"))
write.table(fc_ref,output,sep = "\t",row.names = FALSE,col.names = TRUE)  ##将featureCounts结论与定量结果tpm/rpkm输出为文件，output为自己设定的输出路径

4.4 合并表达矩阵

此时生成的文件为数个进行了定量表达后的文件，我们需要将他们合并为一个表达矩阵以便于后续的差异性分析，此时先退出R语言环境，回到linux并进行表达矩阵的合并。

ls *.count > genes.quant_files.txt #将count矩阵文件名全部列出，并输出合并为一个样本定量信息表。
perl ~/script/abundance_estimates_to_matrix.pl --est_method featureCounts --quant_files genes.quant_files.txt --out_prefix genes 
#第二行命令为使用perl脚本abundance_estimates_to_matrix.pl 统计方式使用featureCounts 输入文件为合并输出的那个文件名，输出的文件名为genes

该步骤主要使用了perl脚本，perl语言是一个文本处理十分强大的古老语言，其结构样式比较复杂，将会在后续阐明其用法。

4.5 差异分析

我们获得了基因的表达矩阵，需要对其的基因表达量的差异进行分析，找出其中差异的表达部分，和表达矩阵一样，都需要使用到文本处理强大的perl语言。

perl /ME4012_Vol0002/bio_soft/public/miniconda3/bin/run_DE_analysis.pl --matrix genes.counts.matrix --method METHOD --samples_file sample2.txt

在该步骤中主要调用perl脚本来调动函数进行我们的差异分析，其主要模块为method模块，该模块可以选择使用两个软件进行分析，分别是DESeq2与edgeR，这两个都是差异分析使用的软件，前后参数都为进行输入的文件，matrix为输入的合并矩阵，samples2.txt需要一些特别说明。

关于样本信息表
假设我们具有10个样本，分别是样本1-10，在寄出时样本1-5我们标记其为对照组，而样本6-10我们标记为实验组。那我们理应可以得到以下数据：
- 对照组：样本1-5，以及其表达数据。
- 实验组：样本6-10，以及其表达数据。

此时我们可以归纳为，这是一个分组实验，每个实验下有5个重复（此为生物学重复）。而每个实验为一个不同的组，这便是样本信息表sample2需要的数据：组和样本名，以便得到实验组之间的表达差异（而非单个样本之间的比较差异）。

什么叫生物学重复?
生物学重复，是指在实验设计中，从不同的生物个体或样品中进行独立的实验测量，以评估样本间的自然生物学变异。其和技术重复不同，对同一个样本进行多次测量，关注的是实验过程中的误差。例如，使用同一份 RNA 样品进行多次测序或 qPCR，旨在评估实验过程中的技术误差。
重点在于是对多个样本重复实验（生物学重复）还是单个样本多次过程重复（技术重复）。
DESeq2和edgeR的选择
- DESeq2
  由于DESeq2软件特性，在生物学重复较多时，DESeq2 更加保守，可能会给出较少的差异表达基因，但结果的可信度更高。
edgeR
而edgeR在生物学重复更少时倾向于更敏感地检测差异表达基因。

因此我们团队在有生物学重复时都会选择DESeq2，而在无生物学重复时会选择edgeR，请根据自己的需要选择适合的方法（METHOD）。

完成上述的步骤后，我们得到了数个matrix矩阵文件。

到该步骤时，我们已经完成了需要在linux服务器上完成的步骤，我们进行了基因比对，定量表达和标准化，差异分析三个大步骤。得到了数个涵盖基因差异的表达矩阵，接下来我们需要将差异分析和定量标准化后产生的文件拷贝到自己的电脑上，进行数据的处理和作图可视化自己的结果。

5 个人计算机上使用R的后续处理

将上述结果拷贝到个人计算机后，还需要使用R语言相关的处理进行文本删减处理，然后对其进行可视化处理。

5.1 IDE

IDE 是“集成开发环境”（Integrated Development Environment）的缩写。它是一个软件应用程序，提供了一系列用于软件开发的工具和功能，以帮助开发者编写、调试和测试代码。R语言的IDE我们这次选择使用Rstudio。

5.2 下载Rstudio

首先需要先下载R到个人计算机，点此跳转至R语言清华镜像站下载网址R
如果是windows系统请点击R for windows
如果是MacOS系统请点击R for macOS
随后请在点此进入Rstudio下载，还是同理安装并打开。

5.3 前期准备

我们需要安装一些必需的库，让我们进行数据分析，还是按照linux版的R的下载指令进行安装

package.install("BiocManager")
BiocManager::install("dplyr")
BiocManager::install("tidyr")
BiocManager::install("tidyverse")
BiocManager::install("clusterProfilter")
#以上四个包是我们进行表达数据的提取与挖掘需要用到的一些函数包，包含了使用的表格，数据提取，表格编辑等功能

5.4 进行分析与可视化

首先我们需要对这些矩阵进行处理，突出我们需要的重要的信息，挖掘异常表达的基因，找到我们需要的基因。

5.4.1 表格导入与前期处理

在R语言中，并不自带linux中的|这样的管道，极其不便于我们连续的操作，tidyverse中提供了%>%作为管道功能，将这一步骤的运行结果传递给下一步进行。
如何更方便的导入tibble表格？点击右上的Import Dataset，并点击from text(readr)，在FILE中提供路径，将Delimiter改为制表符（tab），先观察是否可以列名对应上（即是否存在少列名的情况，如果少列名请先在原matrix里面在缺失那列中间打一个空格），即可形成一张tibble表格，列名为导入表格的第一行。

使用tidyverse包中附带的readr函数包进行表格读取的操作
#准备TMM标准化后的信息表
sample_info <-read.table("~/genes.TPM.not_cross_norm.TMM_info.txt",header = T, row.names = 1)

#准备表达矩阵，导入read counts数据
gene_counts <- read.table("~/genes.counts.matrix",header = T, row.names = 1)

#导入标准化表达矩阵
gene_exp <- read.table("~/genes.TMM.EXPR.matrix",header = T, row.names = 1)

#此处额外使用delim函数额外读取一个TMM表达矩阵，使用read_delim生成的表格并非R语言中的数据框dataframe格式，而是tibble表格格式，可以有更好的行列名以及行列转换的操作，该表格用于后续的画图操作，前三张表用于后续表格处理，以根据自己的不同需求画出不同的数据图
genes_TMM_EXPR <- read_delim("~/Downloads/result/3.Merge_result/genes.TMM.EXPR.matrix", 
                             delim = "\t", escape_double = FALSE, 
                             trim_ws = TRUE) %>%
                dplyr::rename(id=...1)  #将第一列名改为id，此处应该先观察自己的表格，第一列的列名为什么，如果没有列名则tibble默认列名为...1，需要改为id方便后续表格统计

#差异分析结果导入，同样使用delim函数
DE_result <- read_delim("~/genes.counts.matrix.0minG_vs_90minG.DESeq2.DE_results", 
                        delim = "\t", escape_double = FALSE, 
                        trim_ws = TRUE)
DE_result <- DE_result %>%
        dplyr::rename(id=...1)    #将第一列名改为id

该步骤只是导入了原表格数据，我们还需要进行一些处理,将复杂的DE差异分析表格简介明了的表现差异信息

#差异表达结果处理
DE_result <- 
  # 添加上下调信息
  mutate(DE_result, direction = if_else(
    padj > 0.05, 'ns', if_else(        
      abs(log2FoldChange) < 1, 'ns', if_else(
        log2FoldChange >= 1, 'up', 'down') #将log2foldchange的绝对值小于1且padj大于0.05的定义为结果不显著，log2foldchange绝对值大于1的定义为上调或下调
    ))) %>%
  # 关联基因信息，将样品id和基因id关联
  left_join(gene_info, by = c('id' = 'GID')) %>%  
  # 关联表达量信息将表达矩阵与DE结果合并，合并以id列为依据，同样品id的行合并
  left_join(rownames_to_column(gene_exp, var = 'id'), by = 'id') %>%
  # 按 log2FoldChange 绝对值降序排列
  arrange(desc(abs(log2FoldChange)))
#lef_jion函数指将某个表与原表进行合并，合并的依据由两表提供，可以是相同的列，根据该列中每行对应的名字合并。
#mutate表示添加列，添加的列的信息符合direction里的要求

经过了上述这个表格的处理，我们可以得知基因的上调与下调，并选择关注表达量差异更大（log2foldchange）且误差更小（padj）的基因数据。
经过以上步骤后生成了五个表格，可以通过print命令将其打印出来。

如果以后我们需要快速读取这五个表格，我们还需要重新运行一次上面代码吗？
我们可以通过将其储存为rdata的方式在r中快速导入表格信息使用save命令保存为rdata

save(gene_counts, gene_exp, genes_TMM_EXPR,
   sample_info,DE_result,
   file = '~/*.rdata')

此时我们就可以作出火山图和热图了
如果需要作图，我们需要加载额外的绘图库

5.4.2 火山图

#导入前文处理好的文本
load('~/*.rdata')
BiocManager::install('EnhancedVolcano')   #简单的火山图绘制软件
#作图（火山图）
library(EnhancedVolcano)
EnhancedVolcano(DE_result,
                x='log2FoldChange',      #X轴
               y='padj' ,               #Y轴
               lab = DE_result$id,    #选择散点
             FCcutoff = 1,           #log2切除线
            pCutoff = 0.01)         #p值切除线

此时我们就已经生成了一个简单的热图，自己可以尝试一下。
火山图很简单对吧？原理就是读取DE表格中的Log2foldchange作为X轴，左右差距体现基因的表达量差异，Y轴则是p值的-log10，上下差距体现基因可信度，越高越可信。
由此可以直观的看出哪些基因表达差异又高又准确。

5.4.3 热图

热图的绘制需要将前面的数个合为一个表，我们需要得到的是样本与基因之间的表达关系，因此需要数表关联，将基因，样本id，表达量，p值合并进一个表，并且列名与行名分别是样本与基因。

#三表关联,用于制作热图
TNP <- dplyr::select(DE_result,id,log2FoldChange,pvalue,padj) %>%    #选择出DE列，筛选出id，log2，pvalue，padj列
  mutate(direction=if_else(abs(log2FoldChange)<1 | padj>0.01,'ns',if_else(log2FoldChange >= 1 ,'up','down'))) %>% #设置不显著 上调和下调
  left_join(genes_TMM_EXPR, by = c('id'='id'))       #合并TMM表格和DE表格，用id列对照
#统计差异基因,该步骤可以反映有多少上调或下调的基因，如果基因过多可以将log2Foldchange的值调大，如if_else(abs(log2FoldChange)<2 | padj>0.01,'ns'，如基因过少则可以将值调小，p值同理
group_by(TNP,direction) %>% 
  summarise(count=n())

获得了一个关联表后，我们就可以进行热图绘制，本次使用更为简单的pheatmap绘图包进行绘图。

BiocManager::install("pheatmap")
library(pheatmap)
temp <- arrange(TNP,desc(abs(log2FoldChange))) %>%              #选择表格，以及用于制作热值的列
  dplyr::slice(1:15) %>%                                        #选择基因
  dplyr::select(-log2FoldChange,-pvalue, -padj,-direction) %>%    #选择计算的列
  column_to_rownames(var='id')                                 #行名和列名
pheatmap(log2(temp+1))                                        #出图

这是一个示例图

## # A tibble: 3 × 2
##   direction count
##   <chr>     <int>
## 1 down         13
## 2 ns           22
## 3 up           20

5.5 ggplot2

绘图软件是有非常多的选择，如果需要一个自定义属性更强大，绘图能力也更强大的程序包，推荐使用ggplot2。

ggplot2的强在哪里？ ggplot2的作图逻辑为先建立底层xy轴图层，后将按照图层的顺序添加散点，柱状图等，因此ggplot2可以绘画生信方面几乎所有图案（但是流程和图层会更加繁琐），并且可以自定义主题颜色，图点形状，大小，渐变等。
如何使用ggplot开始绘图？以火山图为例 请记得安装ggplot2 BiocManager::install("ggplot2")

library(ggplot2)
ggplot(DE_result, aes(x = log2FoldChange, y = -log10(padj))) +  #此处构建xy轴，选定列
  geom_point(aes(color = ifelse(log2FoldChange < -1 & padj<0.05 , "blue", ifelse(log2FoldChange > 1 & padj <0.05, "red", "gray"))), 
             alpha = 0.7, size = 2) +  #根据log2FoldChange 值设置颜色，同时设定散点透明度与大小
  scale_color_manual(
    name = "significant" ,    #点意义命名
    values = c("blue" = "blue", "red" = "red", "gray" = "gray")) +  #定义颜色
  geom_text(
    data = subset(DE_result, abs(log2FoldChange) > 10 & padj < 0.05),  #根据阈值筛选外面的点
    aes(label = id), #点名取自DE_result的id列
    vjust = 1, hjust = 1, size = 3, check_overlap = TRUE) +  #控制标签位置和大小
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "black") +  #添加p值参考线
  geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "black") +  #添加FoldChange阈值参考线
  theme_classic() +  #使用经典绘图，没有背景图
  labs(title = "sample", x = "Log2 Fold Change", y = "p-value") + #标题名，x轴名，y轴名
  theme(plot.title = element_text(hjust = 0.5))  #使得标题居中

是不是ggplot图看起来更好看一些呢，他提供大量自定义操作，可以更换主题颜色等，使得你的图更美观。

热图可以参考这个进行绘图

#加载ggplot2
library(ggplot2)
#将之前处理好的矩阵转换为长格式
temp_long <- temp %>%
  rownames_to_column("gene") %>%
  tidyr::pivot_longer(cols = -gene,names_to = "sample",values_to = "Expression")
temp_long$Expression <- as.numeric(temp_long$Expression)
#使用ggplot2绘制热图
ggplot(temp_long, aes(x = sample, y = gene, fill = Expression)) +   #设置xy轴
  geom_tile(color="black") +  #创建热图
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 100) +  # 设置渐变色
  theme_minimal() +  #使用简约主题
  labs(title = "Heatmap of sample", x = "Samples", y = "Genes", fill = "Expression") #设置标题名字，XY轴名

此处提供了ggplot基础作图用法，再后续中还会有使用到。

5.6 关于一些R包的解释

5.6.1 Rsubread

featurecounts就是Rsubread的核心函数，其函数包主要用于RNA-seq的比对，单端-双端reads的计数，并且支持多线程并行操作，速度得到了极大提升，目前主要是使用与定量分析。

5.6.2 edgeR

它基于负二项分布（Negative Binomial Distribution）建模，非常适合处理有生物学重复的数据集，尤其是小样本数据集。edgeR 能有效分析基因表达、ChIP-Seq、miRNA 数据以及其他涉及计数数据的实验。
原理是先以DGElist对象储存测序数据，而后使用”calcNormFactors”函数进行TMM归一标准化，随后使用estimateDisp估计样本离散程度使用exactTest或glmFit进行差异表达检验，得到每个基因的p值和倍数变化。最后得出结论。因此他更适合少量的生物学重复的计数。

5.6.3 DESeq2

不同于edgeR，DESeq2使用收缩的方式进行技术统计，使用经验贝叶斯方法估计基因的离散度，并通过方差缩减来稳定基因的倍数变化估计，尤其是对低表达基因，减少不可靠的波动。
因此DESeq2更适合处理生物学重复较多的实验，特别是由于其稳健的方差缩减方法，能够有效减少假阳性结果。

5.6.4 DESeq2 vs edgeR

这个时候大家心里一定有一个疑问，这个时候我们如果在生物学重复处于不多不少的数量的时候，我到底选择哪一个软件进行差异分析呢，我的建议是都可以，取决于你想得到什么样的结果，以什么样的统计方法获得什么样的数据对于事实上具有显著差异的基因而言，不会因为软件的更换，而变得不显著，显著基因始终是显著基因，而整体数据的呈现效果则取决于你使用的软件。
点此返回

生信分析教程转录组分析简要说明

向旸（遗传多样性与进化团队）

2024