每日小结

  1. 2021-03-31
    完成
    看了点纠错知识
    服务器跑商三代数据 align,sort,merge
    明日计划
    三代数据比对结果(3970000-39870000)
    二代比对结果到(3970000-39870000)上深度信息有点奇怪等二代数据比对到chr13结果出来后再看看,
    看看三代自纠错和二代数据纠三代

  2. 2021-04-02

    • 三代数据比对到chr13中完成,结果也是比较怪,有部分长度的序列深度特别高,结果在2_chr13_pacbio_align.sort.bam中
    • 二代的在服务器上出了点问题
    • 放假前把二代数据比对到全基因组上,看看效果
    • 之后应该看纠错的知识了
  3. 2021-04-05

    • 将实验室电脑装成ubuntu系统,并配置相应的实验环境。
    • 看了看比对到全基因组上的结果,发现2.JY中是两个拷贝,不是给的4个拷贝,对于不同的拷贝,其只是有偶尔一个碱基不同
  4. 2021-04-06

    • 得到三代数据的比对结果后,发现拷贝位置的深度确实是其他位置的四倍,重新查看二代数据的比对信息时,发现比对的二代数据是选择错误的H25,重新比对二代数据
  5. 2021-04-07

    • 获得二代数据的比对结果,使用IGV软件进行查看,发现,二代数据比对深度不平稳,但也明显可以看出是其他正常区域的四倍,三代数据深度平稳,是其他位置的四倍,其中四个拷贝,不同的地方仅仅是一个碱基,考虑简单的排列组合进行组装是否可行,以及此项目是否有继续进行的必要。
    • 部分拷贝区域内的二代数据深度低的部分,三代数据中的错误率也较高,但并不绝对,这一点不仅仅在多个拷贝的区域这样,具体原因仍需要考虑(组装的参考基因组的问题?)
    • 几个拷贝之间基本上相同,主要应该解决的问题是否是几个拷贝间的连接问题。以及不同拷贝中不同部分的顺序
    • 明天看看两个拷贝接触位置的覆盖深度,从而可以确定,拷贝之间是否是相连的。以及确定三代数据的长度范围。
  6. 2021-04-08

    • 三代数据的长度从 几k 到 一百k 左右都存在。查看了拷贝接触位置的覆盖深度,没有什么结论,疑问:最后一个拷贝是否不是和之前的拷贝的长度相同,也就是最后一个拷贝的中间位置连接之后的部分
  7. 2021-04-09

    • 开会 仔细注意两边连接处的详细信息,
    • 将74k拷贝后接的于它开始相同的到nnnnnnn之间的碱基去掉,进行全基因组的比对,查看比对结果,
    • 多个拷贝连接失败的原因,假设拷贝是-AXAXAXAX-,AX是一个拷贝,第一个X后边接A的信号远强于接-的信号,但是接A就会形成的一个环,然后无限循环无法停住,所以原有的所有组装软件都拼接不起来。其次,最后一个X也并非和前面几个X的长度相同。
    • 下一步,查看删除后的比对结果,然后对三代数据进行纠错,提取出所有于三代数据有相关的reads,设计算法进行组装
  8. 2021-04-12

    • 早上查看比对结果时,发现使用minimap2比对三代数据时意外终止,猜测原因可能是自己使用多线程导致内存不够用,或者代码对多线程有bug,原来使用14个线程跑,使用8个线程时恢复正常,继续跑三代数据比对
    • 下午去健身了~~~
    • 晚上回到实验室查看比对结果,发现有一条约3k的read比对到了拷贝,以及拷贝后的序列,中间有约1k的deletion,将这条read添加到参考基因组中重新比对,虽然只有一条read偶然性较大,若比对结果理想,仍能表明这个课题基本上能有个好的结果,这可以算是一个重大进展了。
    • 第一次体会到了科研的乐趣,看到这条read时,心中有股难以言喻的喜悦之情,跑出去坐了10分钟才平息下来,哈哈。
  9. 2021-04-13

    • 查看了比对结果,coverage只有二十几,并不是很完美,但感觉也够用了,
    • 重新修改了一个参考基因组,重新比对看看结果会不会好一点
    • 修正三代数据
    • 开始看sga源码,论文
  10. 2021-04-14

    • coverage还是只有20多,而且二代数据的coverage很多位置是0,猜测原因开始是添加的三代read中错误率较高,可以等之后纠错后再找到read添加测试
    • 纠错,开始时发现软件经常意外终止,开始还以为是软件有bug,后面才知道是自己内存不够大,而且swap分区也比较小,把swap分区调大(64G)后就可以了,之前minimap2的意外终止也可能是这个原因
  11. 2021-04-15

    • ubuntu出问题了,进不去桌面,命令行正常,重装了一遍显卡的驱动后就好了,自带的显卡驱动也太辣鸡了呀
    • 查看纠错后的read,输入的read数量少了很多。。。
    • 手动将纠错之后的read添加到参考基因组上进行比对,
  12. 2021-04-16

    • 手动将纠错之后的read添加到参考基因组上进行比对后,在三代数据中,添加的和拷贝的连接处有一段的coverage比较高,有500左右,暂不清楚是什么原因
    • 将原始的参考数据中的断开的区域中的nnnn去掉后,在进行比对,查看是否能有其他的数据能够连接上这两段
    • 进行了三代数据的纠错,不知道两天多能不能跑完
    • 对两个拷贝的进行的比对
    • 继续看组装算法的论文

资料

生物信息学博客

  1. HOMOLOG.US - FRONTIER IN BIOINFORMATICS
  2. DE BRUIJN GRAPHS - I
  3. DE BRUIJN GRAPHS - II
  4. DE BRUIJN GRAPHS - III
  5. STRING GRAPH ASSEMBLER
  6. STRING GRAPH OF A GENOME

论文

  1. Long-read error correction: a survey and qualitative comparison
  2. The fragment assembly string graph
  3. Efficient construction of an assembly string graph using the FM-index
  4. Efficient de novo assembly of large genomesusing compressed data structures
  5. Fast and accurate short read alignment with Burrows–Wheeler transform

知识点

  1. mapping:只提供read比对的参考序列上的位置信息;alignment:不仅提供位置信息,还有他们的实际位点和参考位点之间的匹配、不匹配、插入或删除的关系

Fast and accurate short read alignment with Burrows–Wheelertransform

prefix trie

the prefix trie of googol
prefix trie
字符串X的前缀树中,每条边代表一个字符,每条从叶子节点到根节点的路径上的边上字符形成的字符串都是X的前缀,每个节点到根节点中的路径上的字符形成的字符串都是X的字串

Burrows–Wheeler transform

X表示字符串,X[i]表示字符串中的第i个字符,X[n-1]= $
suffix array

  • 上图左边的数组是suffix array,用 XiX_i 表示,X0X_0=‘googol$\$’,X1X_1=‘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]
    对于精确的匹配问题,只有一个解;对于不精确的问题,可能会有多个解
  • a:字符表中的一个字符
  • C(a):在X[0,n-2]中小于a(字典排序在a前面)的字符的数量
  • O(a,i):在B[0:i]中a出现的次数
    W是X的一个字串,R(aW)=C(a)+O(a,R(W)1)+1\underline R(aW)=C(a)+O(a,\underline R(W)-1)+1 , R(aW)=C(a)+O(a,R(W)\overline R(aW)=C(a)+O(a,\overline R(W) ,当且仅当R(aW)R(aw)\underline R(aW) \le \overline R(aw) 时,aW是X的字串

The fragment assembly string graph

Building string graph