基因组组装
每日小结
-
2021-03-31
完成
看了点纠错知识
服务器跑商三代数据 align,sort,merge
明日计划
三代数据比对结果(3970000-39870000)
二代比对结果到(3970000-39870000)上深度信息有点奇怪等二代数据比对到chr13结果出来后再看看,
看看三代自纠错和二代数据纠三代 -
2021-04-02
- 三代数据比对到chr13中完成,结果也是比较怪,有部分长度的序列深度特别高,结果在2_chr13_pacbio_align.sort.bam中
- 二代的在服务器上出了点问题
- 放假前把二代数据比对到全基因组上,看看效果
- 之后应该看纠错的知识了
-
2021-04-05
- 将实验室电脑装成ubuntu系统,并配置相应的实验环境。
- 看了看比对到全基因组上的结果,发现2.JY中是两个拷贝,不是给的4个拷贝,对于不同的拷贝,其只是有偶尔一个碱基不同
-
2021-04-06
- 得到三代数据的比对结果后,发现拷贝位置的深度确实是其他位置的四倍,重新查看二代数据的比对信息时,发现比对的二代数据是选择错误的H25,重新比对二代数据
-
2021-04-07
- 获得二代数据的比对结果,使用IGV软件进行查看,发现,二代数据比对深度不平稳,但也明显可以看出是其他正常区域的四倍,三代数据深度平稳,是其他位置的四倍,其中四个拷贝,不同的地方仅仅是一个碱基,考虑简单的排列组合进行组装是否可行,以及此项目是否有继续进行的必要。
- 部分拷贝区域内的二代数据深度低的部分,三代数据中的错误率也较高,但并不绝对,这一点不仅仅在多个拷贝的区域这样,具体原因仍需要考虑(组装的参考基因组的问题?)
- 几个拷贝之间基本上相同,主要应该解决的问题是否是几个拷贝间的连接问题。以及不同拷贝中不同部分的顺序
- 明天看看两个拷贝接触位置的覆盖深度,从而可以确定,拷贝之间是否是相连的。以及确定三代数据的长度范围。
-
2021-04-08
- 三代数据的长度从 几k 到 一百k 左右都存在。查看了拷贝接触位置的覆盖深度,没有什么结论,疑问:最后一个拷贝是否不是和之前的拷贝的长度相同,也就是最后一个拷贝的中间位置连接之后的部分
-
2021-04-09
- 开会 仔细注意两边连接处的详细信息,
- 将74k拷贝后接的于它开始相同的到nnnnnnn之间的碱基去掉,进行全基因组的比对,查看比对结果,
- 多个拷贝连接失败的原因,假设拷贝是-AXAXAXAX-,AX是一个拷贝,第一个X后边接A的信号远强于接-的信号,但是接A就会形成的一个环,然后无限循环无法停住,所以原有的所有组装软件都拼接不起来。其次,最后一个X也并非和前面几个X的长度相同。
- 下一步,查看删除后的比对结果,然后对三代数据进行纠错,提取出所有于三代数据有相关的reads,设计算法进行组装
-
2021-04-12
- 早上查看比对结果时,发现使用minimap2比对三代数据时意外终止,猜测原因可能是自己使用多线程导致内存不够用,或者代码对多线程有bug,原来使用14个线程跑,使用8个线程时恢复正常,继续跑三代数据比对
- 下午去健身了~~~
- 晚上回到实验室查看比对结果,发现有一条约3k的read比对到了拷贝,以及拷贝后的序列,中间有约1k的deletion,将这条read添加到参考基因组中重新比对,虽然只有一条read偶然性较大,若比对结果理想,仍能表明这个课题基本上能有个好的结果,这可以算是一个重大进展了。
- 第一次体会到了科研的乐趣,看到这条read时,心中有股难以言喻的喜悦之情,跑出去坐了10分钟才平息下来,哈哈。
-
2021-04-13
- 查看了比对结果,coverage只有二十几,并不是很完美,但感觉也够用了,
- 重新修改了一个参考基因组,重新比对看看结果会不会好一点
- 修正三代数据
- 开始看sga源码,论文
-
2021-04-14
- coverage还是只有20多,而且二代数据的coverage很多位置是0,猜测原因开始是添加的三代read中错误率较高,可以等之后纠错后再找到read添加测试
- 纠错,开始时发现软件经常意外终止,开始还以为是软件有bug,后面才知道是自己内存不够大,而且swap分区也比较小,把swap分区调大(64G)后就可以了,之前minimap2的意外终止也可能是这个原因
-
2021-04-15
- ubuntu出问题了,进不去桌面,命令行正常,重装了一遍显卡的驱动后就好了,自带的显卡驱动也太辣鸡了呀
- 查看纠错后的read,输入的read数量少了很多。。。
- 手动将纠错之后的read添加到参考基因组上进行比对,
-
2021-04-16
- 手动将纠错之后的read添加到参考基因组上进行比对后,在三代数据中,添加的和拷贝的连接处有一段的coverage比较高,有500左右,暂不清楚是什么原因
- 将原始的参考数据中的断开的区域中的nnnn去掉后,在进行比对,查看是否能有其他的数据能够连接上这两段
- 进行了三代数据的纠错,不知道两天多能不能跑完
- 对两个拷贝的进行的比对
- 继续看组装算法的论文
资料
生物信息学博客
- HOMOLOG.US - FRONTIER IN BIOINFORMATICS
- DE BRUIJN GRAPHS - I
- DE BRUIJN GRAPHS - II
- DE BRUIJN GRAPHS - III
- STRING GRAPH ASSEMBLER
- STRING GRAPH OF A GENOME
论文
- Long-read error correction: a survey and qualitative comparison
- The fragment assembly string graph
- Efficient construction of an assembly string graph using the FM-index
- Efficient de novo assembly of large genomesusing compressed data structures
- Fast and accurate short read alignment with Burrows–Wheeler transform
知识点
- mapping:只提供read比对的参考序列上的位置信息;alignment:不仅提供位置信息,还有他们的实际位点和参考位点之间的匹配、不匹配、插入或删除的关系
Fast and accurate short read alignment with Burrows–Wheelertransform
prefix trie
the prefix trie of googol
字符串X的前缀树中,每条边代表一个字符,每条从叶子节点到根节点的路径上的边上字符形成的字符串都是X的前缀,每个节点到根节点中的路径上的字符形成的字符串都是X的字串
Burrows–Wheeler transform
X表示字符串,X[i]表示字符串中的第i个字符,X[n-1]= $
- 上图左边的数组是suffix array,用 表示,=‘googol’,=‘oogol’
- 上图右边的数组是经过字典排序后的suffix array
- S[i]: 经过字典排序后的第i个字符串在排序前(suffix array)的位置
- B[i]: 经过字典排序后的第i个字符串的最后一个字符
Burrows–Wheeler transform 转换后的结果是S和B
Suffix array interval and sequence alignment
对于一条字符串,SA interval 就是这条字符串在 经过字典排序后的suffix array 中作为子串的位置区间,SA interval 就是 prifix trie 中节点上的数字
例如:
gol
只在字典排序后的suffix array位置为1的字符串中作为前缀,它所对应的就是[1,1]o
在字典排序后的suffix array位置为4,5,6的字符串中作为前缀,它所对应的就是[4,6]
对于精确的匹配问题,只有一个解;对于不精确的问题,可能会有多个解
Exact matching: backward search
- a:字符表中的一个字符
- C(a):在X[0,n-2]中小于a(字典排序在a前面)的字符的数量
- O(a,i):在B[0:i]中a出现的次数
W是X的一个字串, , ,当且仅当 时,aW是X的字串
The fragment assembly string graph
Building string graph
All articles in this blog are licensed under CC BY-NC-SA 4.0 unless stating additionally.
Comment
GitalkValine