利用purge_haplotigs去杂合 - 简书

1. 安装软件
   1. conda安装的部分

conda create -n purha -y

conda install purge_haplotigs -y
conda install bedtools -y
conda install samtools -y
conda install minimap2 -y
conda install MUMmer -y

1. conda解决不了的部分

$ sudo apt install r-base r-base-dev

# on a new install we wont have the required R library 'ggplot2' installed
$ sudo su - -c "R -e \"install.packages('ggplot2', repos='http://cran.rstudio.com/')\""

sudo apt-get install sra-toolkit #下载处理sra数据的

1. 数据准备
   1. 组装得到的基因组文件（genome.fasta）
   2. 用于组装的三代测序文件（subreads.fasta.gz）或者是二代Ilumilina数据也可以
   3. 用上面两个文件比对的到的aligned.bam文件
2. 比对得到bam文件
   1. hifi比对
      1. a：sam格式的输出
      2. x：指定比对文件的格式，这里的是PacBio HiFi 数据

nohup minimap2 -ax map-hifi pulchrum_asm_third.hic.hap1.p_ctg.fa \
    leaf8.ccs.fasta.gz \
    m84128_s4.fq.gz \
    m84288_s1.fq.gz \
    -o hifi.sam \
    -t 10 &

1. 过滤

nohup samtools view -hF 256 hifi.sam -@ 20 -o hifi.filtered.sam &

1. 排序

nohup samtools sort -@ 20 -m 1G -o aligned.bam -T tmp.ali hifi.filtered.sam &

1. 通过运行readhist脚本生成覆盖率直方图。
   1. 该脚本将生成一个直方图png图像文件供您查看，并生成一个BEDTools 'genomecov'输出文件供您在步骤2中使用。
   2. -b aligned.bam：指定输入的 BAM 文件，该文件是通过 minimap2 和 samtools 处理过的排序后比对文件，包含 PacBio 子读与参考基因组的比对结果。
   3. -g genome.fasta：指定参考基因组的 FASTA 文件。
   4. 你应该会得到一个双峰直方图——一个单倍体水平的峰值，一个二倍体水平的峰值。 如果你使用的是 phased assembly，二倍体峰值可能会非常小。

nohup purge_haplotigs readhist -b aligned.bam -g pulchrum_asm_third.hic.hap1.p_ctg.fa -t 20 &

1. 运行contigcov脚本，使用上一步的 cutoffs值(就是那三个点的数值)来分析重叠的覆盖。
   1. hap1和ctg的相同

nohup purge_haplotigs contigcov \
  -i aligned.bam.200.gencov \
  -l 15 \
  -m 75 \
  -h 185 \
  -o coverage_stats.csv \
  -j 80 \
  -s 80 &

1. 运行purge脚本， 这个脚本将自动运行一个BEDTools窗口覆盖分析，一个blast搜索，并执行lastz比对等，以评估哪些contigs重新分配及哪些contigs保留。
   1. -a 70：设置 70% 的对比度，用于识别并标记 haplotigs（次要异质基因组）。

nohup purge_haplotigs purge \
  -g pulchrum_asm_third.hic.hap1.p_ctg.fa \
  -c coverage_stats.csv \
  -b aligned.bam \
  -t 20 \
  -a 70 & 

1. 最终生成5个文件

.artefacts.fasta：无用的contig，也就是没有足够覆盖度的contig.
.fasta：新的单倍型组装
.haplotigs.fasta：从原本组装分出来的haplotigs
.reassignments.tsv: 单倍型的分配信息
.contig_associations.log: 运行日志。

1. quast评估
   1. 参数m：设置 contig （基因组片段）的 最小长度阈值。，低于这个数值的会被忽略

nohup quast pulchrum_asm_third.hic.hap1.p_ctg_curated.fa pulchrum_asm_third.hic.p_ctg_curated.fa -o quast -m 0 -t 20 &

1. busco评估

nohup busco -i curated.fasta \
  -c 20 \
  -m geno \
  --out curated_hap1 \
  -l /home/server/database/busco_database/embryophyta_odb12 \
  --offline \
  -f &