病原微生物扩增子数据分析实战（一）：bcl2fastq软件完成数据拆分

懂得了那么多道理，却依然过不好这一生。

所以理论归理论，最终要落实到分析代码上，咱们从这一篇开始，介绍一套扩增子数据分析流程。

闲话少叙，首先介绍一下项目背景：

测序平台: Illumina MiSeq，双端（Paired-end）测序
服务器：CentOS 7操作系统，128G内存，CPU 48线程
编程语言：Bash script、Python

NGS测序项目大致分为二个部分：

湿实验部分，即取样、提取、建库以及上机测序
干实验部分，即测序数据的生物信息分析

测序仪运行结束，就完成了生物样本的数字化过程，但此时核酸的信息，具体来说是DNA的序列信息（RNA样本的话会先反转录成DNA再测序），是保存在称为BCL（Binary Base Call）的文件中的，而生物信息分析工作即是从这里开始。

首先，要将BCL文件转化成FASTQ文件格式，过程并不复杂，一条命令就够了。

1
2
3
4
5
6


/path-to/bcl2fastq_V219 \

-R /path-to/200821_M06862_0006_000000000-C9TL7 \

--sample-sheet /path-to/SampleSheet.csv \

--output-dir /path-to/Demultiplexed \

--barcode-mismatches 1 \

--use-bases-mask Y75n,I8,I8,Y75n

bcl2fastq_V219，是illumina官方提供的BCL转FASTQ格式的软件，这里用的版本是v2.19。

-R参数，是下机数据所在的目录，内容如下，其中的RunInfo.xml文件记录了测序的一些重要信息，如总共运行了多少个循环等。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27


200821_M06862_0006_000000000-C9TL7
├── AnalysisError.txt
├── AnalysisLog.txt
├── Basecalling_Netcopy_complete_Read1.txt
├── Basecalling_Netcopy_complete_Read2.txt
├── Basecalling_Netcopy_complete_Read3.txt
├── Basecalling_Netcopy_complete_Read4.txt
├── Basecalling_Netcopy_complete.txt
├── CompletedJobInfo.xml
├── Config
├── Data
├── GenerateFASTQRunStatistics.xml
├── ImageAnalysis_Netcopy_complete_Read1.txt
├── ImageAnalysis_Netcopy_complete_Read2.txt
├── ImageAnalysis_Netcopy_complete_Read3.txt
├── ImageAnalysis_Netcopy_complete_Read4.txt
├── ImageAnalysis_Netcopy_complete.txt
├── InterOp
├── Logs
├── QueuedForAnalysis.txt
├── Recipe
├── RTAComplete.txt
├── RunCheckDetail.txt
├── RunInfo.xml
├── runParameters.xml
├── SampleSheet.csv
└── Thumbnail_Images

--sample-sheet参数，样本清单，内容如下：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26


[Header],,,,,,,,,
IEMFileVersion,5,,,,,,,,
Date,7/30/2020,,,,,,,,
Workflow,GenerateFASTQ,,,,,,,,
Application,FASTQ Only,,,,,,,,
Instrument Type,MiSeq,,,,,,,,
Assay,Nextera DNA,,,,,,,,
Index Adapters,"Nextera Index Kit (24 Indexes, 96 Samples)",,,,,,,,
Description,,,,,,,,,
Chemistry,Amplicon,,,,,,,,
,,,,,,,,,
[Reads],,,,,,,,,
76,,,,,,,,,
76,,,,,,,,,
,,,,,,,,,
[Settings],,,,,,,,,
ReverseComplement,0,,,,,,,,
Adapter,CTGTCTCTTATACACATCT,,,,,,,,
,,,,,,,,,
[Data],,,,,,,,,
Sample_ID,Sample_Name,Sample_Plate,Sample_Well,I7_Index_ID,index,I5_Index_ID,index2,Sample_Project,Description
wenku1,,,,D701,ATTACTCG,D504,GGCTCTGA,,
wenku2,,,,D702,TCCGGAGA,D504,GGCTCTGA,,
wenku3,,,,D703,CGCTCATT,D504,GGCTCTGA,,
wenku4,,,,D704,GAGATTCC,D504,GGCTCTGA,,
wenku5,,,,D705,ATTCAGAA,D504,GGCTCTGA,,

这是一个每行都用逗号,隔开的SampleSheet.csv文件，可以用Excel或WPS软件编辑成新的，也可以用Illumina Experiment Manager软件来制作清单文件。

由于测序仪的通量比较大，如果一次只测一个样本，会造成非常大的浪费，因此实际测序过程中都是多个样本混合在一起测，那么后期怎么把各个样本的数据（即大量的reads）分开呢，这就要涉及到建库的原理了，需要较大的篇幅才能说清，在此不便展开。

简单来说，就是各个样本都加了自己特有的条码（barcode），就是SampleSheet文件中的index和index2序列，如wenku1的index序列ATTACTCG，以及index2序列GGCTCTGA，它们的组合与其他所有文库的都不一样，依据这些条码就能实现数据的拆分，[Data]部分一行是一个样本，每行最低限度只需要填写Sample_ID和index就可以了，如果是双端index测序，再填上index2即可。

--output-dir参数，结果文件保存的目录。

--barcode-mismatches参数，拆分时允许上述index序列错配的碱基个数，通常设为1。

--use-bases-mask参数，用于指导哪些测序循环的数据要保留，哪些要丢弃，示例中Y75n,I8,I8,Y75n的含义是：第一轮测序的前75个循环要保留，之后的循环丢弃，index的8个循环全部保留，index2的8个循环全部保留，第四轮测序的前75个循环要保留，之后的循环丢弃，这样我们拆分出的数据是双端75bp，当然这个参数并非必需。

以上，就是illumina数据拆分的基本过程，结果得到FASTQ文件（每个样本2个），后续分析都是基于这些文件的。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11


Demultiplexed/
├── wenku1_S1_L001_R1_001.fastq.gz
├── wenku1_S1_L001_R2_001.fastq.gz
├── wenku2_S2_L001_R1_001.fastq.gz
├── wenku2_S2_L001_R2_001.fastq.gz
├── wenku3_S3_L001_R1_001.fastq.gz
├── wenku3_S3_L001_R2_001.fastq.gz
├── wenku4_S4_L001_R1_001.fastq.gz
├── wenku4_S4_L001_R2_001.fastq.gz
├── wenku5_S5_L001_R1_001.fastq.gz
└── wenku5_S5_L001_R2_001.fastq.gz

数据分析任重道远，正如莎士比亚的咏叹：

自然啊，你是充满无穷神秘的书！而我，只能读懂些许部分！

文章目录