WGS数据分析流程学习与开发过程全纪录（2）

最新推荐文章于 2024-06-11 11:28:56 发布

littleblack201608

最新推荐文章于 2024-06-11 11:28:56 发布

阅读量1.4k

点赞数 1

分类专栏：学习总结文章标签： WGS 数据质控

本文链接：https://blog.csdn.net/littleblack201608/article/details/87270581

版权

学习总结专栏收录该内容

4 篇文章 2 订阅

订阅专栏

上一篇主要是跟着别的教程走了一遍大肠杆菌的数据分析流程，这一篇开始就一块一块地具体学习。

本篇主要学习的是数据质控。

1. FastQC

a. 下载并安装、配置FastQC

(wgs) [xuquan@zhswlcserver:conda]$cd ~/conda/softwares/
(wgs) [xuquan@zhswlcserver:softwares]$wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.8.zip
(wgs) [xuquan@zhswlcserver:softwares]$unzip fastqc_v0.11.8.zip
(wgs) [xuquan@zhswlcserver:softwares]$cd FastQC/
(wgs) [xuquan@zhswlcserver:FastQC]$chmod 755 fastqc
(wgs) [xuquan@zhswlcserver:FastQC]$echo "export PATH=/home/xuquan/conda/softwares/FastQC:$PATH" >> ~/.bashrc
(wgs) [xuquan@zhswlcserver:FastQC]$source ~/.bashrc

至此，已经可以在系统的任何工作目录下直接使用“fastqc”调用FastQC这个软件了！

b. 测试使用FastQC

为规范文件路径，在系统的存储硬盘上创建一个名为condaHomeConda的目录，并在其下分别创建如下三个文件夹：
databases
output
rawdata
分别用来存储数据库文件、输出文件以及原始文件。
其中，output以及rawdata下面又分别根据不同的分析项目，创建相应的子目录。
然后，将该名为condaHomeConda的目录软连接到/home/xuquan/conda路径下，这样就能通过“/home/xuquan/conda/condaHomeConda”看到这三个文件夹了：

(wgs) [xuquan@zhswlcserver:condaHomeConda]$pwd
/home/xuquan/conda/condaHomeConda
(wgs) [xuquan@zhswlcserver:condaHomeConda]$ll
总用量 2
drwxrwsr-x 1 xuquan sx04 4096 2月  14 12:03 databases
drwxrwsr-x 1 xuquan sx04 4096 2月  14 12:04 output
drwxrwsr-x 1 xuquan sx04 4096 2月  14 12:08 rawdata

接着，使用前一篇文章中下载的E.coli的测序数据进行质控练习：

(wgs) [xuquan@zhswlcserver:condaHomeConda]$cd output/
(wgs) [xuquan@zhswlcserver:output]$ll
总用量 0
(wgs) [xuquan@zhswlcserver:output]$mkdir testFastqc01
(wgs) [xuquan@zhswlcserver:output]$cd testFastqc01/
(wgs) [xuquan@zhswlcserver:testFastqc01]$fastqc /home/xuquan/conda/rawdata/SRR1770413_* -o /home/xuquan/conda/condaHomeConda/output/testFastqc01/
Started analysis of SRR1770413_1.fastq.gz
Approx 5% complete for SRR1770413_1.fastq.gz
Approx 10% complete for SRR1770413_1.fastq.gz
Approx 15% complete for SRR1770413_1.fastq.gz
Approx 20% complete for SRR1770413_1.fastq.gz
Approx 25% complete for SRR1770413_1.fastq.gz
Approx 30% complete for SRR1770413_1.fastq.gz
Approx 35% complete for SRR1770413_1.fastq.gz
Approx 40% complete for SRR1770413_1.fastq.gz
Approx 45% complete for SRR1770413_1.fastq.gz
Approx 50% complete for SRR1770413_1.fastq.gz
Approx 55% complete for SRR1770413_1.fastq.gz
Approx 60% complete for SRR1770413_1.fastq.gz
Approx 65% complete for SRR1770413_1.fastq.gz
Approx 70% complete for SRR1770413_1.fastq.gz
Approx 75% complete for SRR1770413_1.fastq.gz
Approx 80% complete for SRR1770413_1.fastq.gz
Approx 85% complete for SRR1770413_1.fastq.gz
Approx 90% complete for SRR1770413_1.fastq.gz
Approx 95% complete for SRR1770413_1.fastq.gz
Analysis complete for SRR1770413_1.fastq.gz
Started analysis of SRR1770413_2.fastq.gz
Approx 5% complete for SRR1770413_2.fastq.gz
Approx 10% complete for SRR1770413_2.fastq.gz
Approx 15% complete for SRR1770413_2.fastq.gz
Approx 20% complete for SRR1770413_2.fastq.gz
Approx 25% complete for SRR1770413_2.fastq.gz
Approx 30% complete for SRR1770413_2.fastq.gz
Approx 35% complete for SRR1770413_2.fastq.gz
Approx 40% complete for SRR1770413_2.fastq.gz
Approx 45% complete for SRR1770413_2.fastq.gz
Approx 50% complete for SRR1770413_2.fastq.gz
Approx 55% complete for SRR1770413_2.fastq.gz
Approx 60% complete for SRR1770413_2.fastq.gz
Approx 65% complete for SRR1770413_2.fastq.gz
Approx 70% complete for SRR1770413_2.fastq.gz
Approx 75% complete for SRR1770413_2.fastq.gz
Approx 80% complete for SRR1770413_2.fastq.gz
Approx 85% complete for SRR1770413_2.fastq.gz
Approx 90% complete for SRR1770413_2.fastq.gz
Approx 95% complete for SRR1770413_2.fastq.gz
Analysis complete for SRR1770413_2.fastq.gz
(wgs) [xuquan@zhswlcserver:testFastqc01]$ll -lh
总用量 903K
-rw-rw-r-- 1 xuquan sx04 218K 2月  14 2019 SRR1770413_1_fastqc.html
-rw-rw-r-- 1 xuquan sx04 224K 2月  14 2019 SRR1770413_1_fastqc.zip
-rw-rw-r-- 1 xuquan sx04 225K 2月  14 2019 SRR1770413_2_fastqc.html
-rw-rw-r-- 1 xuquan sx04 236K 2月  14 2019 SRR1770413_2_fastqc.zip

得到的read1和read2的“Per base sequence quality”图如下所示：
read1 read2 除了这个图之外，还有其他的一些图，这里不一一上贴。FastQC的主要作用是为了让你认识你的测序数据，但是其本身并不能够对数据进行任何的处理操作。后续还需要根据你所了解到的数据的质量，通过其他质控软件进行相应的数据质控。

当我们理解了fq数据之后，做这些过滤就不会很难，你也完全可以自己编写工具来进行个性化的过滤。目前也已有很多工具用来切除接头序列和低质量碱基，比如SOAPnuke、cutadapt、untrimmed等不下十个，但这其中比较方便好用的是Trimmomatic（也是一个java程序）、sickle和seqtk。Trimmomatic的好处在于，它不但可以用来***切除illumina测序平台的接头序列***，还可以去除由我们自己指定的***特定接头序列***，而且同时也能够***过滤read末尾的低质量序列***，sickle和seqtk只能去除低质量碱基。具体的原理就是通过***滑动一定长度的窗口***，计算窗口内的碱基平均质量，如果过低，就直接往后全部切除，注！意！不是挖掉read中的这部分低质量序列，而是像切菜一样，直接从低质量区域开始把这条read后面的所有其它碱基***全！部！剁！掉！***否则就是在人为改变实际的基因组序列情况。

如果下机的fq数据中不含有这些测序接头，那么我们除了trimmomatic之外，也可以直接使用sickle（同时支持PE和SE数据）或者seqtk（仅支持SE），这两个处理起来会更快，消耗的计算资源也更少。

2. Trimmomatic

a. 下载并安装、配置Trimmomatic

(wgs) [xuquan@zhswlcserver:~]$cd ~/conda/softwares/
(wgs) [xuquan@zhswlcserver:softwares]$wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.38.zip
(wgs) [xuquan@zhswlcserver:softwares]$unzip Trimmomatic-0.38.zip
(wgs) [xuquan@zhswlcserver:softwares]$cd Trimmomatic-0.38/
(wgs) [xuquan@zhswlcserver:Trimmomatic-0.38]$ll
总用量 168
drwxr-xr-x 2 xuquan xuquan   4096 5月  16 2018 adapters
-rw-r--r-- 1 xuquan xuquan  35147 5月  16 2018 LICENSE
-rw-r--r-- 1 xuquan xuquan 127618 5月  16 2018 trimmomatic-0.38.jar

目录下的“trimmomatic-0.38.jar”就是执行程序，可以直接使用java来运行：

(wgs) [xuquan@zhswlcserver:Trimmomatic-0.38]$java -jar trimmomatic-0.38.jar
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

b. 尝试使用Trimmomatic

下面用上一篇文章中下载的大肠杆菌测试数据来测试使用一下该软件：

(wgs) [xuquan@zhswlcserver:adapters]$java -jar /home/xuquan/conda/softwares/Trimmomatic-0.38/trimmomatic-0.38.jar PE -phred33 -trimlog /home/xuquan/conda/output/trimmomatic.logfile /home/xuquan/conda/rawdata/SRR1770413_1.fastq.gz /home/xuquan/conda/rawdata/SRR1770413_2.fastq.gz /home/xuquan/conda/rawdata/out.SRR1770413_1.fastq.gz /home/xuquan/conda/rawdata/trim.SRR1770413_1.fastq.gz /home/xuquan/conda/rawdata/out.SRR1770413_2.fastq.gz /home/xuquan/conda/rawdata/trim.SRR1770413_2.fastq.gz ILLUMINACLIP:/home/xuquan/conda/softwares/Trimmomatic-0.38/adapters/TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:50
TrimmomaticPE: Started with arguments:
 -phred33 -trimlog logfile /home/xuquan/conda/rawdata/SRR1770413_1.fastq.gz /home/xuquan/conda/rawdata/SRR1770413_2.fastq.gz /home/xuquan/conda/rawdata/out.SRR1770413_1.fastq.gz /home/xuquan/conda/rawdata/trim.SRR1770413_1.fastq.gz /home/xuquan/conda/rawdata/out.SRR1770413_2.fastq.gz /home/xuquan/conda/rawdata/trim.SRR1770413_2.fastq.gz ILLUMINACLIP:/home/xuquan/conda/softwares/Trimmomatic-0.38/adapters/TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:50
Using PrefixPair: 'TACACTCTTTCCCTACACGACGCTCTTCCGATCT' and 'GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT'
ILLUMINACLIP: Using 1 prefix pairs, 0 forward/reverse sequences, 0 forward only sequences, 0 reverse only sequences
Input Read Pairs: 643253 Both Surviving: 599379 (93.18%) Forward Only Surviving: 36513 (5.68%) Reverse Only Surviving: 1485 (0.23%) Dropped: 5876 (0.91%)
TrimmomaticPE: Completed successfully

同时需要明确指明质量值体系是Phred33还是Phred64，默认是Phred64，这需要特别注意，因为我们现在的测序数据基本都是Phred33的了，所以一定要指定这个参数。

OK，执行成功，可以在指定的输出路径下看到两端分别被丢弃的和保留下来的序列：

(wgs) [xuquan@zhswlcserver:adapters]$cd /home/xuquan/conda/rawdata/
(wgs) [xuquan@zhswlcserver:rawdata]$ll -lh
总用量 629M
-rw-rw-r-- 1 xuquan xuquan  95M 2月  16 10:38 out.SRR1770413_1.fastq.gz
-rw-rw-r-- 1 xuquan xuquan  87M 2月  16 10:38 out.SRR1770413_2.fastq.gz
-rw-rw-r-- 1 xuquan xuquan 120M 2月  14 09:34 SRR1770413_1.fastq.gz
-rw-rw-r-- 1 xuquan xuquan 134M 2月  14 09:36 SRR1770413_2.fastq.gz
-rw-rw-r-- 1 xuquan xuquan 190M 2月  14 09:02 SRR1770413.sra
-rw-rw-r-- 1 xuquan xuquan 5.3M 2月  16 10:38 trim.SRR1770413_1.fastq.gz
-rw-rw-r-- 1 xuquan xuquan 173K 2月  16 10:38 trim.SRR1770413_2.fastq.gz

质控后的数据再用fastqc来看一下：

(wgs) [xuquan@zhswlcserver:rawdata]$fastqc /home/xuquan/conda/rawdata/out.SRR1770413_* -o /home/xuquan/conda/condaHomeConda/output/testFastqc02

结果分别如下：

reads1 reads2 对比质控之前，差别还是蛮大的！

c. Trimmomatic参数说明

ILLUMINACLIP，接头序列切除参数。ILLUMINACLIP:/home/xuquan/conda/softwares/Trimmomatic-0.38/adapters/TruSeq3-PE.fa:2:30:10 意思分别是：
1.1 TruSeq3-PE.fa是接头序列；
1.2 2是比对时接头序列时所允许的最大错配数；
1.3 30指的是要求PE的两条read同时和PE的adapter序列比对，匹配度加起来超30%，那么就认为这对PE的read含有adapter，并在对应的位置需要进行切除【注】。
1.4 10和前面的30不同，它指的是，我就什么也不管，反正只要这条read的某部分和adpater序列有超过10%的匹配率，那么就代表含有adapter了，需要进行去除；

【注】测序的时候一般只会测到一部分的adapter，因此read和adaper对比的时候肯定是不需要要求百分百匹配率的，上述30%和10%其实是比较推荐的值。

SLIDINGWINDOW，滑动窗口长度的参数，SLIDINGWINDOW:5:20代表窗口长度为5，窗口中的平均质量值至少为20，否则会开始切除；
LEADING，规定read开头的碱基是否要被切除的质量阈值；
TRAILING，规定read末尾的碱基是否要被切除的质量阈值；
MINLEN，规定read被切除后至少需要保留的长度，如果低于该长度，会被丢掉。

参考资料：从零开始完整学习全基因组测序（WGS）数据分析：第3节数据质控

littleblack201608

关注

1
点赞
踩
10

收藏

觉得还不错? 一键收藏
0
评论
WGS数据分析流程学习与开发过程全纪录（2）

上一篇主要是跟着别的教程走了一遍大肠杆菌的数据分析流程，这一篇开始就一块一块地具体学习。本篇主要学习的是数据质控。1. FastQCa. 下载并安装、配置FastQC(wgs) [xuquan@zhswlcserver:conda]$cd ~/conda/softwares/(wgs) [xuquan@zhswlcserver:softwares]$wget http://www.bio...
复制链接

扫一扫

专栏目录