生信技能树-day17 转录组上游分析-数据质控、过滤

生信菜鸟团

发布于 2024-06-25 20:20:12

960

发布于 2024-06-25 20:20:12

文章被收录于专栏：生信菜鸟团生信菜鸟团

从我们生信技能树历年的几千个马拉松授课学员里面募集了一些优秀的创作者，某种意义来说是传承了我们生信技能树的知识整理和分享的思想！

今天的是三周合计15天的数据挖掘授课学员一点一滴整理的授课知识点笔记哦，还有互动练习题哈，欢迎大家点击文末的阅读原文去关注我们学员的公众号哦！

数据质控

fastqc

目标：使用fastqc对原始数据进行质量评估

如果想查看历史命令记录，可以使用history | less或者history | tail

fastqc常用参数

实际代码：

# 激活conda环境
conda activate rna

# 连接数据到自己的文件夹
cd $HOME/project/Human-16-Asthma-Trans/data/rawdata
ln -s /home/t_rna/data/airway/fastq_raw25000/*gz  ./

# 使用FastQC软件对单个fastq文件进行质量评估，结果输出到qc/文件夹下
# 方式一：直接运行
fastqc -t 6 -o ./ SRR*.fastq.gz
# -t 6指同时处理的文件数6个，-o ./是指输出的目录为当前文件夹，后面没有参数跟着输入文件名称，用>重定向到qc.log日志文件

# 方式二：后台运行
nohup fastqc -t 6 -o ./ SRR*.fastq.gz >qc.log &
# 命令前后加入nohup &指将命令在后台运行

# 方式三：写入sh脚本，用nohup &在后台运行脚本
vim qc.sh
i
fastqc -t 6 -o ./ SRR*.fastq.gz
# 按esc
:wq
nohup bash qc.sh >qc.log &
# 如何知道运行完了？
cat qc.log

# 质控的结果给每个fq文件生成一个zip文件和一个html网页报告
# zip文件里是一些报告图片等，主要需要看的是html网页报告，可以将其下载到本地通过浏览器打开

FastQC Report解释

basic statistics

Filename：原始分析文件名称
File type：文件类型，base cells表示碱基识别，带有fq序列信息文件
Encoding：编码ASCII编码的质量值系统
Total Sequences：测序的总read数
Filtered Sequences：过滤的read数
Sequence Length：提供最长和最短read的length，如果所有read一样长，只有一个数值
%GC：GC含量

生物学中统计测序数据量的大小的单位（# of Bases），和计算机中信息计量单位（Size）不一样：

b：表示base pairs，即碱基对

M：Million，百万

1Mb：即10^6个base pairs

1Gb：即10^9个base pairs

per base sequence quality

fq文件中每一个位置上的所有read的碱基质量值（Q值）的箱图

比如第一个箱子指文件中所有的read第一个位置的碱基（base），一个文件中有25000个reads，就有25000个值产生这一个箱式图

Q值大于30指准确率＞99.9%，箱式图落在绿色的区域指测序准确率较高，随着测序读长变长，碱基质量值会下降

per tile sequence quality

横坐标表示碱基位置，纵坐标表示flowcell中的tile

这张图用于判断测序质量是否是由于flowcell中的特定位置出现了异常，好的数据应该是全蓝的，如果出现红色，说明flowcell中的某一个tile出现了异常，导致数据质量差

per sequence quality scores

横坐标为Q值，纵坐标为打分为该Q值的reads数

如果出现较多的read数对应低Q值，则说明数据质量较差

per base sequence content

横坐标为各碱基位置，纵坐标为碱基百分比

表示每个碱基位置上ATGC含量的分布，用于检查有无ATGC分离

理论上，G和C、A和T的含量每个测序循环上应分别相等，且整个测序过程稳定不变，呈水平线

前面的位置有一些抖动是正常的，因为是引物的位置

per sequence GC content

文件中每个序列的整个长度的GC内容

横坐标为平均GC含量，纵坐标为GC含量对应的序列数量

蓝色线为理论值，红色线为测量值，越接近越好，应为单峰分布

per base N content

N指的是荧光信号识别的时候识别不出ATCG，指未知的没有测出来的值，含量越低越好

横坐标表示碱基位置，纵坐标表示N的含量

sequence length distribution

read长度分布图，横坐标为reads长度，纵坐标为数量

sequence duplication levels

序列的重复率，大部分序列只出现一次，少部分序列出现的数量较多，横坐标是duplicate的次数，纵坐标是重复的百分比

只统计前100000reads，重复率大于10的reads会被统计为>10，重复率过高说明测序失败

overrepresented sequences

指过表达的reads，和上面的报告意义一样，指出重复率较高的reads

adapter content

接头含量，表示序列中两端接头的情况

http://www.bioinformatics.babraham.ac.uk/projects/fastqc/中可以看质量好的测序结果和质量差的测序结果的示例报告

实际操作过程中需要每个fq文件都单独查看一遍质控报告吗？当然不需要

使用MultiQC整合FastQC结果：

multiqc *zip

同样生成一个html文件，下载到本地打开后可以看到所有fq文件整合的质控报告

另外，以上命令可以使用绝对路径运行，好处是可以在任何目录位置运行，并且可以不激活小环境

# 使用绝对路径运行
multiqc=/home/t_rna/miniconda3/envs/rna/bin/multiqc
fastqc=/home/t_rna/miniconda3/envs/rna/bin/fastqc
fq_dir=$HOME/project/Human-16-Asthma-Trans/data/rawdata
outdir=$HOME/project/Human-16-Asthma-Trans/data/rawdata

# 使用绝对路径运行
$fastqc -t 6 -o $outdir ${fq_dir}/SRR*.fastq.gz >${fq_dir}/qc.log

# 报告整合
$multiqc $outdir/*.zip -o $outdir/ >${fq_dir}/multiqc.log

数据过滤

如何判断你的数据需要过滤？

原始序列质量控制的标准为：

去除含接头的reads
过滤去除低质量值的数据，确保数据质量
去除含有N的比例大于5%的reads

trim_galore过滤

官网：http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/

常用参数：

# 激活小环境
conda activate rna
# 新建文件夹trim_galore
cd $HOME/project/Human-16-Asthma-Trans/data/cleandata/trim_galore

# 单个样本的运行命令
trim_galore -q 20 --length 20 --max_n 3 --stringency 3 --fastqc --paired -o ./ ../../rawdata/SRR1039510_1.fastq.gz ../../rawdata/SRR1039510_2.fastq.gz
# 其中-q 20指过滤q值小于20的碱基，--length 20指剪切后小于20长度的reads丢弃掉，--max_n 3指n的数量大于3的reads丢弃，--stringency 3指接头序列比对大于3个碱基的丢弃，--fastqc指对过滤后的数据进行fastqc，--paired指成对的reads如果一条质量不合格，则把另一条也丢弃,-o指设定输出目录，默认为当前目录
# 后面为输出目录和输入文件，输出目录为当前目录，输入文件以相对路径表示，为rawdata文件夹中两个成对的fq文件

#对于多个样本，变的只有样本号，因此可以：
# 先生成一个变量,为样本ID
ls $HOME/project/Human-16-Asthma-Trans/data/rawdata/*_1.fastq.gz | awk -F'/' '{print $NF}' | cut -d'_' -f1 >ID
# 把rawdata文件夹中的fastq.gz文件传到awk中进行切割，-F '/'指以/为分割进行切割，'{print$NF}'指print最后一列，然后传到cut中去掉_，得到样本ID的列表
cat ID
# 看看对不对

# 多个样本 vim trim_galore.sh，以下为sh的内容
rawdata=$HOME/project/Human-16-Asthma-Trans/data/rawdata
cleandata=$HOME/project/Human-16-Asthma-Trans/data/cleandata/trim_galore
ID=$HOME/project/Human-16-Asthma-Trans/data/cleandata/trim_galore/ID
cat $ID | while read id
do
  trim_galore -q 20 --length 20 --max_n 3 --stringency 3 --fastqc --paired -o ${cleandata}    ${rawdata}/${id}_1.fastq.gz     ${rawdata}/${id}_2.fastq.gz
done

#可以创建脚本，把以上命令输入trim_galore.sh并在后台运行
vim trim_galore.sh
i
# 复制黏贴以上命令写入脚本
# 按ESC
:wq

# 提交任务到后台 可以 用bash或者sh都行 
nohup bash trim_galore.sh >trim_galore.log &

# 使用MultiQc整合FastQC结果
multiqc *.zip

## ==============================================
## 补充技巧：使用掐头去尾获得样本ID
ls $rawdata/*_1.fastq.gz | while read id
do
name=${id##*/}
name=${name%_*}
echo "trim_galore -q 20 --length 20 --max_n 3 --stringency 3 --fastqc --paired -o ${cleandata} ${rawdata}/${name}_1.fastq.gz ${rawdata}/${name}_2.fastq.gz "
done

任务管理相关命令总结

任务投递

前台运行：直接运行

后台运行：nohup &

前台转后台：运行后ctrl+z可以暂停，bg %1可以将编号1的任务转为后台

后台转前台：fg

暂停/终止

终止任务：ctrl+c，kill -9 %1

暂停任务：ctrl+z

任务查看

top htop

ps fxww

jobs

检查脚本内容

echo "trim_galore -q 20 --length 20 --max_n 3 --stringency 3 --fastqc --paired -o ${cleandata}

用echo可以将命令中的变量打印出来看是否有误

fastp过滤

常用参数

cd $HOME/project/Human-16-Asthma-Trans/data/cleandata/fastp

# 单样本：
#反斜线后回车，可以换行输入长命令，这样比较整齐，最后一行不用输入反斜线
fastp -l 20 -q 20 --compression=6 \
  -i ../../rawdata/SRR1039510_1.fastq.gz \
  -I ../../rawdata/SRR1039510_2.fastq.gz  \
  -o ./SRR1039510_clean_1.fq.gz  \
  -O ./SRR1039510_clean_2.fq.gz \
  -R SRR1039510 \
  -h SRR1039510.fastp.html \
  -j ./SRR1039510.fastp.json

# 定义样本ID
ls $HOME/project/Human-16-Asthma-Trans/data/rawdata/*_1.fastq.gz | awk -F'/' '{print $NF}' | cut -d'_' -f1 >ID

## 写fastp命令并echo一下看看有没有错误
cat ../trim_galore/ID | while read id
do echo "
fastp -l 20 -q 20 --compression=6 \
  -i ${rawdata}/${id}_1.fastq.gz \
  -I ${rawdata}/${id}_2.fastq.gz \
  -o ${cleandata}/${id}_clean_1.fq.gz \
  -O ${cleandata}/${id}_clean_2.fq.gz \
  -R ${cleandata}/${id} \
  -h ${cleandata}/${id}.fastp.html \
  -j ${cleandata}/${id}.fastp.json 
  "
done

## 实际运行时把以下命令写入fastp.sh脚本中并运行
rawdata=$HOME/project/Human-16-Asthma-Trans/data/rawdata
cleandata=$HOME/project/Human-16-Asthma-Trans/data/cleandata/fastp
id=$HOME/project/Human-16-Asthma-Trans/data/cleandata/fastp/ID
cat ../trim_galore/ID | while read id
do
fastp -l 20 -q 20 --compression=6 \
  -i ${rawdata}/${id}_1.fastq.gz \
  -I ${rawdata}/${id}_2.fastq.gz \
  -o ${cleandata}/${id}_clean_1.fq.gz \
  -O ${cleandata}/${id}_clean_2.fq.gz \
  -R ${cleandata}/${id} \
  -h ${cleandata}/${id}.fastp.html \
  -j ${cleandata}/${id}.fastp.json 
done

# 运行fastp脚本
nohup bash fastp.sh >fastp.log &

# fastp同样会生成html文件，下载到本地进行查看

# 整合成一个文件
multiqc *json

数据过滤前后比较

# 进入过滤目录
cd $HOME/project/Human-16-Asthma-Trans/data/cleandata/trim_galore

# 原始数据
zcat $rawdata/SRR1039510_1.fastq.gz | paste - - - - > raw.txt

#  过滤后的数据
zcat SRR1039510_1_val_1.fq.gz |paste - - - - > trim.txt
awk '(length($4)<63){print$1}' trim.txt > Seq.ID

head -n 100 Seq.ID > ID100
grep -w -f ID100 trim.txt | awk '{print$1,$4}' > trim.sm
grep -w -f ID100 raw.txt | awk '{print$1,$4}' > raw.sm
paste raw.sm trim.sm | awk '{print$2,$4}' | tr ' ' '\n' |less -S