基因组组装（全）

2019-12-06 本文已影响0人麋鹿吃了颗草莓

基因组组装一般分为三个层次，contig, scaffold和chromosomes.

contig表示从大规模测序得到的短读(reads)中找到的一致性序列。组装的第一步就是从短片段(pair-end)文库中组装出contig。进一步基于不同长度的大片段(mate-pair)文库，将原本孤立的contig按序前后连接，这一步会得到scaffolds。最后基于遗传图谱或光学图谱将scaffold合并调整，形成染色体级别的组装(chromosome)

一. 下载短序列

首先到Microbiology Resource Annocements（https://mra.asm.org）上找到需要下载reads的SRA号，比如我们找到两篇文章中的SRA号，分别为 SRR020180 和 SRR028694
prefetch下载序列：下载的SRA文件默认保存在/ncbi/public/sra中

prefetch SRR020180
prefetch  SRR028694

结果

1.PNG

5.PNG

6.PNG

fasterq -dump解压sra文件，将sra文件转化为fastq文件
也可以用fastq -dump命令，但相对于fasterq -dump，fastq -dump的速度太慢了

--split-spot: 将双端测序分为两份,但是都放在同一个文件中
--split-files: 将双端测序分为两份,放在不同的文件,但是对于一方有而一方没有的reads直接丢弃
--split-3 : 将双端测序分为两份,放在不同的文件,但是对于一方有而一方没有的reads会单独放在一个文件夹里

fastq-dump --gzip --split-3 SRR020180.sra
fastq-dump --gzip --split-3 SRR028694.sra

结果：看到两个SRA文件都分别只生成了一个文件，所以两个SRA文件都是单端测序的结果

1.PNG

2.PNG

二. Fastqc质控

FastQC可以快速地对测序数据进行质量评估

输入fastqc -h可以查看fastqc的基本使用参数

10.PNG

-o --outdir 生成的报告文件的存储路径
--(no)extract 是否将生成的报告打包成一个压缩文件
--c contaminant file 污染序列选项
-t --threads 选择程序运行的线程数
-q --quiet 安静运行模式，不设置这个参数时，程序实时报告运行状况

对fastq文件进行质控检验

fastqc SRR020180.fastq.gz
fastqc SRR028694.fastq.gz

结果

11.PNG

红色：数据质量很差
黄色：数据质量一般
绿色：数据质量很好

SRR020180:

12.PNG

SRR028694:

13.PNG

可以看到第二条序列的质量很差，我们接下来需要进行数据的过滤

三. Trimmomatic数据过滤

Trimmomatic 支持多线程，处理数据速度快，主要用来去除 Illumina 平台的 Fastq 序列中的接头，并根据碱基质量值对 Fastq 进行修剪。

运行命令行查看Trimmomatic的使用方法

java -jar ~/Trimmomatic-0.38/trimmomatic-0.38.jar

3.PNG

Trimmomatic有两种过滤模式，分别对应 SE 和 PE 测序数据。SE指单末端测序模式，过滤单端测序产生的数据，PE指双末端测序模式，过滤双端测序产生的数据。

之前在fast-dump解压后两个SAR文件都分别只有一个输出文件，所以2个序列都为单端测序产生的序列，所以我们这里选择SE模式。

 java -jar ~/Trimmomatic-0.38/trimmomatic-0.38.jar SE -phred33 SRR020180.fastq.gz SRR020180_clean.fastq.gz ILLUMINACLIP:/Trimmomatic-0.38/adaptersTruSeq3-SE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:50
java -jar ~/Trimmomatic-0.38/trimmomatic-0.38.jar SE -phred33 SRR028694.fastq.gz SRR028694_clean.fastq.gz ILLUMINACLIP:/Trimmomatic-0.38/adaptersTruSeq3-SE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:50

参数：
—— phred33：指质量值体系为phred33，还有phred64，如果不设置默认为phred64，但因为现在基本都用phred33的了，所以这里一定要设置

—— SRR020180.fastq.gz：输入文件
—— SRR020180_clean.fastq.gz：输出文件

—— ILLUMINACLIP：过滤 reads 中的 Illumina 测序接头和引物序列。TruSeq3-PE.fa是接头序列，2是比对时接头序列时所允许的最大错配数；30指的是要求PE的两条read同时和PE的adapter序列比对，匹配度加起来超30%，那么就认为这对PE的read含有adapter，并在对应的位置需要进行切除。
—— SLIDINGWINDOW：滑动窗口长度的参数，SLIDINGWINDOW:5:20代表窗口长度为5，窗口中的平均质量值至少为20，否则会开始切除；
—— LEADING：规定read开头的碱基是否要被切除的质量阈值；
—— TRAILING：规定read末尾的碱基是否要被切除的质量阈值；
—— MINLEN：规定read被切除后至少需要保留的长度，如果低于该长度，会被丢掉。

数据过滤后我们再用Fastqc质控，查看过滤后的结果如何

fastqc SRR020180_clean.fastq.gz
fastqc SRR028694_clean.fastq.gz

结果：
SRR020180：
数据过滤前：

5.PNG

数据过滤后：

1.png

SRR028694：
数据过滤前：

index.png
数据过滤后：

index.png

四. SPAdes短序列拼接

SPAdes 主要用于进行单细胞测序的细菌与基因组拼接，也能用于非单细胞测序数据。现在的SPAdes版本基本都支持paired-end reads,mate pairs和unpairede reads,多个paired-end和mate pairs可以同时输入。

 spades.py --only-assembler --phred-offset 33 -k 55 --s1 SRR020180_clean.fastq.gz -o ./SPAdes2
 spades.py --only-assembler  --careful --phred-offset 33 -k 33,55,77 --s1 SRR028694_clean.fastq.gz -o ./SPAdes1

参数：
—— only-error-correction：只做数据纠错
—— only-assembler：只组装，不做数据纠错
—— careful：减少错误和插入缺失，添加此选项，会消耗更多的时间

—— phred-offset 33：phred质量体系，在数据纠错中会用到，现在illumina数据一般采用phred 33，并且我们之前数据过滤时采用的就是phred 33

—— k：k值，一次可以输入多个，用逗号分隔，kmer最大为127，并且注意只能是奇数，不设置时会自动计算合适的k值，但运算时间较长

——s1：表明是single reades
——pel：表明是paired-end和mate-pair reades

—— o：输出目录

结果

SRR020180：

6.PNG

SRR028694：

7.PNG

五.Quast评价序列拼接结果

评价结果

#评价SRR020180序列结果
 quast.py ~/ncbi/public/sra/SPAdes2/contigs.fasta
#评价SRR028694序列结果
 quast.py ~/ncbi/public/sra/SPAdes1/contigs.fasta

结果
SRR020180：

9.PNG

SRR028694：

8.PNG

这个序列文件本来就很小，经过了数据过滤，数据纠错等等数据优化后已经没有剩下可以连接的contig了