小鼠转录组测序上游数据分析

一、流程概括

• RNA-seq原始数据的质量评估

• 低质量测序数据的过滤和清除

• 过滤后数据的质量评估

• reads比对基因组和转录组

• 计数和定量

二、准备工作

• 注释文件和基因组文件的准备

• 相关软件的安装

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

• 可以从NCBI或者Ensembl网站下载，以Ensembl网站提供的基因组为例，可以选择Mus_musculus.GRCm39.dna.primary_assembly.fa 或者早期的GRCm38版本
• 注释文件也可以在同样的网站下载，以Ensembl网站提供的注释文件为例，可以选择Mus_musculus.GRCm39.109.gtf 或者早期的GRCm38对应版本

• 鉴于测序物种是小鼠，也可以选择GENCODE网站完成基因组文件和注释文件的下载

  
  基因组
  注释文件
  

  最重要的是，用于比对的基因组文件和注释文件的版本要保持一致。

2.软件安装

使用conda进行软件的安装。

• 质控：fastqc, multiqc, trim-galore

• 比对：hisat2, samtools

• 计数：featurecounts, bedtools, salmon

首先创建名为rna_seq的conda环境

conda create -n rna_seq

• 配置conda，添加清华镜像源

conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/
conda config --set show_channel_urls yes

软件安装：conda install <software>会自动安装软件和对应的依赖环境
值得注意的是需要在rna_seq的环境下安装以上软件，激活rna_seq环境的代码如下：

source activate rna_seq

三、质量报告的生成

fastq质量报告

使用如下命令来生成质量报告，每个fastq文件会获得一个质量分析报告，来描述RNA-seq的测序质量。

fastqc -o output_dir *.fastq.gz

multiqc可以对所有生成的fastqc报告文件进行总结并汇总到一个报告文件中，以便更直观地展示。命令如下：

multiqc <analysis directory>

这一步的目的是根据质量报告结果进行低质量序列的切割。

四、数据处理

数据处理内容：过滤低质量reads，去除测序接头

处理软件：数据处理可以使用多款软件，trim_galore在各文献中表现良好。

1.trim_galore 的使用方法

trim_galore可以处理illumina，nextera3，smallRNA测序平台的双端和单端数据，包括去除adapter和低质量reads。

trim_galore的参数：

trim_galore [options] <filename>
--quality<int>  #设定phred quality阈值。默认20（99%的read质量），如果测序深度较深，可以设定25
--phred33       #设定记分方式，代表Q+33=ASCII码的方式来记分方式。这是默认值。
--paired          # 对于双端结果，一对reads中若一个read因为质量或其他原因被抛弃，则对应的另一个read也抛弃。
--output_dir   #输出目录，需确保路径存在并可以访问
--length        #设定长度阈值，小于此长度会被抛弃。
--strency     #设定可以忍受的前后adapter重叠的碱基数，默认是1。简单来说就是最多能允许末端残留多少个接头序列的碱基
-e<ERROR rate>  #设定默认质量控制数，默认是0.1，即ERROR rate大于10%的read 会被舍弃，如果添加来--paired参数则会舍弃一对reads
<filename>  #如果是采用illumina双端测序的测序文件，应该同时输入两个文件。

命令如下：

trim_galore -q 25 --phred33 --stringency 3 --length 36 --clip_R1 20 --clip_R2 20 --paired ./seq_R1.fastq.gz ./seq_R2.fastq.gz --gzip -o ./clean_data

其中--clip_R1和--clip_R2参数的具体值是根据测序质量报告确定的，通常会去掉5'端前面15-20bp测序质量不稳定的部分。

2. 处理后数据的质量分析

对过滤后的文件进行质量分析，同样使用fastqc和multiqc两个软件进行质量分析。
观察到总reads数减少和总体reads质量变高，测序adapter也被去除。更具体的数据处理情况可以在seq_trimming_report.txt中查看。

SUMMARISING RUN PARAMETERS
==========================
Input filename: ./G1_1_combined_R1.fastq.gz
Trimming mode: paired-end
Trim Galore version: 0.6.2
Cutadapt version: 3.4
Number of cores used for trimming: 1
Quality Phred score cutoff: 25
Quality encoding type selected: ASCII+33
Adapter sequence: 'CTGTCTCTTATA' (Nextera Transposase sequence; auto-detected)
Maximum trimming error rate: 0.1 (default)
Minimum required adapter overlap (stringency): 3 bp
Minimum required sequence length for both reads before a sequence pair gets removed: 36 bp
All Read 1 sequences will be trimmed by 20 bp from their 5' end to avoid poor qualities or biases
All Read 2 sequences will be trimmed by 20 bp from their 5' end to avoid poor qualities or biases (e.g. M-bias for BS-Seq applications)
Output file will be GZIP compressed

This is cutadapt 3.4 with Python 3.6.7
Command line parameters: -j 1 -e 0.1 -q 25 -O 3 -a CTGTCTCTTATA ./G1_1_combined_R1.fastq.gz
Processing reads on 1 core in single-end mode ...
Finished in 706.34 s (23 us/read; 2.59 M reads/minute).

=== Summary ===

Total reads processed:              30,512,439
Reads with adapters:                 2,554,479 (8.4%)
Reads written (passing filters):    30,512,439 (100.0%)

Total basepairs processed: 4,576,865,850 bp
Quality-trimmed:              20,143,870 bp (0.4%)
Total written (filtered):  4,482,066,612 bp (97.9%)

五、比对

概况：使用处理后的fastq文件与基因组/转录组比对，确定在转录组/基因组中的关系。转录组和基因组的比对采取的方案不同，分别是ungapped alignment to transcriptome和Gapped alignment to genome。
软件：hisat2和STAR在比对回帖上都有比较好的表现。有文献显示，hisat2纳伪较少弃真较多，但是速度比较快。STAR就比对而言综合质量比较好，在长短reads回帖上都有良好发挥（Sahraeian SME, Mohiyuddin M, Sebra R, et al. Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis. Nat Commun. 2017; 8(1): 59）。
由于hisat2的速度优势，选择hisat2作为本次比对的软件，在比对之前首先要进行索引文件的获取/制作。

1. 索引文件的获取和生成

• 不同的比对软件构建索引方式不同，所用的索引也不尽相同
• 下载hisat2基因组索引：https://daehwankimlab.github.io/hisat2/download/
• 本次分析采用的是GRCm39参考基因组，网页中没有现成的index文件，使用下面的命令构建索引文件

hisat2_extract_exons.py /ref/annotation/Mus_musculus.GRCm39.109.gtf > /ref/annotation/exons_mouse.txt
hisat2_extract_splice_sites.py /ref/annotation/Mus_musculus.GRCm39.109.gtf > /ref/annotation/splices_sites_mouse.txt
hisat2-build -p 8 --ss /ref/annotation/splices_sites_mouse.txt --exon /ref/annotation/exons_mouse.txt /ref/genome/Mus_musculus.GRCm39.dna.primary_assembly.fa GRCm39

• 索引文件的格式如下，由多个文件构成，要保证索引文件的格式和名称部分一致

hisat2索引文件列表

2. hisat2的比对

使用hisat2比对

hisat2 -p 12 -x /ref/index/GRCm39 -1 ./clean_data/seq_val_1.fq.gz -2 ./clean_data/seq_val_2.fq.gz -S ./align_data/seq.hisat2.sam

参数说明：

-p #多线程数
-x #参考基因组索引文件的目录和前缀
-1 #双端测序中的一端测序文件
-2 #双端测序中的另一端测序文件
-S #输出的sam文件

说明：在比对过程中，hisat2会自动将双端测序匹配同一reads并在基因组中比对，最后两个双端测序生成一个sam文件。比对过程需要消耗大量时间和电脑运行速度和硬盘存储空间，比对完成会生成一个报告。

samtools 软件进行格式转换

sam格式文件由于体量过大，一般都是使用bam文件来进行存储，命令如下：

samtools view -S seq.sam -b > seq.bam  #文件格式转换
samtools sort  -@ 16 seq.bam -o seq_sorted.bam  #将bam文件排序
samtools index seq_sorted.bam  #对排序后的bam文件索引

至此，一个比对到基因组的RNA-seq文件构建完成。

3.对比对后的bam文件进行质量评估

samtools flagstate ：统计bam文件中比对flag信息，然后输出比对结果。

samtools flagstate seq_sorted.bam > seq_sorted.flagstate

六、定量

featurecounts定量

featureCounts -T 16 -t exon -g gene_id -a <gencode.gtf> -o seq_gene_counts.txt <seq.bam>

参数：

-g # 注释文件中提取对Meta-feature 默认是gene_id
-t # 提取注释文件中的Meta-feature 默认是 exon
-p #参数是针对paired-end 数据
-a #输入GTF/GFF 注释文件
-o #输出文件

salmon定量

与hisat2不同，salmon不经过比对计数步骤而直接对基因进行定量，如果不研究新转录本，用salmon方法可以更快更方便得到基因的定量信息。
先下载参考转录本序列cDNA.fa文件，在ensembl官网选择相应文件Index of /pub/release-109/fasta/mus_musculus/cdna
salmon软件推荐根据基因组和转录组参考序列构建decoys文件之后再建立索引，构建decoys的脚本可以在GitHub上下载。构建decoys的步骤时间较长，不要认为是流程卡住而随意中断。
使用salmon index 建立索引文件：

generateDecoyTranscriptome.sh -a /ref/annotation/Mus_musculus.GRCm39.109.gtf -g /ref/genome/Mus_musculus.GRCm39.dna.primary_assembly.fa -t /ref/transcriptome/Mus_musculus.GRCm39.cdna.all.fa -o /ref/index/mm_decoy
salmon index -p 16 -t /ref/transcriptome/Mus_musculus.GRCm39.cdna.all.fa -i /ref/index/salmon --decoys /ref/index/mm_decoy/decoys.txt

使用salmon quant命令对clean fastq文件直接进行基因定量：

salmon quant -i /ref/index/salmon -l A -1 ./clean_data/seq_val_1.fq.gz -2 ./clean_data/seq_R2_val_2.fq.gz -p 16 --output ./salmon/seq_quantity --validateMappings

salmon定量参数

运行结果存放在salmon文件夹下，里面的quant.sf即为下游分析所需要的文件。

至此，转录组测序上游数据分析全部结束。