初探基因组组装——生信原理第四次实验报告

news/2025/2/13 21:49:09/

初探基因组组装——生信原理第四次实验报告

文章目录

  • 初探基因组组装——生信原理第四次实验报告
    • 实验目的
    • 实验内容
    • 实验题目
      • 第一题
        • 题目
        • 用SOAPdenovo 进行基因组组装
        • 评估组装质量
      • 第二题
        • 题目
        • Canu组装
        • Hifiasm组装
        • 基于nucmer的基因组比对
        • 过滤比对结果
        • 转换为可读性强的tab键分隔的文件, -H去除标题行
        • 通过鉴定SNP和Indel,识别两种组装结果的差异
        • 统计SNP和Indel的数量
    • 讨论
      • 调整参数
      • 基因组组装的连续性
      • 组装差异
      • 提高空间

实验目的

1.回顾Linux系统的常用命令的使用

2.掌握常用三代和二代测序组装软件各至少一种的使用,并理解关键参数的含义,熟悉测序数据fastq等格式

3.会编写程序计算N50/L50等组装连续性指标

4.会使用基因组比对工具MUMmer进行序列比对,并寻找SNP等变异

实验内容

1.使用组装软件SOAPdenovo、canu、Hifiasm分别组装大肠杆菌Escherichia coli K12基因组的二代和三代测序数据

2.编写程序计算SOAPdenovo组装contig和scaffold序列的N50/L50等组装质量评估指标

3.使用基因组比对工具MUMmer比较canu组装的contig序列和hifiasm组装的contig序列,并寻找两者之间的序列差异

实验题目

第一题

题目

首先,使用SOAPdenovo软件,组装E.coli基因组二代Illumina测序数据。然后,使用Perl/Python/C++/C/Java…任意一种编程语言编写程序,统计你组装好的E.coli基因组contigscaffold序列的 N50N90L50L90总的碱基数总的序列数目最长序列的长度。注意:contigscaffold序列中长度小于200bp的不要统计在内。

用SOAPdenovo 进行基因组组装

创建好文件夹之后用以下命令进行组装。

SOAPdenovo-63mer all  -s /home/uu01/data/ecoli.cfg  -K 51 -R -o ecoli-soap &

在这里插入图片描述

图1 SOAPdenovo结果

.contig为组装的contig序列,fasta格式

.scafSeq为组装的scaffold序列,fasta格式

评估组装质量

将组装好的E.coli基因组contig和scaffold序列传输到本地,用python计算相应指标。代码如下:

import pandas as pd#两个文件选一个
file_path = "D:/00大三上/生信原理/Data/ecoli-soap.contig"
file_path = "D:/00大三上/生信原理/Data/ecoli-soap.scafSeq"with open(file_path) as fa:fa_dict = {}for line in fa:# 去除末尾换行符line = line.replace('\n', '')if line.startswith('>'):# 去除 > 号seq_name = line[1:]fa_dict[seq_name] = ''else:# 去除末尾换行符并连接多行序列fa_dict[seq_name] += line.replace('\n', '')for name, seq in fa_dict.items():fa_dict[name] = len(seq)fa_df = pd.DataFrame.from_dict(fa_dict, orient='index', columns=['Length'])
# 这一行是筛选scaffold序列用的
# fa_df = fa_df[fa_df.index.str.contains('scaffold')]
fa_df = fa_df[fa_df.Length >= 200]
fa_df.sort_values(by='Length',ascending=False, axis=0, inplace=True)TotalLength=fa_df['Length'].sum()
print(TotalLength)
SeqNumber = fa_df.shape[0]
print(SeqNumber)
MaxLength=fa_df['Length'][0]
print(MaxLength)sum=0
for i in range(SeqNumber):sum+=fa_df.Length[i]if sum/TotalLength>0.5:N50=fa_df.Length[i]L50=i+1break
print(N50)
print(L50)
sum=0
for i in range(SeqNumber):sum+=fa_df.Length[i]if sum/TotalLength>0.9:N90=fa_df.Length[i]L90=i+1break
print(N90)
print(L90)

最终结果如下

评估指标ContigScaffold
N50736110687
N9029227697
L50188014
L90573349
总碱基数44912834890792
总序列数目7603107
最长序列的长度4660326711

第二题

题目

首先,使用canu软件, 组装E.coli基因组的Nanopore测序数据;然后,使用hifiasm软件,组装E.coli基因组的PacBio HiFi测序数据。最后,利用MUMmer软件包,比较canu组装的contig序列和hifiasm组装的contig序列两者之间的差异(需要统计有多少SNP, 多少Indel)。

Canu组装

创建完文件夹后,用如下命令进行组装:

canu -p ecoli-ont -d ./  genomeSize=4.8m  -nanopore /home/uu01/data/ont.fastq & 

Hifiasm组装

这个首先得解压,解压命令如下:

tar zxvf pacbio.fastq.tar.gz

解压到自己文件夹之后进行组装

hifiasm -o ecoli-hifi -t 1 ./hifi.fastq &

接着通过awk进行提取contig序列

awk '/^S/{print ">"$2;print $3}' ecoli-hifi.bp.p_ctg.gfa> ecoli-hifi.p_ctg.fa 

基于nucmer的基因组比对

比对hifiasmcanu的组装序列

nucmer --maxmatch ../hifiasm/ecoli-hifi.p_ctg.fa ../canu-ont/ecoli-ont.contigs.fasta &

过滤比对结果

delta-filter -i 90 -l200 -r -q out.delta >out.filter

转换为可读性强的tab键分隔的文件, -H去除标题行

show-coords -H -T -r out.filter > out.flt.tab

通过鉴定SNP和Indel,识别两种组装结果的差异

show-snps -r -T -x 5 out.filter >out.snps

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-VVdKGlKn-1670180460088)(D:/typora%E5%9B%BE%E7%89%87/image-20221205002123986.png)]

图2 比对结果

统计SNP和Indel的数量

SNP:Single-nucleotide polymorphism,单核苷酸多态性在此数据文件中表现为第二列和第三列都是字母且不一样。

Indel:Insertion-deletion,插入缺失,表现为第二列或第三列是.

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-aTqpuo2M-1670180460088)(D:/typora%E5%9B%BE%E7%89%87/image-20221205002348301.png)]

图3 out.snps

统计SNPIndel的代码如下:

library(data.table)
data<-fread("D:/00大三上/生信原理/Data/out.snps")
Indel <- sum(data[,2]=='.' |  data[,3]=='.')
# indel的个数
Indel
# snp的个数
dim(data)[1] - Indel
SNPIndel
10104700

讨论

调整参数

尝试调整参数,结果有什么改变?为什么会出现这种变化?

我调整了SOAPdenovo-K参数,也即是调整了组装kmer的大小,我调整为该参数的最大值和最小值。

SOAPdenovo-63mer all -s /home/uu01/data/ecoli.cfg -K 63 -R -o ecoli-soap &
SOAPdenovo-63mer all -s /home/uu01/data/ecoli.cfg -K 13 -R -o ecoli-soap &

K为63时,运行时间为1min,k=13时运行时间为6min。

我们在选取K-mer的大小时,应该能够使得每一个K-mer都唯一的比对到基因在上。除了试,还可以用一些工具帮助我们觉得K-mer大小,比如KmerGenie等软件。

KmerGenie (psu.edu)

K=13时,config最长为114,连接config没有得到scaffold
在这里插入图片描述
在这里插入图片描述

图4、5 K=13时的scafflod和config

K=63时评估指标如下

评估指标ContigScaffold
N502421224796
N9067045410
L506298
L90213527
总碱基数50815645006149
总序列数目352963
最长序列的长度17032607906

在组装时,由于机器读长的限制,直接采用overlap进行组装的算法效果并不好,为了提升组装效果,基于kmer的算法流行了起来。

kmer 是一段固定长度的序列,这个长度是自己定义的。

kmer大小越大,就越有可能避免图中相似区域(重复、错误等)之间的歧义。如果kmer在基因组中多次存在,就会产生歧义。因此理论上较大的kmer大小会增加N50。然而,大的kmer尺寸对测序错误、杂合性和覆盖率更敏感。

基因组组装的连续性

哪个基因组组装的连续性最好?为什么它的连续性最好?

我们总共用了三种方法进行组装:SOAPdenovoCanuhifiasm

我们可以用N50/N90L50/L90评估组装连续性,如果N50/N90越大,L50/L90越小,则组装连续性越高,预示着组装质量越好。

组装方法N50/N90L50/L90
SOAPdenovo K=514.000.16
SOAPdenovo K=634.950.09

从连续性上来看,SOAPdenovo 方法中K取51更好。

canu和hifiasm最终的结果都是一条序列,其长度分别为46333714662761

从原理上来看hifiasm连续型最好,Hifiasm使用了graph-binning的策略对此进行了改进。它不预先划分reads,而是在string-graph中对reads进行标记。因此在一个较长的bubble中,即使只有一小部分reads被正确标记,hifiasm也可以正确地将其定相。通过这种方式,可以避免因为reads划分错误而引入的错误位点和组装断裂,从而获得更完整和更准确的单倍体组装结果。

单倍体组装工具Hifiasm简介及基本运行命令(一) - 哔哩哔哩 (bilibili.com)

组装差异

同样的基因组,为什么Canu的组装序列和hifiasm的组装序列会有差异?

二者原理不同。

hifiasm的分析流程如下,主要分为3个阶段

第一阶段:通过所有序列的相互比对,对前在测序错误进行纠正。如果一个位置只存在两种碱基类型,且每个碱基类型至少有3条read支持,那么这个位置会被当作杂合位点,否则,视作测序错误,将被纠正。

第二阶段:根据序列之间的重叠关系,构建分型的字符串图(phased string graph)。其中调整朝向的序列作为顶点(vertex),一致重叠作为边(edge)。字符串图中的气泡(bubble)则是杂合位点。

第三阶段:如果没有额外的信息,hifiasm会随机选择气泡的一边构建primary assembly,另一边则是alternate assembly. 该策略和HiCanu,Falcon-Unzip一样。对于杂合基因组而言,由于存在一个以上的纯合haplotype,因此primary assembly可能还会包含haplotigs。HiCanu依赖于第三方的purge_dups, 而hifiasm内部实现了purge_dups算法的变种,简化了流程。如果有额外的信息,那么hifiasm就可以正确的对haplotype进行分型。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-c2pAcQye-1670180460089)(D:/typora%E5%9B%BE%E7%89%87/image-7b5344ff6ef24f789babd6da65f27cf4.png)]

图5 Hifiasm流程

Canu的流程分为三个步骤,前两步是原始输入的纠错,最后一步是基于纠错后的reads来构建contigs。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-1nI7m226-1670180460089)(D:/typora%E5%9B%BE%E7%89%87/image-20191206100241398-c5b8e2076fcb4120a969dea2ad203eda.png)]

图6 Canu步骤三流程图

hifiasm对HiFi PacBio进行组装(xuzhougeng.top)

Canu的graph和contig生成过程(xuzhougeng.top)

提高空间

SOAPdenovo、Canu或者hifiasm组装序列的质量(连续性、准确率、完整性),是否有进一步提高的空间?有什么办法可以提高?

由于现有的Hi-C或Strand-seq分箱算法从一个折叠装配开始,它们可能不能很好地工作在初始装配中表示的高度杂合区域。而且对多倍体植物仍然具有挑战性。

一种可能的解决办法是将Hi-C或Strand-seq数据映射到hifiasm组装图上,用图的拓扑结构将单元格分组并排序到染色体长的支架上,然后沿着支架将杂合子事件分阶段。

Cheng, Haoyu et al. “Haplotype-resolved de novo assembly using phased assembly graphs with hifiasm.” Nature methods vol. 18,2 (2021): 170-175. doi:10.1038/s41592-020-01056-5


http://www.ppmy.cn/news/865297.html

相关文章

tensorflow实现图像增强:Flip and Rotation

#tf.image.rot90()与tf.image.flip_up_down()增加原始数据到8倍 #python 3.6 import matplotlib.pyplot as plt import tensorflow as tfimage_raw_data tf.io.gfile.GFile(r./data/IMG/TaylorSwift.jpg,rb).read()with tf.compat.v1.Session() as sess:img_data tf.image.…

Qt| There‘s no Qtversion assigned to project... 解决方法

问题&#xff1a; 原因&#xff1a;相同工程在不同电脑下qt配置不一致导致&#xff0c;该项目qt setting设置有误。 解决方法&#xff1a;右键项目打开属性 找到Qt Project Settings->Qt Installation&#xff0c;切换到当前电脑所使用的qt版本即可。

7月14日每日两题

第一题:找到最大岛 小哼通过秘密方法得到一张不完整的钓鱼岛航拍地图。钓鱼岛由一个主岛和一些附属岛屿组成,小哼决定去钓鱼岛探险。下面这个10*10的二维矩阵就是钓鱼岛的航拍地图。图中数字表示海拔,0表示海洋,1~9都表示陆地。现在需要计算出最大岛的面积(即有多少个格子…

潮人专属好物!HCK哈士奇x可口可乐联名限量款小冰吧

工业革命之后&#xff0c;人们就进入了大工厂时代&#xff0c;越来越多的东西走上了流水线&#xff0c;千篇一律、毫无新意的东西正在抹杀人们的想象力和审美体验&#xff0c;于是人们对这样的东西越来越拒绝。冰箱也如此&#xff0c;人们看遍了方方正正的白色抑或灰色的“大铁…

虚拟人春节搞事情!先在央视《对话》,又跟李玉刚组团除夕出道

金磊 发自 凹非寺量子位 | 公众号 QbitAI 一个女孩登上了央视《对话》栏目&#xff0c;仅是浅唱了一首歌&#xff0c;便让全场惊叹连连。 讲真&#xff0c;这种reaction还真没有一点夸张。 话不多说&#xff0c;先来感受下这个feel&#xff1a; 或许你会问了&#xff0c;人美歌…

itchat微信助手,kaggle 电影数据集分析,基于内容的电影推荐

项目的github地址&#xff1a;https://github.com/qihe777/itchatApplication 1.项目运行 配置完成环境后&#xff1a;&#xff08;关于调用的服务接口所需要的秘钥我都没有修改&#xff0c;但为了你程序能一直运行最好自己申请相应的服务&#xff09; 运行项目&#xff1a;…

中国首个虚拟学生入学清华大学!双商在线、颜值出众,你想跟她做同学吗?

点击上方“视学算法”&#xff0c;选择加"星标"或“置顶” 重磅干货&#xff0c;第一时间送达 来源&#xff1a;学校共青团、中国新闻网、新浪财经、科技日报、江苏共青团等 6月1日&#xff0c;清华大学迎来一位特殊的学生&#xff0c;名为华智冰&#xff0c;她是一位…

小冰公司CEO李笛:AI不会江郎才尽,创造力只会持续向上攀升丨MEET2022

编辑部 整理自 MEET 2022量子位 报道 | 公众号 QbitAI 从GPT-3到DALLE&#xff0c;AI在创作这件事上的进步速度&#xff0c;远超我们想象。 但有不少观点认为&#xff0c;相比围棋、电子游戏等知识领域&#xff0c;AI在创作上不仅比不过人类&#xff0c;甚至不具备真正的创造力…