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 SeqRecordIN_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 += 1simple_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 备注基因组加群。

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

引用另一模板的宏_生信人值得拥有的编程模板Shell
前言 "工欲善其事必先利其器",生信工程师每天写代码.搭流程,而且要使用至少三门编程语言,没有个好集成开发环境(IDE,Integrated Development Environ ...
生存分析系列教程（一）使用生信人工具盒进行生存分析
生信人工具盒是生信人团队的开发的一款软件,非常方便.下面我将演示一下如何通过这款软件进行生存分析.为了方便大家理解,形式依然是数据结构-操作-结果解读. 1. 表达矩阵与生存信息矩阵表达矩阵依然 ...
生信c语言,生信人的R使用
接下来介绍R语言: [生信技能树]生信人应该这样学R语言 R语言在你开始R之旅前,建议你看看下面这两个 1. 介绍R语言及Rstudio 了解R,Rstudio及R包;安装的包在packages中检 ...
生信人的20个R语言习题的答案
这是生信技能树关于生信人的20个R语言习题的答案: 1 安装R包数据包: ALL, CLL, pasilla, airway 软件包:limma,DESeq2,clusterProfiler 工具包 ...
python perl 比较生信_科学网—生信人写程序1. Perl语言模板及配置 - 刘永鑫的博文...
科学网对Markdown排版支持较差,对格式不满意的用户请跳转至 CSDN 或微信阅读: 如果感觉文章对您有帮助,想继续阅读同类文章,请扫描下方二维码关注"生信宝典"公众号,每天接 ...
生信人值得拥有的编程模板-Shell
前言 "工欲善其事必先利其器",生信工程师每天写代码.搭流程,而且要使用至少三门编程语言,没有个好集成开发环境(IDE,Integrated Development Environ ...
生信人的自我修养：Linux 命令速查手册
标题:生信人的自我修养:Linux 命令速查手册目标:致力于为生信人打造一个完整的 Linux 命令速查手册作者:简佐义(jianzuoyi@qq.com) 版本:1.0 日期:2020-11-2 ...
2019linux考试题库,生信人的linux考试 2019-05-03
生信人的linux考试(20题) 1.在任意文件夹下面创建形如 1/2/3/4/5/6/7/8/9 格式的文件夹系列 mkdir 生信人的linux考试 ls cd 生信人的linux考试/ mkdi ...
生信人值得拥有的编程模板-Perl
为什么要学编程图1. 重复工作任务量与时间关系[1] 如上图,对于大量重复工作,非编程者(non-geek)工作量和时间是正相关的,就像富士康流水线上的工人,这种工作对于高智商的人是无法忍受(富士康 ...

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