生信马拉松 Day17 转录组RNA-seq-2

今天的主要内容是转录组上游的质控，设计到4个包：fastqc、multiqc、trim_galore、fastp

质控的目的：1.数据质量评估2.过滤低质量

现在大部分数据质量都很好，一般质控得到的结果就是数据质量很不错，然后就直接进行下游的数据比对和定量，不做数据处理，但是质控的内容依然需要掌握

内容1：FastQC

FastQC软件可以对fastq格式的原始数据进行质量统计，评估测序结果，为下一步修剪过滤提供参考

FastQC主页：http://www.bioinformatics.babraham.ac.uk/projects/fastqc/

常见参数

FastQC的常用参数

数据质控-fastqc运行

方法一：直接运行，可以再开一个终端运行多个程序

fastqc -t 6 -o ./ SRR*.fastq.gz

方式二：在命令前后加上nohup &，适用于比较简单的命令

nohup fastqc -t 6 -o ./ SRR*.fastq.gz >qc.log &

方式三：将命令写入sh脚本，使用nohup &运行sh脚本，适用于比较长和复杂的命令

nohup bash qc.sh > qc.log &

示例代码

# 激活conda环境
conda activate rna

fastqc -t 1 -o  ./qc/  ./fastq_25000/SRR1039510_1.fastq.gz

# 多个数据质控
fastqc -t 1 -o  ./qc/  ./fastq_25000/SRR*.fastq.gz

fastqc报告结果为带有fastqc结尾的文件，包括html和zip，html使用浏览器打开

需要关注的FastQC结果内容：

Filename 文件名

Encoding Illumina 1.9 代表是+33

Total Sequences 测序的总reads数

Sequnce length 原始数据的测序长度是一致的，经过过滤的长度会发生变化，报warning但是是正常的

数据量统计方式

Per base sequence quality

以碱基为单位：每个位置为对象上所有read的碱基质量值的箱线图

一共25000个read，第一个位置上的base就有25000个质量值，这个分布情况画成一个箱线图，依次往后推

绿色区域是碱基质量比较高的，红色区域就是出错率>1%，都分布在绿色区域说明出错率比较低，质量高，大量分布在红色就是大量碱基出错率相当高

Per tile sequence quality

Illumina测序，并且数据保留了原始的序列标识符，tile可以裂解为测序仪上的小单元，blue代表质量好，一个好的结果应该是全蓝色的

Per sequence quality scores

查看数据中每条序列质量分布情况，可判断是否有普遍质量较低的质量的序列，quality score distribution在>34位置分布多，就是质量比较高的状态，可以作为筛选方向

Per sequence GC content

以每条序列GC碱基比例为对象

红色和蓝色的拟合线接近，GC含量是不能作为筛选靶标

Sequence length distribution

sequence length没有筛选标准，一般都是63bp左右的长度

Sequence duplicattion levels

重复数为1的含量特别高，代表read差异性比较好

Per base sequence content

以ATGC碱基比例，理论上，ATGCN的含量每个序列循环上应分别相等，且整个测序过程稳定不变，呈水平线，是紧紧缠绕在一起的曲线，这个结果受测序影响

Per base N content

N含量分布图，N碱基比例越小越好

Overrepresented sequence

一般认为没有单个序列占整体的一小部分，发现单个序列在集合中的代表性非常高，要么意味着它具有高度的生物学意义，要么表明文库被污染，或者没有预期的那么多样化

Adapter content

接头含量，横坐标是碱基上的位置，纵坐标是占的百分比，-a参数可以自定义接头序列

内容2：multiqc

使用multiqc整合fastqc结果

multiqc不生产数据，只是数据的搬运工

示例代码

multiqc  ./qc/*.zip  -o ./qc/

内容3：trim_galore

测序得到的原始序列含有接头序列或低质量序列，需要对原始数据进行质控，得到高质量序列（即clean reads）

建议标准为：

（1）接头：去除含接头的reads

（2）序列质量：过滤去除低质量值数据，确保数据质量，因人而异，什么都不知道的情况用软件的默认值就好

（3）N碱基比例：去除含有N（无法确定碱基信息）的比例大于5%的reads

trim_galore：是Cutadapt和FastQC的集合，前者去除接头问题，后者生成报告，因此可以用multiqc在质控后得到整合数据

-j --cores 使用线程数，默认为1

-a --adapter 输出adapter序列，也可以不输入，自动搜索可能性最高的平台对应的adapter，或者可以直接输入接头的平台

--stringency 限定最少与adapter序列重叠的碱基数，默认接头序列：AGA

-q --quality 切除质量得分低于设置值的序列，默认为20

--max_n去除含有N碱基数>n的序列

--phred33/--phred64

-o 设定输出目录，必须是已存在的目录，否则运行将报错

#可以先用单个样本试一下
trim_galore --phred33 -q 20 --length 36    \
--stringency 3 --fastqc --paired --max_n 3 \
-o ./trim_galore/                          \
./fastq_25000/SRR1039510_1.fastq.gz        \
./fastq_25000/SRR1039510_2.fastq.gz
#用\分隔同一行，使代码看起来更美观

ls ./fastq_25000/*_1.fastq.gz | cut -d"/" -f3 | cut -d "_" -f1 > sample.ID
cat sample.ID | while read id
do
  trim_galore --phred33 -q 20 --length 36    \
  --stringency 3 --fastqc --paired --max_n 3 \
  -o ./trim_galore/                          \
  ./fastq_25000/${id}_1.fastq.gz             \
  ./fastq_25000/${id}_2.fastq.gz
done

ls  -lh  ./trim_galore/
#查看是否有新文件生成，文件大小、时间戳是否正常

multiqc  ./trim_galore/*.zip  -o ./trim_galore/
#使用multiqc整合，可以对比质控前后的结果

内容4：Fastp

fastp过滤的优点是基于C语言，运行速度快

-i -I（大写的i）后接需要过滤的fastq文件

-o -O（大写的o）后接过滤完输出的fastq文件名，小o大o更改位置不影响，但是记得输入输出相应更改位置

-h html文件名

-j json报告文件名

-n 限制N碱基的个数

-q设置碱基质量阈值，默认阈值为15

-u 设置允许不合格碱基阈值，超过比例时，整条read将被丢弃，默认为40

-l 设置read的最小长度，默认是15，和trim_galore保持一致，可以设置为36

-w 可用的线程，默认为2

-z 设置文件压缩比

示例脚本

#单一样本尝试
fastp -i ./fastq_25000/SRR1039510_1.fastq.gz \
      -I ./fastq_25000/SRR1039510_2.fastq.gz \
      -o ./fastp/SRR1039510_1.fastp.fq.gz    \
      -O ./fastp/SRR1039510_2.fastp.fq.gz    \
      -h ./fastp/SRR1039510.html             \
      -j ./fastp/SRR1039510.json             \
      -l 36 -q 20


cat sample.ID | while read id
do
fastp -i ./fastq_25000/${id}_1.fastq.gz \
      -I ./fastq_25000/${id}_2.fastq.gz \
      -o ./fastp/${id}_1.fastp.fq.gz    \
      -O ./fastp/${id}_2.fastp.fq.gz    \
      -h ./fastp/${id}.html             \
      -j ./fastp/${id}.json             \
      -l 36 -q 20
done

生信技能树，生信马拉松，火龙果老师