NGS Lab2 (数据预处理)
2024.10.28 追加: PolyA 好像没必要切的那么干净......具体什么情况需要切干净有待商榷
安装 Trimmomatic
( 挑个你喜欢的方法安装就好～ )
使用 anaconda 安装
- 新建虚拟环境
trimmomatic
conda create -n trimmomatic
# 等待一会儿后输入 y , 然后回车- 进入虚拟环境
trimmomatic
conda activate trimmomatic- 安装
trimmomatic
# 注意！如果不添加 -c conda-forge, 则默认无法下载最新版本的 trimmomatic
conda install -c bioconda -c conda-forge trimmomatic
# 等待一会儿后输入 y , 然后回车
# 检验是否安装成功
# 注意! 这里的 version 前面只有一个 -，不是 --
trimmomatic -version
# 我的输出是:
# 0.39使用 micromamba 安装
- 新建虚拟环境
trimmomatic
# 和 anaconda 不同, micromamba 在新建虚拟环境的时候需要加上 env
micromamba env create -n trimmomatic- 进入虚拟环境
trimmomatic
micromamba activate trimmomatic- 安装
trimmomatic
# 注意！如果不添加 -c conda-forge, 则默认无法下载最新版本的 sra-tools
micromamba install -c bioconda -c conda-forge trimmomatic
# 等待一会儿后输入 Y , 然后回车
# 检验是否安装成功
# 注意! 这里的 version 前面只有一个 -，不是 --
trimmomatic -version
# 我的输出是:
# 0.39安装 cutadapt
( 挑个你喜欢的方法安装就好～ )
使用 anaconda 安装
- 新建虚拟环境
cutadapt
conda create -n cutadapt
# 等待一会儿后输入 y , 然后回车- 进入虚拟环境
cutadapt
conda activate cutadapt- 安装
cutadapt
# 注意！如果不添加 -c conda-forge, 则默认无法下载最新版本的 cutadapt
conda install -c bioconda -c conda-forge cutadapt
# 等待一会儿后输入 y , 然后回车
# 检验是否安装成功
cutadapt --version
# 我的输出是:
# 4.9使用 micromamba 安装
- 新建虚拟环境
cutadapt
# 和 anaconda 不同, micromamba 在新建虚拟环境的时候需要加上 env
micromamba env create -n cutadapt- 进入虚拟环境
cutadapt
micromamba activate cutadapt- 安装
cutadapt
# 注意！如果不添加 -c conda-forge, 则默认无法下载最新版本的 sra-tools
micromamba install -c bioconda -c conda-forge cutadapt
# 等待一会儿后输入 Y , 然后回车
# 检验是否安装成功
cutadapt --version
# 我的输出是:
# 4.9安装 fastp
( 挑个你喜欢的方法安装就好～ )
使用 anaconda 安装
- 新建虚拟环境
fastp
conda create -n fastp
# 等待一会儿后输入 y , 然后回车- 进入虚拟环境
fastp
conda activate fastp- 安装
fastp
# 注意！如果不添加 -c conda-forge, 则默认无法下载最新版本的 fastp
conda install -c bioconda -c conda-forge fastp
# 等待一会儿后输入 y , 然后回车
# 检验是否安装成功
fastp --version
# 我的输出是:
# fastp 0.23.4使用 micromamba 安装
- 新建虚拟环境
fastp
# 和 anaconda 不同, micromamba 在新建虚拟环境的时候需要加上 env
micromamba env create -n fastp- 进入虚拟环境
fastp
micromamba activate fastp- 安装
fastp
# 注意！如果不添加 -c conda-forge, 则默认无法下载最新版本的 sra-tools
micromamba install -c bioconda -c conda-forge fastp
# 等待一会儿后输入 Y , 然后回车
# 检验是否安装成功
fastp --version
# 我的输出是:
# fastp 0.23.4任务
- 数据:
/public/workspace/shaojf/Course/NGS/DataSets/Lab2/Con_sequence_{1,2}.fastq.gz
- 质量评估：
FastQC、MultiQC - 质量控制：
fastqc后可发现原始数据adapter含量较高, 请用trimmomatic、cutadapt或fastp中的任意一个工具（或所有）完成切除 adapter 任务，并在完成后再次运行fastqc检查切除效率。（提示，可以用MultiQC将处理前后、不同工具预处理后的结果汇总在同一个报告中)
参考代码
0. 准备工作
# 新建文件夹
mkdir -p ~/myNGS/Lab2
cd ~/myNGS/Lab2
# 原始数据软链接
ln -s /public/workspace/shaojf/Course/NGS/DataSets/Lab2/Con_sequence_*.fastq.gz .1. 质量评估
- 寻找 adapter 类型
# 激活 fastqc 虚拟环境
# 假如你使用 anaconda
conda activate fastqc
# 假如你使用 micromamba
micromamba activate fastqc
# 新建文件夹 raw_fastqc, 用于保存原始数据的质检结果
mkdir raw_fastqc
# FastQC 质检 (自动根据 *.fastq.gz 个数来调整线程数)
fastqc --threads $(ls *.fastq.gz | wc -w) *.fastq.gz -o ./raw_fastqc从以上两个示例文件中，我们可以发现需要切的 adapter 有三种，为
Nextera Transposase Sequence、PolyG、PolyA查看 adapter 具体序列
# adapter 序列来源: fastqc 的配置文件
# (具体路径根据你 fastqc 所在的虚拟环境而定, 这里使用 Lab1 中所安装的路径)
# 假如你使用的是 anaconda
cat $(find ~/anaconda3/envs/fastqc/ -name "adapter_list.txt")
# 假如你使用的是 micromamba
cat $(find ~/micromamba/envs/fastqc/ -name "adapter_list.txt")
# 最后几行即为 fastqc 可检测到的 adapter：
# Illumina Universal Adapter AGATCGGAAGAG
# Illumina Small RNA 3' TGGAATTCTCGG
# Illumina Small RNA 5' Adapter GATCGTCGGACT
# Nextera Transposase Sequence CTGTCTCTTATA
# PolyA AAAAAAAAAAAA
# PolyG GGGGGGGGGGGG- 我们可以建立 adapter 序列文件
adapter_Nextera_Transposase_Sequence_Poly_G_A.fa，方便后续的 adapter 剪切工作
echo -e ">Nextera_Transposase_Sequence\nCTGTCTCTTATA\n>Poly_G\nGGGGGGGGGGGG\n>Poly_A\nAAAAAAAAAAAA" > adapter_Nextera_Transposase_Sequence_Poly_G_A.fa2. 质量控制
使用 fastp 切除 adapter
- 先尝尝甜头,
fastp, 启动!
# 激活 fastp 虚拟环境
# 假如你使用 anaconda
conda activate fastp
# 假如你使用 micromamba
micromamba activate fastp- 新建文件夹
fastp用于保存fastp的输出结果
mkdir fastp
cd fastp- 要开始切咯～
# fastp 效率最高, 16线程下仅需20秒左右
fastp -i ../Con_sequence_1.fastq.gz -I ../Con_sequence_2.fastq.gz \
-o fastp_Con_sequence_1.fastq.gz -O fastp_Con_sequence_2.fastq.gz \
--detect_adapter_for_pe \
-5 -c \
--cut_mean_quality 20 \
--length_required 36 \
--adapter_fasta ../adapter_Nextera_Transposase_Sequence_Poly_G_A.fa \
-j fastp_Con_sequence.json \
-h fastp_Con_sequence.html \
-w 16
# 参数解析:
# -5 (--cut_front) 需要指定滑动窗口方向 (5'->3'), 否则会跳过切除低质量碱基
# -c (--correctionenable) 加上 -c 才会进行碱基校正 (仅适用于双末端数据, 默认不启用)
# -w 线程数 (别贪心哦！目前 fastp 最多就支持16线程)FastQC再来质检一下
# 需切换自己虚拟环境下的 fastqc，可点击网页右侧 "1.质量评估" 查看，不再赘述
# FastQC 质检 (自动根据 *.fastq.gz 个数来调整线程数)
fastqc --threads $(ls *.fastq.gz | wc -w) *.fastq.gz- 可以发现 adapter 们全都是一条直线啦！
使用 cutadapt 切除 adapter
- 激活
cutadapt虚拟环境
# 假如你使用 anaconda
conda activate cutadapt
# 假如你使用 micromamba
micromamba activate cutadapt- 新建文件夹
cutadapt用于保存cutadapt的输出结果
# 先回到 Lab2
cd ~/myNGS/Lab2
mkdir cutadapt
cd cutadapt- 开切！(个人最喜欢
cutadapt，因为它的进度条是一把小剪刀！)
# cutadapt 比 fastp 慢点, 20线程下需70秒左右
cutadapt --quality-base=33 -q 20 -m 36 \
-a file:../adapter_Nextera_Transposase_Sequence_Poly_G_A.fa \
-A file:../adapter_Nextera_Transposase_Sequence_Poly_G_A.fa \
-o cutadapt_Con_sequence_1.fastq.gz \
-p cutadapt_Con_sequence_2.fastq.gz \
--cores 20 \
../Con_sequence_1.fastq.gz ../Con_sequence_2.fastq.gzFastQC再来质检一下
# 需切换自己虚拟环境下的 fastqc，可点击网页右侧 "1.质量评估" 查看，不再赘述
# FastQC 质检 (自动根据 *.fastq.gz 个数来调整线程数)
fastqc --threads $(ls *.fastq.gz | wc -w) *.fastq.gz- 可以发现 adapter 们又全都是一条直线啦！
使用 trimmomatic 切除 adapter
(最苦不过 trimmomatic)
- 激活
trimmomatic虚拟环境
# 假如你使用 anaconda
conda activate trimmomatic
# 假如你使用 micromamba
micromamba activate trimmomatic- 新建文件夹
trimmomatic用于保存trimmomatic的输出结果
# 先回到 Lab2
cd ~/myNGS/Lab2
mkdir trimmomatic
cd trimmomatic- 切！切！切！
# trimmomatic 慢的嘞... 20线程下耗时将近两分半 (一坤钟 (bushi))
# 而且还不稳定, 同样的命令有的时候跑一次还不一定成功... 可多尝试几次
trimmomatic PE -phred33 ../Con_sequence_1.fastq.gz ../Con_sequence_2.fastq.gz \
trimmomatic_Con_sequence_1_paired.fastq.gz \
trimmomatic_Con_sequence_1_unpaired.fastq.gz \
trimmomatic_Con_sequence_2_paired.fastq.gz \
trimmomatic_Con_sequence_2_unpaired.fastq.gz \
ILLUMINACLIP:../adapter_Nextera_Transposase_Sequence_Poly_G_A.fa:0:15:5 \
LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 \
-threads 20
# 参数解析:
# 未完待续...- 麻烦
FastQC再来质检一下
# 需切换自己虚拟环境下的 fastqc，可点击网页右侧 "1.质量评估" 查看，不再赘述
# FastQC 质检 (自动根据 *.fastq.gz 个数来调整线程数)
fastqc --threads $(ls *.fastq.gz | wc -w) *.fastq.gz- 可以发现 adapter 们又双叒叕全都是一条直线啦！
3. 汇总报告
- 使用
MultiQC汇总FastQC报告
# 先回到 Lab2
cd ~/myNGS/Lab2
# 激活 multiqc 虚拟环境
# 假如你使用 anaconda
conda activate multiqc
# 假如你使用 micromamba
micromamba activate multiqc
# 一键汇总剪切过 adapter 的报告
multiqc ./fastp/ ./cutadapt/ ./trimmomatic/ ./raw_fastqc/- 可以观察到，三种软件的 adapter 剪切结果大差不差