近期做hic掛載發(fā)現(xiàn)一個物種的染色體有問題垒在，需要把對應(yīng)兩條染色體的三代測序數(shù)據(jù)挑出，比對到兩條染色體上查看深度场躯，以及挑出這部分序列重新拼接，拿到完整的染色體旅挤。
#!/bin/bash
# 設(shè)置比對基因組和HiFi數(shù)據(jù)的路徑
contaminant_genome="/path2fasta"
hifi_reads="/path2fastq"
output_dir="/home/output"
line="HIFI"
# 創(chuàng)建輸出目錄
mkdir -p $output_dir/map2out
mkdir -p $output_dir/clean_hifi
# 1. 創(chuàng)建索引
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 文件轉(zhuǎn)換為 FASTQ 格式
samtools bam2fq $output_dir/clean_hifi/${line}.scaffold_specific.bam > $output_dir/clean_hifi/${line}_scaffold_specific_clean.fastq
