宏基因组流程-质控（Trimmomatic）

这几期我们将专注与宏基因组的分析流程，在这个过程中我们将会详细拆分和解析每一个步骤。此外，还会对每一个步骤所使用的主流生信软件做一个详细的教程。（这一期涉及到一些原理，文字比较多，谨慎食用）

---

宏基因组学[1]（Metagenomics）是一种研究微生物群落中所有基因的分析方法，而不是针对单一生物体的基因组研究。它涵盖了环境中的多种微生物，包括细菌、古菌、真菌和病毒等。宏基因组学通过直接从环境样品中提取DNA并进行测序，以了解微生物群落的遗传多样性、功能和生态学特性。宏基因组学可以用于研究土壤、水体、肠道、海洋、湖泊等多种生态系统中的微生物群落。随着测序成本的降低，宏基因组测序在多个领域收到青睐，使其逐渐成为炙手可热的组学分析手段之一。其流程如下：

---

质 控

宏基因组流程的第一个步骤就是质控，首先我们要搞清楚什么叫做质控。FASTQ序列格式为如下，以@开头，同时使用ASCII编码来代表每一个碱基的质量评分（“+”后的符号，与碱基一一对应）。质控其实就是对FASTQ序列尾端的低质量碱基做一个切除处理，提高数据的准确性。此外，质控对序列数据的处理还包括：①去除序列中的接头Adapter；②去除重复序列；③去除由于测序错误等原因导致序列产生的噪音、杂质等；④去除残留的P5或P7序列。

@HISEQ:669:HW5Y7BCXY:1:1101:10017:75689
ACTCGCACAGGACCGTGCTCCCCCGCCAATTCCTTT
+
DDDDCIIIHIIIIHIIIIIIIIHIIHHIIIIIIIII

质控软件-Trimmomatic

Trimmomatic[2]是一款基于illumina平台的质控软件，它既能去除FASTQ序列的Adapter，也能剪切FASTQ序列中的较低质量的碱基。Trimmomatic支持单端SE、双端PE数据的处理。

安装

#conda安装
conda install -c bioconda trimmomatic
#查看帮助文档
trimmomatic --help

#下载
wget https://github.com/usadellab/Trimmomatic/files/5854859/Trimmomatic-0.39.zip
#解压
unzip Trimmomatic-0.39.zip
cd Trimmomatic
tree
#解压后的文件如下，分别是存放接头文件的文件夹adapters，包含单端SE和双端PE接头文件。
#conda安装方式的接头文件一般目录为conda下的share文件夹的trimmomatic文件夹内，比如我的存放位置为：/miniconda3/share/trimmomatic-0.39-2/adapters；
#Java编译的运行程序trimmomatic-0.39.jar

Trimmomatic数据过滤原理

Trimmomatic对序列的质控包含两个模式：Simple mode和Palindrome mode。 图中四个颜色方块分别为：
浅蓝色：人工序列，可以理解为识别Adapter序列的引物（下文称引物序列）；
蓝色：adapter接头序列；
绿色：有效序列，即测序的生物序列；
红色：切除掉的不合格序列。\

• Sample mode模式：该模式是对单端序列测序的质控过程，如下图：
  A模式下，引物序列在read的起始位置（即5'端）识别到接头序列，根据测序原理该read往后部分均不含有可用序列，应舍去。很多的文章也只是提到这么一句，但并未做过多解释，很容易让人一头雾水。这里我详细讲一下，首先我们知道DNA的合成 方向是5'-3'的方向进行的，那么对应的模板链顺序应该是从3'-5'的顺序。在测序中，接头adapter的作用是包括将文库的序列固定在Flowcell（流动池，是样品进行测序过程发生的位置，为一块玻璃板，通常有2或8个lane（泳道）），从模板链3’端开始一边合成碱基一边测序，测序出来的序列即read，如果从一开始就是识别到接头序列，那么表示模板链均为Adapter序列，不可用。
  B模式与A模式相同；
  C模式与D模式一样，read末端包含部分接头序列，这种情况的出现是因为“测穿”，什么意思呢，就是说文库的片段短于测序的长度（如测序长度250bp，而目的片段只有120bp），最后测序的时候就会超过目的片段的长度，测到接头序列。trimmomatic会识别并切除，但如果接头序列过短，比接头匹配参数设置的最短长度还短，那么就无法去除。\

• Palindrome mode：该模式是使用测序的Forward read和Revered read比对来去除污染的接头序列。
  A模式：正向测序和反向测序有部分完全反向互补，这表明两条序列从起始位置均包含了接头序列。这样的 reads 不包含任何有用序列，正反向测序 reads 都被丢弃。
  B模式：如果两条序列中间部分比对不上，前后对上了，那么中间的序列就是有用的序列，保留这一部分；
  C模式：与B模式相同，只不过是更长的有用序列；
  D模式：仅双端的序列比对上，接头位置并对有overlap，只需去除接头。

Simple mode

Palindrome mode

使用

#运行
java -jar Trimmomatic-0.39/trimmomatic-0.39.jar --help
#帮助文档如下
Usage:
       PE [-version] [-threads <threads>] [-phred33|-phred64] [-trimlog <trimLogFile>] [-summary <statsSummaryFile>] [-quiet] [-validatePairs] [-basein <inputBase> | <inputFile1> <inputFile2>] [-baseout <outputBase> | <outputFile1P> <outputFile1U> <outputFile2P> <outputFile2U>] <trimmer1>...
   or:
       SE [-version] [-threads <threads>] [-phred33|-phred64] [-trimlog <trimLogFile>] [-summary <statsSummaryFile>] [-quiet] <inputFile> <outputFile> <trimmer1>...
   or:
       -version
#示例
#双端模式下
java -jar Trimmomatic-0.39/trimmomatic-0.39.jar PE \
-phred33  Sample_1.fq.gz Sample_2.fq.gz \
Sample_1_paired.fq.gz Sample_1_unpaired.fq.gz Sample_2_paired.fq.gz Sample_2_unpaired.fq.gz \
ILLUMINACLIP:Trimmomatic-0.39/adapters/TruSeq3-PE.fa:2:30:10  \
LEADING:3  TRAILING:3  SLIDINGWINDOW:5:20  MINLEN:50 -threads 24

#单端模式下
java -jar Trimmomatic-0.39/trimmomatic-0.39.jar SE \
-phred33 sample.fq.gz \
sample_clean.fq.gz \
ILLUMINACLIP:Trimmomatic-0.39/adapters/TruSeq3-SE.fa:2:30:10  \
LEADING:3  TRAILING:3  SLIDINGWINDOW:5:20  MINLEN:50 -threads 24

#多个样品进行质控时，调用for循环
for filename in *_1.fq.gz
do i=$(basename $filename _1.fq.gz)
echo $i
java -jar Trimmomatic-0.39/trimmomatic-0.39.jar PE -phred33 \
${i}_1.fastq.gz ${i}_2.fastq.gz \
${i}_1_paired.fastq  ${i}_1_unpaired.fastq ${i}_2_paired.fastq  ${i}_2_unpaired.fastq \
ILLUMINACLIP:Trimmomatic-0.39/adapters/TruSeq3-PE.fa:2:30:10 \
LEADING:3 TRAILING:3 SLIDINGWINDOW:5:20 MINLEN:50 -threads 24 \
done

Trimmomatic参数

PE模式下，Trimmomatic运行生成的文件"paired"和"unpaired"的文件，其中paired使我们需要的文件，而unpaired是正向序列和反向序列没有匹配上的序列，即两条序列不互补，应舍去。
Trimmomatic具体的参数如下：\

• ILLUMINACLIP:该参数包含多个参数，
  ILLUMINACLIP:<fastaWithAdaptersEtc>:<seed mismatches>:<palindrome clip threshold>:<simple clip threshold>:<minAdapterLength>:<keepBothReads>
  如:ILLUMINACLIP:Trimmomatic-0.39/adapters/TruSeq3-PE.fa:2:30:10
  fastaWithAdaptersEtc:Trimmomatic-0.39/adapters/TruSeq3-PE.fa指定包含接头和引物序列的 fasta 文件路径,可以使用其自带的接头文件，见上文；
  seed mismatches:2为seed 搜索时允许的错配碱基个数；
  palindrome clip threshold:指定针对 PE 的 palindrome clip 模式下，正向序列和反向序列至少多少比对分值（上图中Pailedrome C 模式），才会进行接头切除，这里是30。这个值计算方法是每匹配一个碱基，评分增加0.6，而每错配一次，比对评分减少Q/10
  simple clip threshold:指定切除接头序列的最低比对分值，通常 7-15 之间。
  minAdapterLength: PE模式有效，指定最小比对到的adapter长度，如果未指定，则由于历史原因，默认为8个碱基，但可以设置到1，切除1bp的Adapter序列。
  keepBothReads:PE模式有效，Pailedrome D 模式，删除适配器序列，Forward read与Reverse read反向互补。Trimmomatic默认删除Reverse read。指定该参数true，Reverse read也将被保留，这可能会很有用，例如，如果下游工具不能处理成对和非配对读取的组合。一般设置为true。

• SLIDINGWINDOW: 从read的5'端开始滑动窗口，当窗口内的平均质量小于阈值时，切除后边的碱基。上述示例代码中5：20表示5个bp长度的滑窗，过滤质量阈值低于20的碱基。

• LEADING: 切除read的开头质量低于阈值的碱基，这里为3。

• TRAILING: 切除read的末尾质量低于阈值的碱基，这里为3。

• CROP: 按照该指定参数长度，将所有序列切除为相同长度。

• HEADCROP: 从read的5'端开始切除指定长度的碱基数。

• MINLEN:最小保留长度，低于该长度的序列将被舍弃 。

• AVGQUAL:舍弃低于该质量阈值的序列。

• TOPHRED33: 将质量分数转换为Phred-33

• TOPHRED64: 将质量分数转换为Phred-64

参考文献

[1] Handelsman J, Rondon MR, Brady SF, Clardy J, Goodman RM. Molecular biological access to the chemistry of unknown soil microbes: a new frontier for natural products. Chem Biol. 1998 Oct;5(10):R245-9. doi: 10.1016/s1074-5521(98)90108-9. PMID: 9818143.

[2] Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014 Aug 1;30(15):2114-20. doi: 10.1093/bioinformatics/btu170. Epub 2014 Apr 1. PMID: 24695404; PMCID: PMC4103590.

结语

这一期的内容比较多，涉及了实操和原理部分。部分只是想使用软件的同行来说，可以只看实际操作的代码部分。
最后再推销一下我们的公众号，大家动动小手，点击下方关注我们，您的关注是我们最大的更新动力。

宏基因组流程-质控（Trimmomatic）

这几期我们将专注与宏基因组的分析流程，在这个过程中我们将会详细拆分和解析每一个步骤。此外，还会对每一个步骤所使用的主流生信软件做一个详细的教程。（这一期涉及到一些原理，文字比较多，谨慎食用）

---

宏基因组学[1]（Metagenomics）是一种研究微生物群落中所有基因的分析方法，而不是针对单一生物体的基因组研究。它涵盖了环境中的多种微生物，包括细菌、古菌、真菌和病毒等。宏基因组学通过直接从环境样品中提取DNA并进行测序，以了解微生物群落的遗传多样性、功能和生态学特性。宏基因组学可以用于研究土壤、水体、肠道、海洋、湖泊等多种生态系统中的微生物群落。随着测序成本的降低，宏基因组测序在多个领域收到青睐，使其逐渐成为炙手可热的组学分析手段之一。其流程如下：

---

质 控

宏基因组流程的第一个步骤就是质控，首先我们要搞清楚什么叫做质控。FASTQ序列格式为如下，以@开头，同时使用ASCII编码来代表每一个碱基的质量评分（“+”后的符号，与碱基一一对应）。质控其实就是对FASTQ序列尾端的低质量碱基做一个切除处理，提高数据的准确性。此外，质控对序列数据的处理还包括：①去除序列中的接头Adapter；②去除重复序列；③去除由于测序错误等原因导致序列产生的噪音、杂质等；④去除残留的P5或P7序列。

@HISEQ:669:HW5Y7BCXY:1:1101:10017:75689
ACTCGCACAGGACCGTGCTCCCCCGCCAATTCCTTT
+
DDDDCIIIHIIIIHIIIIIIIIHIIHHIIIIIIIII

质控软件-Trimmomatic

Trimmomatic[2]是一款基于illumina平台的质控软件，它既能去除FASTQ序列的Adapter，也能剪切FASTQ序列中的较低质量的碱基。Trimmomatic支持单端SE、双端PE数据的处理。

安装

#conda安装
conda install -c bioconda trimmomatic
#查看帮助文档
trimmomatic --help

#下载
wget https://github.com/usadellab/Trimmomatic/files/5854859/Trimmomatic-0.39.zip
#解压
unzip Trimmomatic-0.39.zip
cd Trimmomatic
tree
#解压后的文件如下，分别是存放接头文件的文件夹adapters，包含单端SE和双端PE接头文件。
#conda安装方式的接头文件一般目录为conda下的share文件夹的trimmomatic文件夹内，比如我的存放位置为：/miniconda3/share/trimmomatic-0.39-2/adapters；
#Java编译的运行程序trimmomatic-0.39.jar

Trimmomatic数据过滤原理

Trimmomatic对序列的质控包含两个模式：Simple mode和Palindrome mode。 图中四个颜色方块分别为：
浅蓝色：人工序列，可以理解为识别Adapter序列的引物（下文称引物序列）；
蓝色：adapter接头序列；
绿色：有效序列，即测序的生物序列；
红色：切除掉的不合格序列。\

• Sample mode模式：该模式是对单端序列测序的质控过程，如下图：
  A模式下，引物序列在read的起始位置（即5'端）识别到接头序列，根据测序原理该read往后部分均不含有可用序列，应舍去。很多的文章也只是提到这么一句，但并未做过多解释，很容易让人一头雾水。这里我详细讲一下，首先我们知道DNA的合成 方向是5'-3'的方向进行的，那么对应的模板链顺序应该是从3'-5'的顺序。在测序中，接头adapter的作用是包括将文库的序列固定在Flowcell（流动池，是样品进行测序过程发生的位置，为一块玻璃板，通常有2或8个lane（泳道）），从模板链3’端开始一边合成碱基一边测序，测序出来的序列即read，如果从一开始就是识别到接头序列，那么表示模板链均为Adapter序列，不可用。
  B模式与A模式相同；
  C模式与D模式一样，read末端包含部分接头序列，这种情况的出现是因为“测穿”，什么意思呢，就是说文库的片段短于测序的长度（如测序长度250bp，而目的片段只有120bp），最后测序的时候就会超过目的片段的长度，测到接头序列。trimmomatic会识别并切除，但如果接头序列过短，比接头匹配参数设置的最短长度还短，那么就无法去除。\

• Palindrome mode：该模式是使用测序的Forward read和Revered read比对来去除污染的接头序列。
  A模式：正向测序和反向测序有部分完全反向互补，这表明两条序列从起始位置均包含了接头序列。这样的 reads 不包含任何有用序列，正反向测序 reads 都被丢弃。
  B模式：如果两条序列中间部分比对不上，前后对上了，那么中间的序列就是有用的序列，保留这一部分；
  C模式：与B模式相同，只不过是更长的有用序列；
  D模式：仅双端的序列比对上，接头位置并对有overlap，只需去除接头。

Simple mode

Palindrome mode

使用

#运行
java -jar Trimmomatic-0.39/trimmomatic-0.39.jar --help
#帮助文档如下
Usage:
       PE [-version] [-threads <threads>] [-phred33|-phred64] [-trimlog <trimLogFile>] [-summary <statsSummaryFile>] [-quiet] [-validatePairs] [-basein <inputBase> | <inputFile1> <inputFile2>] [-baseout <outputBase> | <outputFile1P> <outputFile1U> <outputFile2P> <outputFile2U>] <trimmer1>...
   or:
       SE [-version] [-threads <threads>] [-phred33|-phred64] [-trimlog <trimLogFile>] [-summary <statsSummaryFile>] [-quiet] <inputFile> <outputFile> <trimmer1>...
   or:
       -version
#示例
#双端模式下
java -jar Trimmomatic-0.39/trimmomatic-0.39.jar PE \
-phred33  Sample_1.fq.gz Sample_2.fq.gz \
Sample_1_paired.fq.gz Sample_1_unpaired.fq.gz Sample_2_paired.fq.gz Sample_2_unpaired.fq.gz \
ILLUMINACLIP:Trimmomatic-0.39/adapters/TruSeq3-PE.fa:2:30:10  \
LEADING:3  TRAILING:3  SLIDINGWINDOW:5:20  MINLEN:50 -threads 24

#单端模式下
java -jar Trimmomatic-0.39/trimmomatic-0.39.jar SE \
-phred33 sample.fq.gz \
sample_clean.fq.gz \
ILLUMINACLIP:Trimmomatic-0.39/adapters/TruSeq3-SE.fa:2:30:10  \
LEADING:3  TRAILING:3  SLIDINGWINDOW:5:20  MINLEN:50 -threads 24

#多个样品进行质控时，调用for循环
for filename in *_1.fq.gz
do i=$(basename $filename _1.fq.gz)
echo $i
java -jar Trimmomatic-0.39/trimmomatic-0.39.jar PE -phred33 \
${i}_1.fastq.gz ${i}_2.fastq.gz \
${i}_1_paired.fastq  ${i}_1_unpaired.fastq ${i}_2_paired.fastq  ${i}_2_unpaired.fastq \
ILLUMINACLIP:Trimmomatic-0.39/adapters/TruSeq3-PE.fa:2:30:10 \
LEADING:3 TRAILING:3 SLIDINGWINDOW:5:20 MINLEN:50 -threads 24 \
done

Trimmomatic参数

PE模式下，Trimmomatic运行生成的文件"paired"和"unpaired"的文件，其中paired使我们需要的文件，而unpaired是正向序列和反向序列没有匹配上的序列，即两条序列不互补，应舍去。
Trimmomatic具体的参数如下：\

• ILLUMINACLIP:该参数包含多个参数，
  ILLUMINACLIP:<fastaWithAdaptersEtc>:<seed mismatches>:<palindrome clip threshold>:<simple clip threshold>:<minAdapterLength>:<keepBothReads>
  如:ILLUMINACLIP:Trimmomatic-0.39/adapters/TruSeq3-PE.fa:2:30:10
  fastaWithAdaptersEtc:Trimmomatic-0.39/adapters/TruSeq3-PE.fa指定包含接头和引物序列的 fasta 文件路径,可以使用其自带的接头文件，见上文；
  seed mismatches:2为seed 搜索时允许的错配碱基个数；
  palindrome clip threshold:指定针对 PE 的 palindrome clip 模式下，正向序列和反向序列至少多少比对分值（上图中Pailedrome C 模式），才会进行接头切除，这里是30。这个值计算方法是每匹配一个碱基，评分增加0.6，而每错配一次，比对评分减少Q/10
  simple clip threshold:指定切除接头序列的最低比对分值，通常 7-15 之间。
  minAdapterLength: PE模式有效，指定最小比对到的adapter长度，如果未指定，则由于历史原因，默认为8个碱基，但可以设置到1，切除1bp的Adapter序列。
  keepBothReads:PE模式有效，Pailedrome D 模式，删除适配器序列，Forward read与Reverse read反向互补。Trimmomatic默认删除Reverse read。指定该参数true，Reverse read也将被保留，这可能会很有用，例如，如果下游工具不能处理成对和非配对读取的组合。一般设置为true。

• SLIDINGWINDOW: 从read的5'端开始滑动窗口，当窗口内的平均质量小于阈值时，切除后边的碱基。上述示例代码中5：20表示5个bp长度的滑窗，过滤质量阈值低于20的碱基。

• LEADING: 切除read的开头质量低于阈值的碱基，这里为3。

• TRAILING: 切除read的末尾质量低于阈值的碱基，这里为3。

• CROP: 按照该指定参数长度，将所有序列切除为相同长度。

• HEADCROP: 从read的5'端开始切除指定长度的碱基数。

• MINLEN:最小保留长度，低于该长度的序列将被舍弃 。

• AVGQUAL:舍弃低于该质量阈值的序列。

• TOPHRED33: 将质量分数转换为Phred-33

• TOPHRED64: 将质量分数转换为Phred-64

参考文献

[1] Handelsman J, Rondon MR, Brady SF, Clardy J, Goodman RM. Molecular biological access to the chemistry of unknown soil microbes: a new frontier for natural products. Chem Biol. 1998 Oct;5(10):R245-9. doi: 10.1016/s1074-5521(98)90108-9. PMID: 9818143.

[2] Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014 Aug 1;30(15):2114-20. doi: 10.1093/bioinformatics/btu170. Epub 2014 Apr 1. PMID: 24695404; PMCID: PMC4103590.

结语

这一期的内容比较多，涉及了实操和原理部分。部分只是想使用软件的同行来说，可以只看实际操作的代码部分。
最后再推销一下我们的公众号，大家动动小手，点击下方关注我们，您的关注是我们最大的更新动力。