itr流程总共包含多少个l2子流程_我要自学生信之生信基础-转录组：分析流程大全解，看这一篇就够了...

【写在前面】本文旨在通过一篇文献带大家走一遍RNA-seq的流程。对于整个转录组有个宏观的认识，如果有其他文章写的又全又好的，欢迎私信。说明我这篇文章就没有什么作用了，可以退出历史舞台了。

在文章开始之前，如果你对转录组、对生信没有任何的概念，可以先看下我的这篇文章

生物信息学入门必看的85个名词：每个都要记牢

【参考文献】Transcriptome analysis of an apple (Malus × domestica) yellow fruit somatic mutation identifies a gene network module highly associated with anthocyanin and epigenetic regulation

影响因子：5.83

发表时间：2015年9月

【实验背景】作者发现一块地里面两棵苹果上的果皮颜色是不同的，作者想探究一下是什么引起了果皮颜色的变化，作者通过文献查找发现果皮颜色常常与花青素有关系，于是作者选用了苹果的果皮，把红苹果标注为KID，黄苹果标注为BLO，分别取4个时期，S1-S4，进行了三个生物学重复，一共24个sample。首先进行了花青素含量的测定，发现果然红色苹果的花青素比较高。为了探到底是什么影响了花青素的变化，作者进行了RNA-seq

【转录组流程】文献中标注了测序后数据的来源：SRP062637 我也跟着跑了一遍流程，总结如下

这个过程是一个非常常规的过程，并没有太多新奇的地方。这是一篇典型的有参转录组的文章，所以并不需要进行组装。下面我来具体蜻蜓点水的科普下这个流程中的问题。

【原始数据：raw data】

首先我们需要探究这个数据是怎么来的，大家都会说测序得到的，是的，这个是二代测序的数据，而且是单端测序，我讲一下PE(paired-end)测序，这里有一个概念很多人都一直很模糊，目前通常双端测序测380-400bp是没有问题的，你会不会产生疑惑了一次不是只有150bp吗？这是因为双端测序所以可以测300bp，但中间的80-100bp在这一个reads里面是没有测出来的，但是不要忘记我们测的序列基本是30X的，reads数量相当的多，所以并不会影响mapping结果，这个reads没有测到，下一个reads总有测到的，而且还不少。

那么测序之前的准备有哪些呢？主要是建库，详情如下：

RNA-seq：测序原理之文库构建

下面我们开始进行测序，那么测序的原理以及流程是怎么样的呢？可以参考下面的文章

https://zhuanlan.zhihu.com/p/190757472​zhuanlan.zhihu.com

为了更形象的理解，可以再看下illumina官方的视频。

illumina测序原理简介_哔哩哔哩 (゜-゜)つロ 干杯~-bilibili​b23.tv

现在我们拿到了测序数据，也就是raw data 那么我们得到的数据是什么格式的呢？fastq格式，关于fastq格式我也详细写过一篇文章，参考如下：

我要自学生信之生信基础：FASTA 与 FASTQ

这样我们就可以开始对数据进行处理了

【数据的质控与过滤】

推荐使用现在比较流行的fastp，质控和过滤都可以一键操作。

通过fastp主要做到了以下的几件事情：

1、去接头：分为index+adapter 为什么会有接头呢？

• 文库的插入片段太短，导致片段被测通，还测到了adapter的位置，好像没理解？给你看下我们的测序片段的组成：adapter+index+插入片段+adapter
• 测序引物的位点仍然在接头（index）的左边，这样就会把index也测到

2、去掉N多的测序片段

一个测序片段如果有很多碱基无法识别(就会出现很多的N)，说明测序质量不高或者测序时出现了问题，这个需要注意，即便去掉了N，也要看下是什么原因引起的，不然治标不治本

3、去质量评分低的片段

4、去末端一定数量的片段

随着读长的增加，酶活性下降，荧光强度也在下降，因此测序质量下降是自然趋势

【比对】这一步叫做回帖，比对之前需要先对参考基因组建index，这个作用也比较简单，就类似于像建立图书馆的索引，这样你找书的时候速度就会快很多，仅此而已。将 reads比对到参考基因组，推荐使用 Hisat2或 STAR。STAR需要更大存，运行时间也更长。准确性相差不大。

苹果的参考基因组下载地址：https://github.com/moold/Genome-data-of-Hanfu-apple

目前常见的参考基因组下载地址如下：

Ensembl: ensembl. org/pub(最推荐)

EnsemblGenomes: ensemblgenomes . org/pub/(ensembl的子网站)

NCBI: ncbi .nih. gov/ genomes/

为什么不选用BWA这些对软件呢？

主要是因为可变剪切，RNA-Seq reads由于不包含内含子，所以来自外显子边界处的 reads被重新 mapping回基因组时，会被中间的内含子分开，这种情况叫做splice-alignment。

比对之后的会形成一个sam文件，通过IGV大家可以粗略的看下比对的结果是什么样子的

这里需要解释一下sam文件，这个文件里面会把比对的结果列出来，每条read比对到那个染色体，比对的序列是什么，插入片段都会写，一个reads的描述有三个部分，如果你有10000个reads那么这个文件里面就会有30000个部分。所以这个文件是非常大的，而且也可以看到这个reads的顺序是随机的，没有顺序

举个例子：这里没有选用文章中的数据，但其实无所谓，不影响理解

每一个符合的含义：

具体来讲一下：

得到的sam文件是非常大的，需要压缩成bam文件，然后因为这个文件是无序的，所以还u需要进行sort。

注意：其实这里的比对严格来讲是不规范的，比对的时候最好是要加上注释信息的，因为要考虑可变剪切的问题，所以有位置信息比对效果是更好的

【定量】定量的输入文件是加了index的sort后的bam文件，同时需要加上注释信息，这里为什么要加注释信息呢？主要原因是为了计算基因的表达量，bam文件只有序列比对的情况，而注释信息会标注哪里是基因。比如bam文件里面标注了比对到哪里，注释信息会告诉你这块是什么基因。我们用的是featurecounts来计算表达量，这里的表达量只是说这个基因上比对了多少条reads，这个概念很重要。

【合成成矩阵】现在是有24个样本，没有合并矩阵之前是计算了单个样本的基因表达量，现在合并到一起。行是样本，列是基因。然后进行TPM标准化，这里面还涉及一个标准化的问题，

那么为什么要进行标准化呢？

reads长度是一样的，如果你基因很长，当然很容易比对到这个基因上，样本深度越深也是一样的，测序深度越深，别人1个G，你10个G的数据，reads数量本身就高了。因此一定要标准化。这样的比较才相对公平

很多人会提到FPKM、RPKM，这些你怎么没有说？

这个很多大佬讲的很清晰了，现在比较常用的就是TPM了，另外两个都是错误的。

【差异表达分析】

差异表达分析其实可以在R里面做，软件安装

1.DESeq2两种安装方法：

conda install bioconductor-deseq2
BiocManager::install('DESeq2')

2.edgeR两种安装方法：

conda install bioconductor-edger
BiocManager::install('edgeR')

最终会形成两两比较的结果，筛选出差异基因，筛选条件是log2FoldChange 绝对值>1；pvalue < 0.05，筛选出来了然后使用eggNOG数据库进行基因功能注释

实际上到这里转录组的基本流程已经出来了，当然从这篇文章来讲，前面的分析只是一个准备，下面才刚刚开始。

【差异表达基因分析】

作者第一步做了一个韦恩图

表述了不同样本间的差异基因情况，总计有321019个基因在至少一个样本中表达，感觉并没有什么用，应该是凑图的。

【WGCNA分析】

这个地方可以说一下，我写简单写过WGCNA的流程，算是入门篇，写好进阶后再贴上来。

我要自学生信之生信基础-转录组：WGCNA全流程分析（入门篇）

使用3299条差异基因进行了WGCNA分析，结果如图：

WFCNA分析

图B中可以看到，Pink模块和花青素的含量相关度最高，达到了0.95，具体看下这个Pink模块，一共有34个基因，现在要从34基因里面找到强中介的那几个hub基因，通过mapman进行功能注释，讲基因分为三组，其中group3包含24个基因与花青素相关，其中12个进行过相关的报道。

mapman

mapman

MDMYB10和MdGST与花青素的GS最高。

【实验验证】用qPCR进行实验验证RNA-seq的可靠性，目前还没有人会把RNA-seq作为最可靠的证据，一般都是要进行qPCR实验验证的，而且在实验验证时只有十几个基因，可以使用多个重复，结果更加准确。

差异表达水平

结果发现：MdMYB10 有23倍的差异，MdGST有46倍的差异

【导致差异的原因】很多文章到上一步就停止了，这篇文章的影响因子将近6分，我感觉有一半在这个问题上的讨论，作者提出了几个假设？

• DNA变异：未发现序列差异
• 表观修饰：MdGST未发现甲基化水平差异，MdMYB10的MR3和MR7区域存在差异甲基化，MdMYB10的甲基化会导致在黄苹果中MdMYB10的表达受到抑制

甲基化水平

• 转录调控：MdGST上游有19MYB-，11MYC-IF结合位点，很有可能是转录调控，MdMYB10的表达受到抑制又导致了MdGST的表达受到了抑制

这篇文献就算简略的讲完了，讲的没有想象中那么完美，还是能力不太够，转录组的流程也带着大家走了一遍，当然没有一步一步的带着大家跑，可能还是会有些问题，文献中的数据的下载地址我也在文中写了，具体的代码我没有给大家列出来，希望大家能够根据流程自己尝试做一下，有什么问题欢迎一起讨论。

写这篇文章花费了自己很长的时间，算是对自己生日的礼物。祝自己生日快乐，愿不负韶华，送自己一句俞敏洪老师的话:一个人真正优秀的特质，来自于内心想要变得更加优秀的那种强烈的渴望，和对生命追求的那种火热的激情！与君共勉。