mRNA-seq分析
mRNA-seq Analysis   

仪器设备

1. 服务器 (型号：I620-G20；系统：centos7；CPU：Xeon(R) CPU E5-2620 v4 @2.10GHz，3核)，用于运行相关命令
2. 普通个人电脑 (需要已经安装好R，能调用个人终端。本教程是使用Mac系统进行操作的。)
   

软件

1. FastQC (v0.11.8)
2. mulitiQC (v1.6)
3. flexbar (v 3.3.0)
4. STAR (STAR_2.6.0a)
5. samtools (v1.9)
6. HTSeq (0.11.0)
7. R (v3.5.1)
   

实验步骤

注：本教程是基于已经在服务器上和个人电脑上安装好相关软件进行的。如果安装出现问题，请参考以下链接：
FastQC：https://www.bioinformatics.babraham.ac.uk/projects/download.html
mulitiQC：https://multiqc.info/
flexbar：https://github.com/seqan/flexbar
STAR：https://github.com/alexdobin/STAR
Samtools：http://www.htslib.org/download/
HTSeq：https://htseq.readthedocs.io/en/release_0.10.0/

1. 数据获得
   本教程以我们实验室2018的eLIFE文章 (Epigenetic drift of H3K27me3 in aging links glycolysis to healthy longevity in Drosophila) 中的部分数据作为测试数据(Ma等, 2018)。
   注：一般已发表文章所包含的数据可以在NCBI (SRA、GEO等)、EMBL-EBI等相关数据库获得。本文中的原始数据 (raw data) 可以在eLIFE文章 (Ma等, 2018) 中获取。
2. 通过FastQC对数据进行质量评估
   $ fastqc所在目录/fastqc -t 4 -o 输出文件夹 输入文件
   
   在文件夹可以看到每个输入文件有两个输出结果，一个结果集合压缩包和一个网页型结果。一般通过查看网页结果即可知道基本的测序质量。但是具有多个样品逐一打开过于麻烦，所以可以通过multiQC将所有的结果综合到一起查看，输入命令：
   $ multiqc 目标文件夹
   输出结果是网页形式，结果如图1所示：
   
   
   图1. MultiQC输出结果
   
   通过结果得知，测序样品结果稳定，测序质量很高。除了用multiQC查看多个QC结果以外，我们还用名为fastqc_summary.py的一个python脚本看每个样本的reads数量，GC含量，Q20，Q30的比例。脚本如图2所示：
   
   
   图2. fastqc_summary.py脚本内容
   
   这段脚本的逻辑为：
   
   1)

   用python的zipfile模块打开zip文件，读取xx_data.txt的数据。
   2)

   读取每一行的数据，用正则表达式模块re，找到目标行。
   3)

   根据分隔符对每一行进行分割，进行赋值。
   4)

   由于只需要读取到>>Per base sequence quality pass这一部分，所以设置一个>>END_MODULE的计数器，数量超过2，就停止。
   
   运行命令为：
   
   $ python3 ./fastqc_summary.py 目标文件
   
   运行后会输出一个txt文本，结果如图3所示：
   
   
   图3. 测序Q30结果
   
   结果表明Q20指标在1.0，Q30指标在0.97，可以说测序质量非常好了！
3. 对测序文件进行比对
   
   3.1

   flexbar
   在比对文件之前需要先要去掉adapter，得到不含有adapter的fastq文件。我们采用的是flexbar程序，运行代码如下：
   
   $ FLEXBAR -t 输出文件名 -qf i1.8 -r 输入文件 -a Adapters.fa
   
   
   3.2

   STAR
   获得trimmed数据后，就要与参考文件进行比对。目前大家所认同的比对算法有很多，比如Bowtie2，BWA，TopHat，HISAT2和STAR等等。我们采用的是综合指标 (Sahraeian等，2017) 比较高的STAR进行比对。
   注：由于STAR的运算需要较高的配置，当使用个人电脑进行学习分析流程的时候建议使用TopHat进行比对，可参考https://www.jianshu.com/p/0abf87b3aca0进行操作。如果是基因组较小的物种 (果蝇、线虫等) 也可以使用个人电脑运行STAR (测试过的笔记本电脑配置：MacBook Pro，2.3GHz Intel Core i5处理器，8GB内存)。
   使用STAR，首先需要对参考基因组建立索引，代码如下：
   
   $ STAR --runThreadN 线程数 --runMode genomeGenerate --genomeDir 存放参考基因的目录 --genomeFastaFiles 参考基因组的fastq文件 --sjdbGTFfile 参考基因组的gtf文件 --sjdbOverhang 索引长度
   
   建立好索引就可以对数据进行比对了，比对命令如下：
   
   $ STAR --genomeDir 存放参考基因的目录 --readFilesIn 输入文件 --runThreadN 线程数 --outSAMtype SAM --outFileNamePrefix 输出文件前缀
   
   STAR运行结束后会生成一个log.final.out文件，打开该文件可以看到比对结果如何。图4是本篇protocol中用到的示例数据的比对结果，特异比对率在92.9%，说明比对结果很好，样品没有污染，可信度高。
   注：一般情况下，uniquely mapped reads %在90%以上就认为样品良好，没有污染。如果比对比例低于70%，需要慎重对待。可以采取的办法是，只保留比对上的数据，如果数据量满足需求 (比如重测序深度达10X，则可以继续使用)；如果不满足要求，则需要重新准备样本。而对于很多医学样品可能是来源于感染的组织/细胞，那么它们的比对到其中一个基因组的效率低是正常的。
   
   
   图4. STAR比对结果
   
   
   3.3

   Samtools
   注：Samtools的使用主要是为了生成可以在igv中可视化的bw文件，后续htseq-count可以不使用此部分结果。
   在STAR中，它是可以直接将sam文件转换成bam文件并进行排序的。为了后续的分析的方便，我们是采取samtools进行这一操作。
   我们首先使用samtools将STAR比对后生成的sam文件转换成bam文件，命令如下：
   
   $ SAMTOOL view -b -S 输入sam文件 -o 输出bam文件
   注：Samtools的新版本已经可以直接不用转换的sam文件进行Samtools sort的排序操作。
   
   然后对获得的bam文件进行排序，命令如下：
   
   $ SAMTOOL sort -o 输出排序后文件 输入bam文件
   
   最后，我们将已经排序的bam文件生成index，命令如下：
   
   $ SAMTOOL index 输出文件
   
   这时获得的bam文件可以直接在igv中进行可视化等其他操作。
   3.4

   HTseq
   在RNA-seq中，我们最关心的是基因的表达量，所以我们通过HTseq程序获得每个基因的reads数，通过reads数可以进行后面的差异表达分析。当然根据问题需要，还可以通过cufflinks等程序选择计算TPM，RPKM等进行下游分析。使用HTseq的命令是：
   
   $ htseq-count 输入比对后的sam文件GTF参考文件 -m intersection-strict -s no > 输出counts文件
   注：这里输入的STAR比对输出的sam文件，对于双端测序是必须输入排序后的sam/bam文件 (htseq-count中含有排序命令，详情请在终端输入htseq-count -h查阅具体使用方法)，对于单端测序排序与否可以忽略。
   
   
4. 差异表达分析
   HTseq的输出文件可以直接与R包DESeq2对接，得到差异表达分析数据。将R的工作环境设置到含有HTseq输出的counts文件下，运行如下脚本 (图5)：
   
   
   图5. DESeq2差异表达分析R脚本
   
   运行脚本后我们可以获得一个差异表达的表格，标准化了的counts表达矩阵，以及一张PCA图 (图6)。
   
   
   图6. PCA图
   
   
5. 下游分析
   到目前为止，我们已经获得了实验组和对照的差异表达数据集。在这个数据集中包含了差异表达倍数、P值、本底表达量等信息，后面可以基于这组数据进行下游分析。比如对筛选出的差异基因做热图展示、GO分析、KEGG分析等。
   

参考文献

1. Ma, Z., Wang, H., Cai, Y., Wang, H., Niu, K., Wu, X., Ma, H., Yang, Y., Tong, W., Liu, F., Liu, Z., Zhang, Y., Liu, R., Zhu, Z. J. and Liu, N. (2018). Epigenetic drift of H3K27me3 in aging links glycolysis to healthy longevity in Drosophila. Elife 7: e 35368.
2. Sahraeian, S. M. E., Mohiyuddin, M., Sebra, R., Tilgner, H., Afshar, P. T., Au, K. F., Bani Asadi, N., Gerstein, M. B., Wong, W. H., Snyder, M. P., Schadt, E. and Lam, H. Y. K. (2017). Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis. Nat Commun 8(1): 59.