序列组装

序列组装

一、摘要

  1. 加深序列组装原理的理解;
  2. 熟悉已有的基因组序列组装软件的使用;
  3. 掌握常用的序列组装软件的使用;
  4. 能够独立地使用自己所掌握的编程语言编写序列组装的简化程序。

二、材料和方法

1、硬件平台

处理器:Intel(R) Core(TM)i7-4710MQ CPU @ 2.50GHz
安装内存(RAM):16.0GB

2、系统平台

Windows 8.1、Ubuntu

3、软件平台

【1】AllPathsLG
【2】SOAP
【3】GenomeABC http://crdd.osdd.net/raghava/genomeabc/
【4】Python3.5
【5】Biopython

4、数据库资源

NCBI数据库:https://www.ncbi.nlm.nih.gov/

5、研究对象

噬菌体基因组Dickeya phage RC-2014
https://www.ncbi.nlm.nih.gov/genome/?term=NC_025452

6、方法

序列组装软件的概览
进入pubmed搜索序列组装相关的文献,仔细阅读并尝试安装、使用其中的组装软件。
待拼接序列数据源的准备
使用实验1的模拟测序程序,生成待拼接序列数据源
小规模序列拼接软件的使用练习
GenomeABC 在线平台,对第2 步获得的数据源进行拼接和评估。
大规模基因组序列拼接软件的使用练习
练习序列拼接软件AllPathsLG、soapdenovo2、Velvet
序列组装程序的编写
下载安装python,并且安装biopython扩展模块,编写程序,对测序结果进行组装
不同组装结果的对比
把上述软件和自行编程得到的拼接结果进行比较,分析它们之间的覆盖率,以及contigs 的差异等

三、结果

1、序列组装软件

生物信息学中的一个重要问题就是基因组的组装,NGS技术产生了大量的基因组序列的片段,随之而来的是需要大量内存空间的基因组组装。伴随着最近测序技术的进步,可以预见的是内存空间占用将会急剧增加,并且会成为处理基因组测序数据的限制因素。
常见的序列组装软件有:
SSAKE:将short reads严格地聚类组装成contigs
Ray:主要用于组装Illumina, 454, SOLiD测序生成的reads
a5-miseq:生成高质量的细菌基因组组装。
Orione:是一个包含公共的软件和pipline的框架,针对微生物测序的组装
SGA:基于线性图的组装软件,由于使用压缩表示DNA测序reads的方法,节约内存。
在这里插入图片描述
图表 1相关文献
在这里插入图片描述
图表 2汇总组装软件的网站
文献链接:https://www.ncbi.nlm.nih.gov/pubmed/24086547
网站链接:https://omictools.com/genome-assembly-category
PS:其他常见软件如AllPathsLG、soapdenovo2、Velvet,此网站也有涉及,由于下面将要具体应用,此处不一一详述。
这里以安装SOAP为例
安装SOAP,首先unzip解压
在这里插入图片描述
图表 3解压缩
然后进行make进行编译
在这里插入图片描述
图表 4编译
最后,配置环境变量。
在这里插入图片描述
图表 5配置环境变量
测试其运行情况
在这里插入图片描述
图表 6程序测试

2、待拼接序列数据源的准备

使用之前模拟测序的程序,进行噬菌体基因组的模拟测序,由于后面做的是编程分支,python程序效率有限,因此这里N值仅仅为50,覆盖率为0.7251
在这里插入图片描述
图表 7模拟测序

3、小规模序列拼接软件的使用练习

选择模拟测序生成的virusSingleRead.fa文件,上传并且运行
在这里插入图片描述
图表 8上传测序文件
在这里插入图片描述
图表 9 GenomeABC结果下载

4、大规模基因组序列拼接软件的使用练习

首先,下载安装AllPaths-LG,并且下载测试数据。
接着./configure配置
在这里插入图片描述
图表 10配置
make命令,进行gcc编译,此步骤耗时很长
在这里插入图片描述
图表 11编译
make install进行安装
在这里插入图片描述
图表 12安装
使用测试数据进行组装准备,生成data文件夹
在这里插入图片描述
图表 13组装准备脚本
在这里插入图片描述
图表 14组装准备
接着进行组装
在这里插入图片描述
图表 15组装脚本
在这里插入图片描述
图表 16运行组装
组装结果,data文件夹下面都是运行的结果,总计拼出来一条scaffold
在这里插入图片描述
在这里插入图片描述

5、编程模拟组装结果

计算reads间的overlap
在这里插入图片描述
拼装生成Contig
在这里插入图片描述
测序覆盖度(m值)1.29,噬菌体基因组总长度:155.346kbp
拼接出的contig总长度:188.754kbp
拼接覆盖度:1.22%

6、不同组装结果的对比

特点:

自己编写的程序运行效率会低得多,但是由于之前模拟测序并没有引入测序错误,所以拼接出来的contig全都能完全匹配到原来的基因组。
AllPaths-LG运行速度相对较快,开源软件准确率有保证,但是运行结果有一大堆,如果在实际中,需要仔细阅读相关的介绍文档,才能正确使用

四、讨论和结论

程序运行方法

需要调整文件句柄参数,分别是读入的测序文件句柄、输出的组装文件句柄、读入进行后续比较的参考基因组文件句柄。
在is_overlap 方法中调整minoverlaplength,进而调整Reads间的最小overlap值

附录

from Bio import SeqIO
from Bio.SeqRecord import SeqRecord

readsList = []
for seq_record in SeqIO.parse("data/virusSingleRead.fa", "fasta"):
    readsList.append(seq_record)

# 过滤掉互相包含的reads
readsfilterList = []
for i in range(len(readsList)):
    j = 0
    flag = True
    while j < len(readsList):
        # noinspection PyProtectedMember
        if i != j and readsList[i].seq in readsList[j].seq:
            flag = False
            # print(readsList[i].id + "包含于" + readsList[j].id)
        j += 1
    if flag:
        readsfilterList.append(readsList[i])

del readsList


# noinspection PyProtectedMember
# 分别比较两字符串开头、结尾,计算overlap的大小;最后返回最大的overlap值
def is_overlap(id1, id2):
    minoverlaplength = 10
    score1 = score2 = 0
    reads1 = readsfilterList[id1]
    reads2 = readsfilterList[id2]
    seq1 = reads1.seq
    seq2 = reads2.seq
    minlength = len(seq1) if len(seq1) < len(seq2) else len(seq2)
    for i in range(minoverlaplength, minlength):
        if seq1[-i:] == seq2[:i]:
            print(reads1.id + "\t" + reads2.id + "\t" + str(i) + "\tafter")
            score1 = i
    for j in range(minoverlaplength, minlength):
        if seq1[:j] == seq2[-j:]:
            print(reads1.id + "\t" + reads2.id + "\t" + str(j) + "\tbefore")
            score2 = j
    if score1 == score2 == 0:
        return id1, id2, 0, "NoOverlap"
    if score1 > score2:
        return id1, id2, score1, "After"
    else:
        return id1, id2, score2, "Before"


# noinspection PyShadowingNames
# 构造边的类,存储:有overlap关系的两个reads在readsfilterList中的序号;
# score值(overlap的程度);tips作为辅助;记录哪儿个reads在前,哪儿个reads在后
# PS:tips其实最后没有用到,存储时候就已经前后排序过
class Edge:
    def __init__(self, id1, id2, score, tips):
        self.id1 = id1
        self.id2 = id2
        self.score = score
        self.tips = tips

    def __str__(self):
        return str(self.id1) + " " + str(self.id2) + " " + str(self.score) + " " + self.tips

    def show(self):
        print(self)


edgeList = []
for i in range(len(readsfilterList)):
    j = i + 1
    while j < len(readsfilterList):
        # print(readsfilterList[i].id + " " + readsfilterList[j].id)
        (id1, id2, score, tips) = is_overlap(i, j)
        if tips != "NoOverlap":
            # print(str(id1) + " " + str(id2) + " " + str(score) + " " + tips)
            if tips == "Before":
                edgeList.append(Edge(id2, id1, score, tips))
            if tips == "After":
                edgeList.append(Edge(id1, id2, score, tips))
                # edgeList.append(Edge(id1, id2, score, tips))
        j += 1
edgeList = sorted(edgeList, key=lambda x: x.score, reverse=True)

# 将edge边组装成一条条路径
pathList = [[edgeList.pop(0)]]

# 将edge按序号排列好
for edge in edgeList:
    existflag = True
    for path in pathList:
        if edge.id1 == path[-1].id2 and edge.tips == "After":
            path.append(edge)
            existflag = False
            break
        if edge.id2 == path[0].id1 and edge.tips == "Before":
            path.insert(0, edge)
            existflag = False
            break
    if existflag:
        pathList.append([edge])

# 根据一条条路径,将reads拼到一起,生成contig
contigList = []
ContigID = 1
for path in pathList:
    contig = readsfilterList[path[0].id1].seq
    for edge in path:
        contig = contig + readsfilterList[edge.id2].seq[edge.score:]
    contigList.append(SeqRecord(contig, id="Contig" + str(ContigID), description=""))
    print("Contig" + str(ContigID))
    print(contig)
    ContigID += 1
SeqIO.write(contigList, "result/virusResult.fa", "fasta")

# 对生成的contigs进行测试
record_iterator = SeqIO.parse("data/NC_025452.fasta", "fasta")
chromosome = record_iterator.__next__()
contigLength = 0
for contig in contigList:
    flag = str(contigList[0].seq) in str(chromosome.seq)
    if flag:
        contigLength += len(contig)
        # print(contig.id+" in "+chromosome.id)
    else:
        print(contig.id + " not in " + chromosome.id)
print("基因组总长度:" + str(len(chromosome) / 1000) + "kbp")
print("contig总长度:" + str(contigLength / 1000) + "kbp")
print("覆盖度:" + str(round(contigLength / len(chromosome), 2)) + "%")
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值