细菌基因组的组装

一、获取SRA号

登录Genome Announcements网站(https://mra.asm.org/),搜索关键词“bacteria genome SRA”,在搜索到的细菌基因组文章中选择一篇,找到文章记载的SRA号。

以下面文章中的SRR号 SRR6466501 为例:


SRR6466501.JPG

二、下载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,并查看:

完成!!!

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 203,547评论 6 477
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,399评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 150,428评论 0 337
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,599评论 1 274
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,612评论 5 365
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,577评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,941评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,603评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,852评论 1 297
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,605评论 2 321
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,693评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,375评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,955评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,936评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,172评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 43,970评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,414评论 2 342

推荐阅读更多精彩内容