写在前面
最近主要忙一些植物群体基因组数据的项目。前面提过，赶时间，全基因组的 SNP Calling 使用 GATK 流程，还是需要跑上两三天。但这个还是耗时，并不一定能够赶上工期。于是我将目标转移到叶绿体基因组。我们都很清楚，叶绿体基因组在植物中极其保守，尤其是同一科属内，更不提同一物种的不同亚类。我们组装了关注物种的一个叶绿体基因组后，即完全可以进行相关 SNP 鉴定。
我关注的物种，叶绿体基因组组装出来，大小是 170Kb，说实话，相对于 600Mb 的基因组来说，确实很小。SNP Calling 也会很快。
但为了加速流程，我选择了最简单粗暴的方式，使用 samtools + bcftools，具体可参考流程（发现这个流程时，我还是有点激动，毕竟是官方文档）。其实我刚接触生信的之后（2015年），那会不知道是否有 bcftools 或者 SNP calling 的功能，我是直接手动解 samtools 的 mpileup 输出。
http://samtools.sourceforge.net/mpileup.shtml
读段回帖、排序与去重复
安装必要的程序，除非环境中已经有相关程序
conda install bwa
conda install samtools
conda install bcftools
注意，后两者安装后，建议直接再 Update 一下
conda update samtools
conda update bcftools
下载并准备 Picard 软件，用于标记重复，事实上，似乎可以直接使用 sambamba 软件
wget https://github.com/broadinstitute/picard/releases/download/2.25.7/picard.jar
建立索引，需要注意 - GetOrganelle 输出的基因序列 ID 比较复杂 - 应该简化
# cpGenome.fa 即叶绿体基因组
bwa index cpGenome.fa
序列比对到叶绿体基因组上
# 将双端测序数据，比对到 叶绿体基因组
ls *.m1.fq|perl -lpe 's/.m1.fq//'|parallel -j 10 "
bwa mem -t 4 -M cpGenome.fa {}.m1.fq {}.m2.fq 2>{}_map.log | samtools sort-O bam -@ 4 -o {}.sorted.bam 1>{}.sort.log 2>&1
"
标记重复
ls *.m1.fq|perl -lpe 's/.m1.fq//'|parallel -j 10 "
java -jar picard.jar MarkDuplicates \
I={}.sorted.bam \
O={}.sorted.mkdup.bam \
M={}.sorted.mkdup.metrics 1>{}_mkdup.log 2>&1
"
Call SNP
# 发现 bcftools 也有 mplieup 功能，有意思
bcftools mpileup -Ou --threads 40 -f 108.ref.cpGenome.fa *.mkdup.bam|bcftools call -vm --ploidy 1 --threads 40 > direct.vcf
进行 PCA 分析
# 考虑调整 --maf
plink --maf 0.02 --allow-extra-chr --vcf direct.vcf --pca header tabs -out plink.stat
# 提取 PC1 和 PC2，使用 R 或者其他工具进行散点图可视化即可
perl -lpe 's/.sorted.mkdup.bam//g' plink.stat.eigenvec|cut -f1,3,4 > pca.plot.tab
绘制进化树
# 使用 VCF2Dis 软件快速计算样品距离矩阵
wget https://github.com/BGI-shenzhen/VCF2Dis/archive/refs/tags/1.44.zip
unzip 1.44.zip
chmod -R a+x VCF2Dis-1.44/
./VCF2Dis-1.44/bin/VCF2Dis -InPut direct.vcf -OutPut sample.vcf.dist
直接到 FastME 网站工具上，基于距离矩阵，绘制进化树
写在最后
整体流程简单，不做赘述。有时候想想，有些东西还是记录下来，以免后续用到时，还得翻翻旧硬盘，找找流程。
当然，更多时候，似乎有新的软件和工具可以使用了。
