前言
最近分析叶绿体基因组，进行SNP Calling，网上搜到生信药丸的简明教程 " 教程 | 简单粗暴的叶绿体基因组 SNP Calling 流程"，整个流程跑了一遍基本没问题，把流程记录一下，以备后续查询方便。
首先，安装相关软件：
conda install -c bioconda bwa
conda install -c bioconda picard
conda install -c bioconda samtools
conda install -c bioconda bcftools
conda install -c bioconda plink
1. 建立索引 cpGenome.fa 即叶绿体基因组
bwa index cpGenome.fa
2. 将双端测序数据，比对到 叶绿体基因组
nohup bwa mem -t 4 -M cpGenome.fa CRRXXX_f1.fastq.gz CRRXXX_r2.fastq.gz 2>CRRXXX_map.log|samtools sort -O bam -@ 4 -o CRRXXX.sorted.bam 1>CRRXXX.sort.log 2>&1 &
多任务批量执行可以参照大神的简明教程利用parallel 命令运行：
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 "
我在运行parallel的时候出现了报错，可能软件没安装好就换了种方法：首先ls，grep，sed生成列表，利用awk把脚本打印出来，然后执行sh：
ls |grep ‘_1.fq.gz’|sed 's/_1.fq.gz//g'|awk '{print "nohup bwa mem -t 4 -M cpGenome.fa "$1"_1.fq.gz "$1"_2.fq.gz 2>"$1"_map.log|samtools sort -O bam -@ 4 -o "$1".sorted.bam 1>"$1".sort.log 2>&1 &" }'|sh
3. 标记重复
nohup picard MarkDuplicates -I CRRXXX.sorted.bam -O CRRXXX.sorted.mkdup.bam -M CRRXXX.sorted.mkdup.metrics 1>CRRXXX_mkdup.log 2>&1 &
同样，首先ls，grep生成输入列表，利用awk把脚本打印出来，然后执行sh：
$ ls|grep sorted.bam|awk -F "." '{print "nohup picard MarkDuplicates -I "$0" -O "$1".sorted.mkdup.bam -M "$1".sorted.mkdup.metrics 1>"$1"_mkdup.log 2>&1 &"}' |sh
4. SNP Calling
bcftools mpileup -Ou --threads 40 -f cpGenome.fa *.mkdup.bam|bcftools call -vm --ploidy 1 --threads 40 > direct.vcf
注意：命令加nohup &后台运行如果报错，可以把命令放.sh脚本里面然后在加nohup &后台运行
5. 进行 PCA 分析绘图
plink --maf 0.02 --allow-extra-chr --vcf direct.vcf --pca header tabs -out plink.stat
plink运行生成plink.stat.eigenval（可以将对应值添加到坐标轴，如：PC1(55.4%)）和plink.stat.eigenvec两个文件，在eigenvec后面增加一列Group分组信息，在excel里面保存为csv文件：Pca.csv，利用ggplot2 绘图，代码：
Pca <- read.csv("Pca.csv",header = T,row.names = 1)
p = ggplot(Pca, aes(PC1, PC2, color = Group))
p <- p + geom_point(size=2)
p <- p + stat_ellipse(level = 0.95, size = 1)
p <- p + geom_hline(yintercept = 0)
p <- p + geom_vline(xintercept = 0)
p <- p + theme_bw()
p <- p + geom_text(label = Pca$Group)
p <- p + labs(x = "PC1(55.4%)", y = "PC2(22.6%)")
p
PCA.png
6.绘制进化树
使用 VCF2Dis 软件快速计算样品距离矩阵,直接到 FastME网页基于距离矩阵，直接输入sample.vcf.dist文件绘制进化树
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结果可以下载.nwk树文件利用本地软件查看树，也可以点击链接在线查看：
FastME.png
Tree.png
7.admixture分析
plink --vcf direct.vcf --recode --allow-extra-chr --out abc #vcf文件转换成ped文件，生成abc.ped;
plink --file abc --make-bed --allow-extra-chr --out abcd #ped文件转换成bed文件，生成abcd.bed,abcd.bim,abcd.fam，注意要将abcd.bim第一列变成整数，不然admixture会报错;
for K in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15;do admixture --cv direct.vcf.g5mac3.bed.bed $K|tee log${K}.out;done
grep -h CV log*.out #选择合适的K值
R作图：
tbl=read.table("direct.vcf.Q")
barplot(t(as.matrix(tbl)), col=rainbow(7),xlab="Individual #", ylab="Ancestry", border=NA, space=0)
admixture.png
