生信人的一天~HIFI数据+HIC数据组装基因组

沐雨@

已于 2022-08-01 21:29:09 修改

阅读量5.9k

点赞数 15

分类专栏：基因组组装文章标签： python 开发语言生物信息

于 2022-07-21 16:47:29 首次发布

本文链接：https://blog.csdn.net/qq_43337286/article/details/125913258

版权

基因组组装专栏收录该内容

4 篇文章 3 订阅

订阅专栏

HIFI加HIC数据组装基因组遇坑记@TOC

最近有一个大项目(大难题)自学基因组组装

生信入门这么久，一直都是使用别人处理好的数据，何时我才能产出自己的数据呢？？？

干-- --实验想要自己产出数据，那就得学组装

首先我拿到的数据是HIFI数据和HIC数据。简单的做了一下了解。HIFI 数据是准确率高的三代数据，不需要纠错就可以进行组装。HIC则是基于染色体构象捕获技术，利用高通量测序技术…。总之HIFI 数据可以理解为染色体片段的具体序列，而HIC数据则标记了大概空间位置上是哪一串序列。这样我对于我将要组装的路径有了一个大致的了解。首先我需要将我的HIFI 数据稍微组装以下，变成contig或者scaffold。然后根据再利用HIC 数据将我的片段挂成染色体。然后剩下的就是具体的实际操作了。

对于刚刚入行的我就随机挑选幸运软件来组装我的数据了

我所挑选的软件就是hifiasm和juicer外加3d-dna

HIFIASM记录

首先我使用HIFIASM将我的HIFI数据直接结合HIC数据组装成为超大号的contig。

  hifiasm --primary -o name -t 26 --h1 hic_r1.fq --h2 hic_r2.fq hifi.fa>HIFI_HIC.txt

这一步得到了众多文件，其中包括几个分型的基因组，但我记得我的确输入的–primary参数就是不组装分型，但是还是会出现这个问题。只不过这一点也不影响，因为我所需要的结果还是有的。
结果文件中我使用到的文件是xxx.p_ctg.gfa。到这一步我算是得到了我的超长contig文件。这个时候我需要将我的gfa格式转为fasta格式。程序猿的高光时刻到了：

# -*- encoding: utf-8 -*-
'''
@File        :gfatofasta.py
@Time        :2022/07/21 12:44:00
@Author        :charles kiko
@Version        :1.0
@Contact        :charles_kiko@163.com
@Desc        :gfatofasta
'''
import sys
from tqdm import trange
from Bio.Seq import Seq
from Bio import SeqIO
from Bio.SeqRecord import SeqRecord

IN_FILENAME = sys.argv[1]
OUT_FILENAME = sys.argv[2]
print(f"Converting GFA {IN_FILENAME} --> FASTA {OUT_FILENAME}...")
num_seqs = 0
out_put = open(OUT_FILENAME,'w')
in_put = open(IN_FILENAME,'r').readlines()
for i in trange(len(in_put)):
	line = in_put[i].strip('\n').split()
	if line[0] == 'S':
		num_seqs += 1
		simple_seq = Seq(line[2])
		simple_seq_r = SeqRecord(simple_seq, id=line[1])
		SeqIO.write(simple_seq_r, out_put, "fasta")
out_put.close()
print(f"FASTA of {num_seqs} sequences created at {OUT_FILENAME}.")

做完这一切我得到了contig的fasta文件。

JUICER记录

juicer的使用比较繁琐，论坛的反应也不是特别快，而且貌似需要谷歌邮箱才能登录，作为一个内地学生，这个渠道有点难搞哦。
首先呢我需要得到酶切位点文件。我就按照徐大佬的博客选择了DpnII。使用juicer所带的generate_site_positions.py脚本对我的contig文件进行操作。(这里我不太明白这么操作的意义是啥，选择不同的酶对组装有多大影响我也不太清楚，毕竟别人发表基因组的论文也不会说自己组装的具体操作，另一方面老板也急着要结果，我先出一个粗略版)

python generate_site_positions.py DpnII name xxx.p_ctg.fasta

这一步生成一个文件为name_DpnII.txt，后续步骤还需要一个染色体长度文件。这个我们使用awk来解决。

    awk 'BEGIN{OFS="\t"}{print $1, $NF}' name_DpnII.txt >name.chrom.sizes

到这里juicer运行所需要的文件就集齐了。接下来我们来整理目录结构。
我们生成一个文件夹为juicer
其中文件夹及文件放置情况如下：

juicer
- HIC
  - fastq
    - hic_r1.fastq
    - hic_r2.fastq
- restriction_sites
  - name.chrom.sizes
  - name_DpnII.txt
- reference
  - xxx.p_ctg.fsata
    好啦，到这里就有juicer运行的条件了，接下来就run起来。(这次的命令有点长，你要~~!)
juicer.sh -g name -d …/juicer -p ./restriction_sites/name.chrom.sizes -y ./restriction_sites/name_DpnII.txt -z ./reference/xxx.hic.p_ctg.fasta -D /apt/juicer -t 26
这里小编推荐大家都是用绝对路径，相对路径可能报错。
好啦接下来就是重度BUG区。熊猫烧香，诸BUG退避。总的来说错误就一类，缺少结果！！！！！

缺少merged_nodups.txt状态(未解决)；
在我反复运行N次之后出现了第二种错误没有inter.hic、inter_30.hic！！！；

对于问题1，我在论坛上找到一个推荐是运行juicer的下游程序。于是乎我照做了！

 java -jar juicer_tools.3.0.0.jar arrowhead inter.hic out

结果显示0 domain写入文件！

java -jar juicer_tools.3.0.0.jar hiccups inter.hic out1

结果是两个文件，文件内容看不懂，也不知道咋运用。
对于问题2，是出现在问题1的基础之上的。至今没有找到解决办法。

‘’‘我翻开论坛一查，这问题没有年代，歪歪斜斜的每页上都写着‘怎么解决’四个字。
我横竖睡不着，仔细看了半夜，才从字缝里看出字来，满本都写着两个字是‘-help’！’’’

不知道有多少生信人在此迷茫，我曾经在各种生信群里求助最后消息都如同泥牛入海。为此我重新建立一个群，进群的小伙伴备注以下自己使用什么软件。咱尽量做到基因组组装方面有问必答。也欢迎大佬来群视察，要是能开个基因组组装的讲座啥的就再好不过了！
进群请charles_kiko@163.com 备注基因组加群。