NGS_lab2

数据预处理

安装 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
  1. 质量评估：FastQC 、MultiQC
  2. 质量控制：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

📄 Con_sequence_1_fastqc.html

📄 Con_sequence_2_fastqc.html

  • 从以上两个示例文件中，我们可以发现需要切的 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.fa

2. 质量控制

使用 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

📄 fastp_Con_sequence_1_fastqc.html

📄 fastp_Con_sequence_2_fastqc.html

  • 可以发现 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.gz
  • FastQC 再来质检一下
# 需切换自己虚拟环境下的 fastqc，可点击网页右侧 "1.质量评估" 查看，不再赘述
# FastQC 质检 (自动根据 *.fastq.gz 个数来调整线程数)
fastqc --threads $(ls *.fastq.gz | wc -w) *.fastq.gz

📄 cutadapt_Con_sequence_1_fastqc.html

📄 cutadapt_Con_sequence_2_fastqc.html

  • 可以发现 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

📄 trimmomatic_Con_sequence_1_paired_fastqc.html

📄 trimmomatic_Con_sequence_2_paired_fastqc.html

📄 trimmomatic_Con_sequence_1_unpaired_fastqc.html

📄 trimmomatic_Con_sequence_2_unpaired_fastqc.html

  • 可以发现 adapter 们又双叒叕全都是一条直线啦！

3. 汇总报告

  • 使用 MultiQC 汇总 FastQC 报告
# 先回到 Lab2
cd ~/myNGS/Lab2

# 激活 multiqc 虚拟环境
# 假如你使用 anaconda
conda activate multiqc
# 假如你使用 micromamba
micromamba activate multiqc

# 一键汇总剪切过 adapter 的报告
multiqc ./fastp/ ./cutadapt/ ./trimmomatic/

📄 multiqc_report.html

  • 可以观察到，三种软件的 adapter 剪切结果大差不差

Contents