一、获取SRA号
登录Genome Announcements网站(https://mra.asm.org/),搜索关键词“bacteria genome SRA”,在搜索到的细菌基因组文章中选择一篇,找到文章记载的SRA号。
以下面文章中的SRR号 SRR6466501 为例:
二、下载SRA文件
SRA (Sequence ReadArchive)数据库,是用于存储二代测序的原始数据,包括454, Illumina, SOLiD, lonTorrent, Helicos 和CompleteGenomics。除了原始序列数据外,SRA现在也存在raw reads在参考基因的比对信息。
根据SRA数据产生的特点,将SRA数据分为四类:Studies-- 研究课题、Experiments-- 实验设计、Runs-- 测序结果集、Samples-- 样品信息。
SRA文件,SRA数据库用不同的前级加以区分:
●ERP或SRP表示Sudics;
●SRS表示Samplo;
●SRX表示Hxpcritmcrnt;
●SRR表示Runs;
从SRA数据库上用prefetch下载sra文件,输入命令如下:
prefetch SRR6466501 #下载
结果如下:
软件自动建立~/ncbi/public/sra文件夹,查看:
三、 Fastq-dump解压
cd ~/ncbi/public/sra #进入到sra文件夹下
fastq-dump --gzip --split-files SRR6466501.sra #解压
结果如下:
四、Fastqc质控
*准备:Fastqc安装,可参考我的简书https://www.jianshu.com/p/e0659f09288c
使用以下命令进行质控:
fastqc SRR6466501_1.fastq.gz
fastqc SRR6466501_2.fastq.gz
结果:
FastQC报告
打开:SRR6466501_1_fastqc.html
SRR6466501_2_fastqc.html
五、Trimmomatic去接头
切除接头序列和低质量碱基,此处使用数据过滤软件Trimmomatic。
Trimmomatic 是一个广受欢迎的Ilumina平台数据过滤工具。
处理数据速度快,主要用来去除Illumina 平台的Fastq序列中的接头,并根据碱基质量值对Fastq进行修剪。
支持多线程,有两种过滤模式,分别对应SE和PE测序数据。
用法:
单末端测序模式
在SE模式下,只有一个输入文件和一个过滤后的输出文件
java -jar [Trimmomatic所在的绝对路径] SE-phred33 [输入文件] [输出文件] [动作1] [动作2]......
双末端测序模式
在PE模式下,有两个输入文件,正向测序序列和反向测序序列,过滤之后,输出文件有四个两端序列都保留的为paried -端序列过滤后被遗弃另- -端保留为unparied
java jar [Trimmomatic的绝对路径] PE -phred33
第一个输入文件(forward端)正向端
第二个输入文件(reverse端)反向端
输出文件[forward_ paried forward_unpaired reverse_ paired reverse_ _unpaired ] [动作1] [动作2]......
另外也支持phred-33和phred-64格式互相转化,不过现在绝大部分Illumina平台的产出数据都转为使用phred-33格式了。
*准备:Trimmomatic安装:
unzip Trimmomatic-0.38.zip
cd Trimmomatic-0.38
java -jar trimmomatic-0.38.jar
Trimmomatic过滤命令如下:
java -jar /home/gxx/Trimmomatic-0.38/trimmomatic-0.38.jar PE -phred33 SRR6466501_1.fastq.gz SRR6466501_2.fastq.gz ./output_forward_paired.fq.gz ./output_forward_unpaired.fq.gz ./output_reverse_paired.fq.gz ./output_reverse_unpaired.fq.gz ILLUMINACLIP:/home/gxx/Trimmomatic-0.38/adapters/TruSeq2-PE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:20 TRAILING:20 MINLEN:75
结果如下:
动作说明
ILLUMINACLIP:
过滤reads中的Illumina测序街头和引|物序列并决定是否去除反向互补的R1/R2中的R2
参数说明( PE测序要注意最后-一个参数, SE最后两个参数不用设置,参数之间用冒号连接):
<接头和引物序列所在位置( Trimmomatic自带在adapters里面,其中TruSeq2为2代测序, TruSeq3为三代测序)> :
<seed搜索允许多少个个碱基错配)>:<alindrome比对分值阈值为多少>:<simple clip比对分值阈值为多少>:
<palindrome模式允许切除的最短接头序列为多少默认为8bp>:
<palindrome模式去除与R1完全反向的R2>
SLIDINGWINDOW:
从reads的5'端开始,进行滑窗质量过滤,切掉碱基质量平均值低于阈值的滑窗
参数说明: <窗口大小>: <碱基平均质量值阈值>
MAXINFO:
一个自动调整的过滤选项,在保证reads长度的情况下尽量降低测序错误率,最大化reads的使用价值
LEADING:
从reads的开头切除质量低于阈值的碱基
TRAILING:
从reads的末尾切除质量低于阈值的碱基
CROP:
从reads的末尾切掉部分碱基是reads达到指定长度
HEADCROP:
从reads的开头切掉指定数量的碱基
MINLEN:
如果经过剪切后reads的长度低于阈值则丢弃这条reads
AVGQUAL:
如果reads的平均碱基质量值低于阈值则丢弃这条reads
TOPHRED33:
将reads的碱基质量值体系转为phred-33
TOPHRED64:
将reads的碱基质量值体系转化为phred-64
六、SPAdes组装基因组草图
SPAdes:
➢由俄罗斯科学院圣彼得堡理工大学计算生物学实验室开发,是目前评价最好的拼接工具之一。
➢主要用于基因组拼接,也可用于一、二、三代测序的混合组装;还可用于转录组从头组装(rnaSPAdes)和宏基因组拼接(metaSPAdes) 。
*准备:SPAdes的安装
python环境下,安装命令如下:
wget http://cab.spbu.ru/files/release3.12.0/SPAdes-3.12.0-Linux.tar.gz
mkdir ~/Biosofts/spades
tar zvxf /home/gxx/SPAdes-3.12.0-Linux.tar.gz -C ~/Biosofts/spades/
~/Biosofts/spades/SPAdes-3.12.0-Linux/bin/spades.py -h
echo 'export PATH=~/Biosofts/spades/SPAdes-3.12.0-Linux/bin:$PATH' >> ~/.bashrc
source ~/.bashrc
spades.py -h
SPAdes运行:
spades.py --careful --pe1-1 SRR6466501_1.fastq.gz --pe1-2 SRR6466501_2.fastq.gz -o ./SPAdesout_SRR6466501new
遇到错误:out of memory,内存不足。
解决:将虚拟机设置中的系统内存调至合理范围内最大,再次尝试以上命令,但最终还是不行。。
尝试使用seqtk抽取1000条,命令如下:
#解压
gunzip -c output_forward_paired.fq.gz >output_forward_paired.fq
gunzip -c output_reverse_paired.fq.gz >output_reverse_paired.fq
#抽取1000条
seqtk sample -s 60 output_forward_paired.fq 1000 >seqtksample1_1000.fq
seqtk sample -s 60 output_reverse_paired.fq 1000 >seqtksample2_1000.fq
#用wc查看,可对比前后文件,判断是否抽取成功
wc -l output_forward_paired.fq
wc -l seqtksample1_1000.fq
然后,再次尝试SPAdes运行:
spades.py --careful --pe1-1 seqtksample1_1000.fq --pe1-2 seqtksample2_1000.fq -o ./SPAdesout_SRR6466501_1000new
结果如下:
七、Quast评价组装的基因组效果
*准备:Quast的安装,安装命令如下
tar zvxf quast-5.0.0.tar.gz -C ~/Biosofts/
~/Biosofts/quast-5.0.0
python quast.py -h
运行代码:
quast.py ~/ncbi/public/sra/SPAdesout_SRR6466501_1000new/contigs.fasta -o ~/ncbi/public/sra/SPAdesout_SRR6466501_1000new/quast_out
结果如下:
查看输出的文件夹quast_out:
本地下载quast报告 report.html,并查看: