Hybpiper和Paragone-nf的使用
1. 前言
1.1 Hybpiper
官方介绍：HybPiper 专为目标序列捕获而设计，其中 DNA 测序文库针对感兴趣的基因区域进行富集，尤其是针对系统发育。HybPiper 是一套 Python 脚本/模块，可包装和连接生物信息学工具，以便从高通量 DNA 测序读数中提取目标序列。
说人话版本：其实就是结合了比对软件（Bwa Blastx Diamond）+组装软件（Spades），基于二代测序数据把特定的基因组装出来。举个例子，我想获得拟南芥的matK基因，就可以通过Hybpiper把这个基因组装出来，而不再需要组装整个叶绿体基因组然后在组装注释从而提取matK基因。
1.2 Pargone-nf
官方介绍：该流程利用了 Yang 和 Smith 2014 年在此处描述和实施的同源性解析（也称为直系同源性推断）方法。这些方法后来被改编为目标捕获数据集。
说人话版本：利用Yang & Smith中的方法来剔除掉Hybpiper结果中的旁系同源基因。作者将这个方法写成了一个Pipline，可以直接以Hybpiper的结果作为输入文件，得到剔除旁系同源基因的直系同源基因。
2. 软件下载
mamba create -n hybpiper
mamba activate hybpiper
#Hybpiper
mamba install bioconda::hybpiper
# Paragone-nf
singularity pull library://chrisjackson-pellicle/collection/hybpiper-paragone:latest
git clone https://github.com/chrisjackson-pellicle/paragone-nf.git
# 修改配置文件 paragone.config
更改一下内容，
container = ‘library://chrisjackson-pellicle/collection/hybpiper-paragone:latest’
更改为，
container = ‘file://+下载的容器的sif文件的绝对路径’。
例如container = ‘file:///home/majunpeng/Software/hybpiper-nf/hybpiper-paragone_latest.sif’
3. 软件运行
3.1 运行Hybpiper
Hybpiper 分为多个模块，常用的一般就是四个模块，即，assemble（组装）+stats（统计组装情况）+recovery_heatmap（可视化基因恢复情况）+paralog_retriever（收集每个基因的的旁系同源基因并且绘制拷贝数的热图）。
首先要准备输入文件（最麻烦的一步，Hybpiper对输出文件的格式要求很严格）。
（1）二代测序数据文件（可以是全基因组测序的DNA数据，也可以是RNA-seq），文件名中只能出现一次区分双端数据使用的下划线，并且下划线之前的字符串会自动被识别为样本的名字。例如，Sample1_1.fa.gz，样本名会被识别为Sample1。
（2）目标基因（你想组装出来的基因的参考文件，如果想组装出来多个基因，请cat到一个文件，一般选择近缘物种的基因作为参考基因，可以有把多个物种的基因cat到一起，软件会自动选择一个相似性最高的作为参考基因组进行组装）。目标文件的要求同样刁钻，基因ID必须是“Taxon-Genename”，即样本名+基因名字。例如Sample-Gene1，软件会自动把Sample识别为样本名，Gene1被识别为基因名。
hybpiper assemble -r ~/sda2/Hybpiper/Fastq_dir/Chrysanthemum-dichroum-RNA_1_RNA_clean.fq.gz ~/sda2/Hybpiper/Fastq_dir/Chrysanthemum-dichroum-RNA_2_RNA_clean.fq.gz -t_dna Chr_2.fasta –bwa –cpu 100
###
-r 二代测序文件
–bwa 比对软件设置为bwa
–cpu 比对软件使用的线程数
-t_dna 目标基因文件
###
以下为全部的参数
###
必选输入
–readfiles READFILES [READFILES …]：指定输入测序文件（FASTQ格式）。若提供两个文件，视为双端测序数据（R1和R2）。例如：–readfiles sample_R1.fq sample_R2.fq。
–targetfile_dna TARGETFILE_DNA：核苷酸靶标文件（FASTA格式），用于BWA比对。序列头格式需为>TaxonID-geneName（如>Human-ITS）。
–targetfile_aa TARGETFILE_AA：氨基酸靶标文件（FASTA格式），用于DIAMOND比对，适用于跨物种分析或RNA数据。
可选输入
–unpaired UNPAIRED：添加单端读段文件，与双端数据联合使用。例如：–unpaired sample_single.fq。
读段比对选项
–bwa：使用BWA进行核苷酸比对（需–targetfile_dna），适用于高相似度靶标（如种内探针）。
–diamond：使用DIAMOND进行氨基酸比对（需–targetfile_aa），适合跨物种分析。
–diamond_sensitivity：调整DIAMOND比对敏感度，选项包括mid-sensitive到ultra-sensitive，默认中等敏感度。建议高灵敏度模式（如ultra-sensitive）提高比对率。
–evalue：BLASTx/DIAMOND比对的e值阈值，默认0.0001。降低至1e-5可减少假阳性。
–max_target_seqs：保留的比对目标序列数，默认10。
读段分发选项
–distribute_low_mem：启用低内存模式，降低内存占用但速度减慢40-50%，适合大规模数据集。
组装选项
–cov_cutoff：SPAdes覆盖度阈值，默认8。低覆盖数据可降至5，高覆盖数据可提高至10。
–single_cell_assembly：启用单细胞模式（MDA扩增数据），适用于微量DNA样本。
–kvals：手动指定SPAdes的k-mer值（默认自动选择），例如–kvals 21,33,55。
–merged：合并双端读段后组装，提升长片段重建能力。
–timeout_assemble_reads：终止耗时超过平均时间X%的基因组装任务。
Contig提取选项
–not_protein_coding：使用BLASTn代替Exonerate提取非编码基因的contig。
–thresh：保留比对的最低相似度，默认55%。Exonerate基于氨基酸，BLASTn基于核苷酸。
–extract_contigs_blast_evalue：BLASTn的e值阈值，默认10。
–paralog_min_length_percentage：contig需≥参考序列长度的百分比，默认0.75。
–depth_multiplier：主旁系同源覆盖度需≥其他旁系的倍数，默认10。
–chimeric_stitched_contig_check：检测拼接contig是否为嵌合体（默认关闭）。
–bbmap_subfilter：允许的最大错配数，默认7。
–trim_hit_sliding_window_size：滑动窗口大小，Exonerate默认5，BLASTn默认15。
–trim_hit_sliding_window_thresh：窗口相似度阈值，Exonerate默认75%，BLASTn默认65%。
通用流程选项
–start_from/–end_with：指定流程的起始和终止步骤，例如–start_from assemble_reads。
–force_overwrite：强制覆盖已有结果文件（慎用）。
–cpu：限制使用的CPU核心数（默认：总核心数-1）。
–keep_intermediate_files：保留中间文件（默认删除），用于调试。
–verbose_logging：启用详细日志（日志体积可能增大10倍）。
–hybpiper_output：指定输出目录，默认生成样本名文件夹。
###
输出结果
（1）以样本名命名的输出文件夹，例如样本名如果为SampleNale，那么最后恢复的序列位于，SampleName/GeneName/SampleName/sequences，组装的序列就在FAA（蛋白）和FNA（cds）文件夹下。
接下来运行stats模块
准备一个Sample.txt文件，储存了样本的名字，每行一个样本。例如，
Sample1
Sample2
Sample3
# 运行
hybpiper stats -t_dna ref_allgenes.fasta gene namelist.txt
###
-t 目标基因文件
###
输出结果
（1）seq_lengths.tsv ：后续绘制基因的恢复情况热图需要的文件
（2）hybpiper_stats.tsv：作用未知，感觉没什么用
运行recovery_heatma可视化基因恢复情况
hybpiper recovery_heatmap seq_lengths.tsv
可能有用的参数
###
–figure_length：输出的图片长度，如果恢复的基因数量很多，可以通过这个参数来降低高度，从而缩短每个热图块的长度，让图看起来舒服一些
–figure_height：输出的图片高度
–no_xlabels 不显示x轴标签，很有用，如果基因数果冻，缩短了长度会让所有基因的ID挤在一起
–no_ylabels 不显示y轴标签
###
输出结果
recovery_heatmap.png：基因恢复情况热图
运行paralog_retriever收集每个基因的假定的旁系同源基因
hybpiper paralog_retriever -t_dna ../Chr_6.fasta --fasta_dir_all paralogs_all --fasta_dir_no_chimeras paralogs_no_chimeras --paralogs_list_threshold_percentage 0 --heatmap_filename paralog_heatmap Sample.txt --figure_length 10 --no_xlabels
运行这一步后，会就在输出文件中得到每个基因的每个物种组装出来的每个基因，包括直系和假定的旁系，接下来使用Pagone-nf去除旁系同源基因
nextflow run ~/Software/paragone-nf/paragone-nf/paragone.nf -c ~/Software/paragone-nf/paragone-nf/paragone.config -profile standard_singularity --gene_fasta_directory paralogs_all --external_outgroups_file ../Chr_6.fasta --external_outgroups Artemisia_annua,Chrysanthemum_rhombifolium --mo --mi --rt --threads 10 --outdir paragone_result
# --external_outgroups_file 指定外类群序列文件
# -external_outgroups 指定外类群物种
最后的输出结果是一个文件夹，其中23 24 25 是后续分析要用到的数据集
