近期做hic挂载发现一个物种的染色体有问题，需要把对应两条染色体的三代测序数据挑出，比对到两条染色体上查看深度，以及挑出这部分序列重新拼接，拿到完整的染色体。
#!/bin/bash
# 设置比对基因组和HiFi数据的路径
contaminant_genome="/path2fasta"
hifi_reads="/path2fastq"
output_dir="/home/output"
line="HIFI"
# 创建输出目录
mkdir -p $output_dir/map2out
mkdir -p $output_dir/clean_hifi
# 1. 创建索引
minimap2 -d $output_dir/contaminant_genome.mmi $contaminant_genome -t 12
# 2. 将 HiFi reads 比对到基因组
minimap2 -ax map-pb -k 25 -B 5 -m 150 -A 2 -O 5,32 -E 3,2 --secondary=no $output_dir/contaminant_genome.mmi $hifi_reads -t 14 > $output_dir/map2out/${line}.contaminant.sam
# 3. 对 SAM 文件进行排序
#samtools sort -@ 14 -O bam -o $output_dir/map2out/${line}.contaminant_sorted.bam $output_dir/map2out/${line}.contaminant.sam
samtools index $output_dir/map2out/${line}.contaminant_sorted.bam
# 4. 提取比对到 scaffold_5 和 scaffold_6 的 reads
samtools view -b $output_dir/map2out/${line}.contaminant_sorted.bam scaffold_5 scaffold_6 > $output_dir/clean_hifi/${line}.scaffold_specific.bam
# 5. 将提取的 BAM 文件转换为 FASTQ 格式
samtools bam2fq $output_dir/clean_hifi/${line}.scaffold_specific.bam > $output_dir/clean_hifi/${line}_scaffold_specific_clean.fastq
