作业记录：在Linux操作系统上利用各类生物软件练习组装细菌基因组
使用到的一系列工具（提前安装）：SRA-toolkit，Fastqc，trimmomatic，SPAdes，quast
(一) 找到一篇细菌基因组文章及其记载的SRA号
①进入美国微生物所官网，我们需要下载SRA测序数据，选择基因序列的文章分类寻找文章
我们找到了文章：
Draft Genome Sequences of Two Cyanobacteria Leptolyngbya spp. Isolated from Microbial Mats in Miravalles Thermal Spring, Costa Rica | Microbiology Resource Announcements (asm.org)
②在文章中找到数据的SRA号
(二) 利用SRA-toolkit中的prefetch下载SRA文件
需要时可批量下载，在这里作为练习简单下载两个。
prefetch SRR14062498 SRR14062499
(三) 利用Fastqc中的fastq-dump解压数据
fastq-dump可将sra文件解压成不同类型的文件，在这里解压为gz文件以节省空间。
fastq-dump --gzip --split-files SRR14062498/SRR14062498.sra
fastq-dump --gzip --split-files SRR14062499/SRR14062499.sra
#--split-files将解压结果生成新的文件夹，后接文件路径
(四) 利用fastqc数据质量评价
使用fastqc和trimmomatic进行常规的质控和过滤数据，fastqc可以对测序数据的质量进行多方面评价，并输出为html网页格式和zip压缩格式。
mkdir fastqcdata #建立文件夹
fastqc SRR14062498_1.fastq.gz SRR14062498_2.fastq.gz SRR14062499_1.fastq.gz SRR14062499_2.fastq.gz -o fastqcdata
# -o 输出至已存在的文件夹，若文件夹不存在，此选项不会帮你建立新的文件夹
在windows上打开html网页文件可直接查看评测结果，查看基础信息与评测出的其他信息，这里以SRR14062498_1为例作简单说明。
绿色的√表示通过，黄色的！警告和红色的×表示不合格，详细判断信息可参考：测序数据质量控制之FastQC
①× Per base sequence content：
此数据中碱基的内容被判为不合格，理论情况下的AT和GC线应该一致，测序机器开启或关闭时受操作或各种状态影响而不够稳定，这种情况下虽然整体质量通过，想要获取高质量序列时需要对首尾进行一定的裁剪。但注意，裁剪得越多不一定质量越高，因为裁剪意味着失去信息，对于未知序列来说也可能裁剪越严格越好。
②! Per sequence GC content：
GC比值与预测值相差较远而被系统警告。
③! Sequence Length Distribution：
测序长度出的长度不同，被系统警告。
fastqc的质控报告不一定符合需求，需要按照不同的情况判断。对于例子中的嗜热蓝藻菌序列来说，GC含量较高正常，因其结构简单，基因组测序草图总数少（7M和15M），受各因素影响，统计数据波动明显，很正常。
（五）用trimmomatic进行质控
下面根据原文献的处理方法使用trimmomatic进行一定的数据过滤。
mkdir trim_out #创建文件夹
java -jar ~/Biosofts/Trimmomatic-0.38/trimmomatic-0.38.jar PE -phred33 ./fastqcdata/SRR14062498_1_fastqc.zip ./fastqcdata/SRR14062498_2_fastqc.zip ./trim_out/output_forward_paired.fq.gz ./trim_out/output_forward_unpaired.fq.gz ./trim_out/output_reverse_paired.fq.gz ./trim_out/output_reverse_unpaired.fq.gz ILLUMINACLIP:/disk/201931107010248/Biosofts/Trimmomatic-0.38/adapters/TruSeq2-PE.fa:2:30:10 SLIDINGWINDOW:4:20 MINLEN:100
#第一个路径为jar文件路径，第二个和第三个路径为待拼接序列路径，再往后的四个路径为结果输出路径，最后为参数选项
#若找不到文件或出错，检查名称，尝试把路径改为绝对路径
一些参数的含义：
ILLUMINACLIP:接头文件路径:2:30:10 在比对接头序列时允许有2个位置的碱基发生错配，双端测序的两条reads与接头序列匹配率超过30%的话，就会被切除掉，单条reads如果与接头序列的匹配率超过10%，也会被切除掉
SLIDINGWINDOW:4:20 以4bp为窗口进行滑窗统计，切除碱基平均质量低于20的窗口及之后的序列
MINLEN:100 序列片段小于100bp则丢掉这段序列
更多使用信息可参照：NGS 数据过滤之 Trimmomatic 详细说明
生成四个输出结果，分别是正向配对、正向未配对、反向配对和反向未配对。
(六) 用SPAdes组装基因组草图
使用SPAdes软件进行基因组草图的组装，用得到的正向配对和反向配对结果进行序列拼接。
mkdir spades_test #创建存放数据的文件夹
spades.py --careful --phred-offset 33 --pe1-1 ./trim_out/output_forward_paired.fq.gz --pe1-2 ./trim_out/output_reverse_paired.fq.gz -o ./spades_test
一些参数的含义：
--carefull 一种拼接模式，减少错误和插入序列
--phred-offset 33 一种碱基质量体系，在用trimmomatic进行数据过滤时也设置过
--pe1-1为第一条序列信息（正链），后接相对路径，pe1-2为第二条（反链）
-o 后接输出路径
更多使用信息可参照：使用 SPAdes 进行基因组组装
运行时出现了未知错误，询问老师过后，可能是共享服务器的内存不足，我决定将拼接移到我的个人虚拟机上运行（之前在老师提供的SSH运行）。
(七) 用quast评测组装的基因组信息
①quast软件需要的环境版本为python2，首先检查python版本
python --version
若python版本不符合，则需切换到python2环境：
（如果用conda安装，在安装时就会提示环境不同且安装失败。需要切换好python2环境后再安装）
conda create --name python27 python=2.7 -y #下载并命名python2环境
conda activate python27 #加载环境
②在相应环境下，执行quast.py指令，结果输出至quast_out文件夹，可打开查看
mkdir quast_out #创建文件夹
quast.py ./spades_out/contigs.fasta -o quast_out
不同文件格式的报告各有特点，如.txt文本格式纯文本，.pdf图片格式纯图表，.html网页格式文本图表兼并等，可按需查看取用。
基因组基本组装完成，在这里可以查看序列连续contigs、基因片段长度length等总和或连续的信息。其中较为重要的参考值为N50，即由长到短拼接序列并将长度进行累加，当拼接长度达到总序列的一半时，最后拼上的片段长。N50越大，代表长片段的reads更多，拼接的效果越好。
