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
- fastq
- restriction_sites
- name.chrom.sizes
- name_DpnII.txt
- reference
- xxx.p_ctg.fsata
好啦,到这里就有juicer运行的条件了,接下来就run起来。(这次的命令有点长,你要~~!)
- xxx.p_ctg.fsata
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退避。总的来说错误就一类,缺少结果!!!!! - HIC
- 缺少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 备注基因组加群。