CC BY 4.0 （除特别声明或转载文章外）
如果这篇博客帮助到你，可以请我喝一杯咖啡~
A02号二倍体组装实战流程：
Step1：hifiasm初步组装contig
辅助使用hic数据进行组装
name = A02 //设置编号
# 创建文件夹
mkdir ${name}_hic&&cd ${name}_hic
# 链接一下原始的read数据，方便处理
ln -s /data/jiangheling/03_Project/00.RC/01.data/01.hifi/${name}.fastq.gz .
ln -s /data/jiangheling/03_Project/00.RC/01.data/02.hic/${name}/${name}_1.fq.gz .
ln -s /data/jiangheling/03_Project/00.RC/01.data/02.hic/${name}/${name}_2.fq.gz .
hic1=${name}_1.fq.gz #使用临时别名，链接上hic_1、hic_2数据
hic2=${name}_2.fq.gz
# 运行hifiasm
hifiasm -t 20 -o ${name}.hic --h1 ${hic2} --h2 ${hic1} ./${name}.fastq.gz
不使用hic数据辅助
name=DR21 #原始数据文件号
mkdir ${name}&&cd ${name}
# 创建链接
ln -s /data/jiangheling/03_Project/00.RC/01.data/01.hifi/${name}.fastq.gz .
# 直接导入数据运行hifiasm
hifiasm -t 80 -o ${name}.asm ./${name}.fastq.gz
Step2:处理hifiasm的输出文本（以有hic数据辅助为代表）
name=A02 #临时别名
version=asm.hic #临时别名，主要为hifiasm的输出结果为${name}.hic....
cd ${name} #切换到neme文件下
#使用grep过滤提取step1结果文件下的带"S"字符串的信息，然后使用awk使用print命令提取相应数据
grep "S" ${name}.${version}.hap2.p_ctg.gfa | awk '{print ">"$2"\n"$3}' > ${name}.${version}.hap2.p_ctg.gfa.fa
grep "S" ${name}.${version}.hap1.p_ctg.gfa | awk '{print ">"$2"\n"$3}' > ${name}.${version}.hap1.p_ctg.gfa.fa
grep "S" ${name}.${version}.r_utg.gfa | awk '{print ">"$2"\n"$3}' > ${name}.${version}.r_utg.gfa.fa
grep "S" ${name}.${version}.p_ctg.gfa | awk '{print ">"$2"\n"$3}' > ${name}.${version}.p_ctg.gfa.fa
grep "S" ${name}.${version}.p_utg.gfa | awk '{print ">"$2"\n"$3}' > ${name}.${version}.p_utg.gfa.fa
#若不使用hic数据，hifiasm的输出为bp....
PS:
awk：这是一个强大的文本处理工具，常用于模式扫描和处理。
'{' ... '}'：大括号内包含的是awk要执行的命令。
print：这是一个输出命令，用于打印内容。
">"$2：这是要打印的内容之一。$2代表输入文本的每一行中的第二个字段（字段之间通常由空格或制表符分隔）。字段前面加上">"，表示在第二个字段的内容前加上大于号>。
"\n"：这是一个换行符，用于在打印完第一个字段后换行。
$3：这是输入文本每一行中的第三个字段。
综上所述，这条awk命令的作用是：
读取输入文本的每一行。
对于每一行，打印一个大于号>，后跟该行的第二个字段。
然后打印一个换行符。
最后打印该行的第三个字段。
Step3：筛选、BLAST比对以去除线粒体序列后，评估剩余序列的基因组组装质量。
step3.cpmt_clean(注意分别对hap1和hap2操作)
第一步：
name=A02#文件别名
hap=hap1
cd ${name}_hic
mkdir ${hap}.clean && cd ${hap}.clean
#链接上一步处理过的文件
ln -s ../${name}.asm1.hic.${hap}.p_ctg.gfa.fa ./${name}.${hap}.fa
#运行seqkit软件，对文本序列进行提取信息，整个命令的作用是从一个FASTA文件中提取序列的ID和长度，筛选出长度小于或等于1000000的序列，并将这些序列的ID和长度按长度排序后保存到一个新文件中。
seqkit fx2tab -n -g -l ${name}.${hap}.fa | sort -n -k 2 | awk '{if($2 <= 1000000) print $1"\t"$2}' > ${name}.${hap}.1M.len
ps:
seqkit fx2tab -n -g -l ${name}.${hap}.fa
seqkit：这是一个用于处理FASTA/Q文件的命令行工具。
fx2tab：这是seqkit的一个子命令，用于将FASTA文件转换为表格格式。
-n：这个选项表示输出序列的名字（FASTA文件中的">"后面的部分）。
-g：这个选项表示输出序列的完整ID（包括名字和后跟的任何文本，直到空格）。
-l：这个选项表示输出序列的长度。
${name}.${hap}.fa：这是一个变量，表示输入的FASTA文件名，其中${name}和${hap}是之前定义的环境变量。
这部分命令的作用是将FASTA文件转换为表格格式，每一行包含序列的ID、名称和长度。
| sort -n -k 2：
|：管道符号，用于将前一个命令的输出作为下一个命令的输入。
sort：这是对输入数据进行排序的命令。
-n：这个选项表示按数值顺序进行排序。
-k 2：这个选项表示根据每行的第二个字段（即序列长度）进行排序。
这部分命令的作用是将seqkit输出的表格按序列长度进行数值排序。
awk '{if($2 <= 1000000) print $1"\t"$2}'：
awk：这是文本处理工具，用于对输入的文本进行模式扫描和处理。
{if($2 <= 1000000) print $1"\t"$2}'：这部分是awk的命令，它检查每行的第二个字段（即序列长度）是否小于或等于1000000。
$1：第一个字段，即序列的ID或名称。
$2：第二个字段，即序列的长度。
"\t"：制表符，用于在打印的字段之间分隔。
如果条件满足（序列长度小于或等于1000000），则打印序列的ID和长度。
这部分命令的作用是筛选出长度小于或等于1000000的序列，并打印它们的ID和长度。
> ${name}.${hap}.1M.len：
>：这是重定向符号，用于将输出保存到文件中。
${name}.${hap}.1M.len：这是输出文件的名称，其中包含了之前定义的环境变量${name}和${hap}，以及表示筛选条件的后缀.1M.len。
这部分命令的作用是将awk命令的输出保存到一个文件中。
第二步：
#输入上一步筛选出来的文件，使用awk提取打印每一行的第一个字段（序列ID）并且输出到新文件
awk '{print $1}' ${name}.${hap}.1M.len > ${name}.${hap}.1M.len.name.list
#使用序列名称列表从原始FASTA文件中筛选出对应的序列，并保存到新的FASTA文件中
seqkit grep -f ${name}.${hap}.1M.len.name.list ${name}.${hap}.fa > ${name}.${hap}.1m.contig.fa
#将筛选出的序列与指定的数据库进行BLASTN比对，并将结果保存到指定的输出文件中
blastn -query ${name}.${hap}.1m.contig.fa -db /data/jiangheling/03_Project/00.RC/01.data/00.ref/OB.cp.mt.fa -num_threads 30 -max_target_seqs 5 -outfmt "7 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qcovs" -out ${name}.${hap}.1m.contig.blast.out
ps：
blastn：使用BLASTN进行核酸序列比对。
-query ${name}.${hap}.1m.contig.fa：指定查询序列的FASTA文件。
-db /data/jiangheling/03_Project/00.RC/01.data/00.ref/OB.cp.mt.fa：指定比对数据库。
-num_threads 30：指定使用的线程数。
-max_target_seqs 5：每个查询序列最多返回5个比对结果。
-outfmt "7 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qcovs"：指定输出格式。
-out ${name}.${hap}.1m.contig.blast.out：指定输出文件。
第三步：
#从BLAST输出中筛选出满足特定比对标准的序列名称，并保存到一个列表文件中
grep -v "#" ${name}.${hap}.1m.contig.blast.out | awk '{if($3 >=95 && $13 >=90)print $1}' | uniq > get.95.seq.name.final.list
#是提取原始FASTA文件中所有序列的名称，并保存到一个列表文件中。
grep ">" ${name}.${hap}.fa | sed 's/>//g' > ${name}.${hap}.fa.name.list
ps：
grep -v "#" ${name}.${hap}.1m.contig.blast.out | awk '{if($3 >=95 && $13 >=90)print $1}' | uniq > get.95.seq.name.final.list：
grep -v "#"：排除包含"#"的行，通常是BLAST输出的注释行。
awk '{if($3 >=95 && $13 >=90)print $1}'：使用awk筛选出比对身份（pident）大于等于95%且序列覆盖度（qcovs）大于等于90%的序列名称。
uniq：去除重复的行。
> get.95.seq.name.final.list：将筛选出的序列名称保存到一个列表文件中。
这条命令的作用是从BLAST输出中筛选出满足特定比对标准的序列名称，并保存到一个列表文件中。
grep ">" ${name}.${hap}.fa | sed 's/>//g' > ${name}.${hap}.fa.name.list：
grep ">"：从FASTA文件中提取序列的名称（不包括">"符号）。
sed 's/>//g'：使用sed命令删除">"符号。
> ${name}.${hap}.fa.name.list：将序列名称保存到一个列表文件中。
这条命令的作用是提取原始FASTA文件中所有序列的名称，并保存到一个列表文件中。
第四步：
#找出原始FASTA文件中的序列名称列表与通过BLAST筛选出的序列名称列表的交集，并将这些交集的序列名称保存到一个新文件中。
comm ${name}.${hap}.fa.name.list get.95.seq.name.final.list | awk -F '\t' '{print $1}' > get.${name}.${hap}.fa.final.name.list
#这条命令的作用是使用之前得到的交集序列名称列表从原始FASTA文件中筛选出对应的序列，并将这些序列保存到一个新的FASTA文件中，这个文件将不包含线粒体序列（假设.rm.cp.mt.fa表示移除了线粒体序列）。
# comm first row A-B # second row B-A
seqkit grep -f get.${name}.${hap}.fa.final.name.list ${name}.${hap}.fa > ${name}.${hap}.rm.cp.mt.fa
#这条命令的作用是使用QUAST工具评估移除了线粒体序列后的基因组组装质量，并将结果输出到指定的目录和文件中。
python /data/jiangheling/00_Biosoft/01_bin/quast.py ${name}.${hap}.rm.cp.mt.fa -o ${hap}.genome.report
ps：
comm ${name}.${hap}.fa.name.list get.95.seq.name.final.list | awk -F '\t' '{print $1}' > get.${name}.${hap}.fa.final.name.list：
comm：这是一个比较两个已排序文件的命令，它将显示两个文件的交集、只在第一个文件中的行和只在第二个文件中的行。默认情况下，输出分为三列，第一列是只在第一个文件中的行，第二列是只在第二个文件中的行，第三列是在两个文件中都有的行（交集）。
${name}.${hap}.fa.name.list：包含原始FASTA文件中所有序列名称的列表。
get.95.seq.name.final.list：包含通过BLAST筛选后满足条件的序列名称的列表。
awk -F '\t' '{print $1}'：使用awk打印第一列，即两个列表文件的交集部分。这里假设comm命令的输出使用制表符作为字段分隔符。
> get.${name}.${hap}.fa.final.name.list：将交集部分的序列名称保存到一个新列表文件中。
这条命令的作用是找出原始FASTA文件中的序列名称列表与通过BLAST筛选出的序列名称列表的交集，并将这些交集的序列名称保存到一个新文件中。
seqkit grep -f get.${name}.${hap}.fa.final.name.list ${name}.${hap}.fa > ${name}.${hap}.rm.cp.mt.fa：
seqkit grep -f：seqkit的一个子命令，用于根据提供的序列名称列表从FASTA文件中筛选序列。
get.${name}.${hap}.fa.final.name.list：包含要筛选的序列名称的列表文件。
${name}.${hap}.fa：原始的FASTA文件。
> ${name}.${hap}.rm.cp.mt.fa：将筛选出的序列保存到一个新的FASTA文件中。
这条命令的作用是使用之前得到的交集序列名称列表从原始FASTA文件中筛选出对应的序列，并将这些序列保存到一个新的FASTA文件中，这个文件将不包含线粒体序列（假设.rm.cp.mt.fa表示移除了线粒体序列）。
python /data/jiangheling/00_Biosoft/01_bin/quast.py ${name}.${hap}.rm.cp.mt.fa -o ${hap}.genome.report：
python：调用Python解释器。
/data/jiangheling/00_Biosoft/01_bin/quast.py：QUAST工具的Python脚本路径，QUAST是一个用于评估组装基因组质量的工具。
${name}.${hap}.rm.cp.mt.fa：移除了线粒体序列的FASTA文件。
-o ${hap}.genome.report：指定输出结果的目录和前缀。
这条命令的作用是使用QUAST工具评估移除了线粒体序列后的基因组组装质量，并将结果输出到指定的目录和文件中。
总结来说，这些命令组合在一起，用于从原始FASTA文件中筛选出特定长度的序列，进行BLAST比对以去除线粒体序列，最后评估剩余序列的基因组组装质量。
Step4：
第一步：
name=A02
seq=hap1
cd ${name}_hic/${seq}.clean
mkdir ${seq}.clean.3ddna&&cd ${seq}.clean.3ddna
mkdir reference&&cd reference
ln -s ../../${name}.${seq}.rm.cp.mt.fa .
cd ..
mkdir fastq&&cd fastq
ln -s /data/jiangheling/03_Project/00.RC/01.data/02.hic/${name}/${name}_1.fq.gz read_R1.fastq.gz
ln -s /data/jiangheling/03_Project/00.RC/01.data/02.hic/${name}/${name}_2.fq.gz read_R2.fastq.gz
cd ..
genome=./reference/${name}.${seq}.rm.cp.mt.fa
#为指定的参考基因组${genome}文件创建bwa比对所需的索引 。
bwa index ${genome}
#生成HindIII限制性内切酶在参考基因组上的切割位点信息，并将结果保存到文件中。
python /data/jiangheling/00_Biosoft/02_3d-dna/juicer-1.6/misc/generate_site_positions.py HindIII ${seq}.genome ${genome}
#是从HindIII切割位点信息文件中提取染色体名称和大小，并保存到一个新文件中。
awk 'BEGIN{OFS="\t"}{print $1, $NF}' ${seq}.genome_HindIII.txt > ${seq}.genome.chrom.sizes
#运行Juicer工具进行Hi-C数据分析，使用指定的参考基因组、HindIII切割位点和染色体大小信息。
/data/jiangheling/00_Biosoft/02_3d-dna/juicer-1.6/scripts/juicer.sh -g ${name}_hap1 -D /data/jiangheling/00_Biosoft/02_3d-dna/juicer-1.6/ -y ./${seq}.genome_HindIII.txt -z ${genome} -p ./${seq}.genome.chrom.sizes -t 20 -s HindIII
#运行3D-DNA工具的组装管道，使用Juicer工具的输出结果来组装和分析Hi-C数据，并将日志输出到文件中，同时命令在后台执行。
/data/jiangheling/00_Biosoft/02_3d-dna/3d-dna/run-asm-pipeline.sh -r 2 ${genome} aligned/merged_nodups.txt &> 3d.log &
第二步：
seq=hap2
cd ../../${seq}.clean
mkdir ${seq}.clean.3ddna&&cd ${seq}.clean.3ddna
mkdir reference&&cd reference
ln -s ../../${name}.${seq}.rm.cp.mt.fa .
cd ..
mkdir fastq&&cd fastq
ln -s /data/jiangheling/03_Project/00.RC/01.data/02.hic/${name}/${name}_1.fq.gz read_R1.fastq.gz
ln -s /data/jiangheling/03_Project/00.RC/01.data/02.hic/${name}/${name}_2.fq.gz read_R2.fastq.gz
cd ..
genome=./reference/${name}.${seq}.rm.cp.mt.fa
cd ..
echo 'sadq'
echo ;';'
bwa index ${genome}
python /data/jiangheling/00_Biosoft/02_3d-dna/juicer-1.6/misc/generate_site_positions.py HindIII ${seq}.genome ${genome}
awk 'BEGIN{OFS="\t"}{print $1, $NF}' ${seq}.genome_HindIII.txt > ${seq}.genome.chrom.sizes
/data/jiangheling/00_Biosoft/02_3d-dna/juicer-1.6/scripts/juicer.sh -g ${name}_hap2 -D /data/jiangheling/00_Biosoft/02_3d-dna/juicer-1.6/ -y ./${seq}.genome_HindIII.txt -z ${genome} -p ./${seq}.genome.chrom.sizes -t 20 -s HindIII
/data/jiangheling/00_Biosoft/02_3d-dna/3d-dna/run-asm-pipeline.sh -r 2 ${genome} aligned/merged_nodups.txt &> 3d.log &
ps:
bwa index ${genome}：
bwa index：这是bwa软件的一个命令，用于为给定的参考基因组创建索引。索引是进行序列比对所必需的。
${genome}：这是一个变量，代表参考基因组的FASTA文件路径。
这条命令的作用是为指定的参考基因组文件创建bwa比对所需的索引。
python /data/jiangheling/00_Biosoft/02_3d-dna/juicer-1.6/misc/generate_site_positions.py HindIII ${seq}.genome ${genome}：
python：调用Python解释器。
/data/jiangheling/00_Biosoft/02_3d-dna/juicer-1.6/misc/generate_site_positions.py：Juicer工具集中用于生成限制性内切酶切割位点的Python脚本。
HindIII：这是使用的限制性内切酶类型。
${seq}.genome：一个变量，代表基因组序列的标识。
${genome}：参考基因组的FASTA文件路径。
这条命令的作用是生成HindIII限制性内切酶在参考基因组上的切割位点信息，并将结果保存到文件中。
awk 'BEGIN{OFS="\t"}{print $1, $NF}' ${seq}.genome_HindIII.txt > ${seq}.genome.chrom.sizes：
awk：文本处理工具，用于模式扫描和处理。
'BEGIN{OFS="\t"}{print $1, $NF}'：awk命令，设置输出字段分隔符为制表符，并打印每一行的第一个字段和最后一个字段。
${seq}.genome_HindIII.txt：包含HindIII切割位点信息的文件。
> ${seq}.genome.chrom.sizes：将awk的输出重定向到新文件，该文件包含染色体的大小信息。
这条命令的作用是从HindIII切割位点信息文件中提取染色体名称和大小，并保存到一个新文件中。
/data/jiangheling/00_Biosoft/02_3d-dna/juicer-1.6/scripts/juicer.sh -g ${name}_hap1 -D /data/jiangheling/00_Biosoft/02_3d-dna/juicer-1.6/ -y ./${seq}.genome_HindIII.txt -z ${genome} -p ./${seq}.genome.chrom.sizes -t 20 -s HindIII：
/data/jiangheling/00_Biosoft/02_3d-dna/juicer-1.6/scripts/juicer.sh：Juicer工具集中用于执行Hi-C数据分析的脚本。
-g ${name}_hap1：设置基因组标识。
-D /data/jiangheling/00_Biosoft/02_3d-dna/juicer-1.6/：设置Juicer工具的路径。
-y ./${seq}.genome_HindIII.txt：指定包含HindIII切割位点信息的文件路径。
-z ${genome}：指定参考基因组文件的路径。
-p ./${seq}.genome.chrom.sizes：指定染色体大小文件的路径。
-t 20：设置使用的线程数。
-s HindIII：指定使用的限制性内切酶。
这条命令的作用是运行Juicer工具进行Hi-C数据分析，使用指定的参考基因组、HindIII切割位点和染色体大小信息。
/data/jiangheling/00_Biosoft/02_3d-dna/3d-dna/run-asm-pipeline.sh -r 2 ${genome} aligned/merged_nodups.txt &> 3d.log &：
/data/jiangheling/00_Biosoft/02_3d-dna/3d-dna/run-asm-pipeline.sh：3D-DNA工具集中用于组装和分析Hi-C数据的脚本。
-r 2：设置迭代次数。
${genome}：参考基因组的FASTA文件路径。
aligned/merged_nodups.txt：Juicer工具输出的合并和去重后的比对结果文件。
&> 3d.log &：将标准输出和标准错误重定向到日志文件3d.log，并在后台运行该命令。
这条命令的作用是运行3D-DNA工具的组装管道，使用Juicer工具的输出结果来组装和分析Hi-C数据，并将日志输出到文件中，同时命令在后台执行。
1. hifiasm软件使用教程：
1.1 基础语法：
$ hifiasm
Usage: hifiasm [options] <in_1.fq> <in_2.fq> <...>
Options:
Input/Output:
-o STR prefix of output files [hifiasm.asm]//输出文件 [***.asm]
-i ignore saved read correction and overlaps//
-t INT number of threads [1]//使用cpu线程数
-z INT length of adapters that should be removed [0]
--version show version number
Overlap/Error correction:
-k INT k-mer length (must be <64) [51]
-w INT minimizer window size [51]
-f INT number of bits for bloom filter; 0 to disable [37]
-D FLOAT drop k-mers occurring >FLOAT*coverage times [5.0]
-N INT consider up to max(-D*coverage,-N) overlaps for each oriented read [100]
-r INT round of correction [3]
Assembly:
-a INT round of assembly cleaning [4]
-m INT pop bubbles of <INT in size in contig graphs [10000000]
-p INT pop bubbles of <INT in size in unitig graphs [0]
-n INT remove tip unitigs composed of <=INT reads [3]
-x FLOAT max overlap drop ratio [0.8]
-y FLOAT min overlap drop ratio [0.2]
-u disable post join contigs step which may improve N50
--lowQ INT
output contig regions with >=INT% inconsistency in BED format; 0 to disable [70]
--b-cov INT
break contigs at positions with <INT-fold coverage; work with '--m-rate'; 0 to disable [0]
--h-cov INT
break contigs at positions with >INT-fold coverage; work with '--m-rate'; -1 to disable [-1]
--m-rate FLOAT
break contigs at positions with <=FLOAT*coverage exact overlaps;
only work with '--b-cov' or '--h-cov'[0.75]
--primary output a primary assembly and an alternate assembly
Trio-partition:
-1 FILE hap1/paternal k-mer dump generated by "yak count" []
-2 FILE hap2/maternal k-mer dump generated by "yak count" []
-c INT lower bound of the binned k-mer's frequency [2]
-d INT upper bound of the binned k-mer's frequency [5]
-3 FILE list of hap1/paternal read names []
-4 FILE list of hap2/maternal read names []
--t-occ INT
force remove unitigs with >INT unexpected haplotype-specific reads;
ignore graph topology; [60]
Purge-dups:
-l INT purge level. 0: no purging; 1: light; 2/3: aggressive [0 for trio; 3 for unzip]
-s FLOAT similarity threshold for duplicate haplotigs [0.75 for -l1/-l2, 0.55 for -l3]
-O INT min number of overlapped reads for duplicate haplotigs [1]
--purge-cov INT
coverage upper bound of Purge-dups [auto]
--n-hap INT
number of haplotypes [2]
Hi-C-partition:
--h1 FILEs file names of Hi-C R1 [r1_1.fq,r1_2.fq,...]
--h2 FILEs file names of Hi-C R2 [r2_1.fq,r2_2.fq,...]
--seed INT RNG seed [11]
--n-weight INT
rounds of reweighting Hi-C links [3]
--n-perturb INT
rounds of perturbation [10000]
--f-perturb FLOAT
fraction to flip for perturbation [0.1]
1.2 输出结果：
1.2.1对于非三重组装：
prefix.r_utg.gfa(Haplotype-resolved raw unitig graph in GFA format)：保留了组装生成的所有单体型信息，包括体细胞突变和重复的测序错误。
prefix.p_utg.gfa(Haplotype-resolved processed unitig graph without small bubbles)：无小气泡的单倍型解析；去掉由于体细胞突变和数据背景噪音引起的small bubbles（这个并不是真正的单体型信息），对于高度杂合基因组物种优先选择这个结果。
prefix.p_ctg.gfa(Primary assembly contig graph)：对于低杂合度物种来说，优先选择该文件；对于高杂合度物种，该结果代表其中一个单倍型。
prefix.a_ctg.gfa(Alternate assembly contig graph)：组装出来的另一套单体型基因组结果。
1.2.2 对于三重组装：
prefix.r_utg.gfa(Haplotype-resolved raw unitig graph in GFA format):保存了所有的单倍型信息。
prefix.hap1.p_ctg.gfa(Phased paternal/haplotype1 contig graph):保留了阶段性父系/单倍型1组装。
prefix.hap2.p_ctg.gfa(Phased maternal/haplotype2 contig graph):保留了阶段性母系/单倍型2组装。