RNA-seq：转录组数据分析处理（上）

RNA-seq：转录组数据分析处理（上）

一、流程概括

• RNA-seq的原始数据（raw data）的质量评估
• linux环境和R语言环境
• raw data的过滤和清除不可信数据（clean reads）
• reads回帖基因组和转录组（alignment）
• 计数（count ）
• 基因差异分析（Gene DE）
• 数据的下游分析

二、准备工作

• 学习illumina公司测序原理
• 测序得到的fastq文件
• 注释文件和基因组文件的准备

1. fastq测序文件

在illumina的测序文件中，采用双端测序（paired-end），一个样本得到的是seq_1.fastq.gz和seq_2.fastq.gz两个文件，每个文件存放一段测序文件。在illumina的测序的cDNA短链被修饰为以下形式（图源见水印）：

两端的序列是保护碱基（terminal sequence）、接头序列（adapter）、索引序列（index）、引物结合位点（Primer Binding Site）：其中 adapter是和flowcell上的接头互补配对结合的；index是一段特异序列，加入index是为了提高illumina测序仪的使用率，因为同一个泳道可能会测序多个样品，样品间的区分就是通过index区分。参考：illumina 双端测序（pair end）、双端测序中read1和read2的关系。

在illumina公司测得的序列文件经过处理以fastq文件协议存储为*.fastq格式文件。在fastq文件中每4行存储一个read。
第一行：以@开头接ReadID和其他信息，分别介绍了
第二行：read测序信息
第三行：规定必须以“+”开头，后面跟着可选的ID标识符和可选的描述内容，如果“+”后面有内容，该内容必须与第一行“@”后的内容相同
第四行：每个碱基的质量得分。记分方法是利用ERROR P经过对数和运算分为40个级别分别与ASCII码的第33号!和第73号I对应。用ASCII码表示碱基质量是为了减少文件空间占据和防止移码导致的数据损失。fastq文件预览如下：

2.注释文件和基因组文件的获取

• 基因组获取方式：可以从NCBI、NCSC、Ensembl网站或者检索关键词“hg38 ftp UCSC” 人类基因组hg38.fa.gz大概是938MB左右。文件获取可以点击网站下载。可以通过云盘的离线下载来加速下载进程
• 基因组的选择：以Ensembl网站提供的基因组为例，比对用基因组应该选择Homo_sapiens.GRCh38.dna.primary_assembly.fa
• Ensembl基因组的不同版本详见README和高通量测序数据处理学习记录（零）：NGS分析如何选择合适的参考基因组和注释文件

三、软件安装

• 安装方式：软件安装可以通过例如apt-get、miniconda等方式来安装。由于miniconda的便捷行，使用conda进行如下软件的安装。
• 软件列举
  质控：fastqc ,multiqc , trimmomatic, cutadapt, trim-galore
  比对：star , hisat2 , tophat , bowtie2 , bwa , subread
  计数：htseq , bedtools, salmon, featurecount
• miniconda的安装：

1. 可以通过点击清华大学开源软件站或者检索“清华大学 conda”访问镜像网站（清华镜像站因为服务器在中国访问速度比较快），点击Anoconda界面，选择Miniconda下载安装，windows在安装好需要设置环境变量。
2. linux测试Miniconda的安装：conda -v
3. 创建名为rna的环境变量：conda create -n rna python=2（许多软件依赖python2环境）环境退出：conda deactivate
4. 配置conda，添加镜像源头：输入如下命令（更新：2019年05月06日）

conda config --add channels https://mirrors.cloud.tencent.com/anaconda/pkgs/free/
conda config --add channels https://mirrors.c