寫在前面
最近主要忙一些植物群體基因組數(shù)據(jù)的項目。前面提過，趕時間，全基因組的 SNP Calling 使用 GATK 流程，還是需要跑上兩三天。但這個還是耗時，并不一定能夠趕上工期。于是我將目標轉(zhuǎn)移到葉綠體基因組。我們都很清楚，葉綠體基因組在植物中極其保守，尤其是同一科屬內(nèi)，更不提同一物種的不同亞類。我們組裝了關注物種的一個葉綠體基因組后，即完全可以進行相關 SNP 鑒定。
我關注的物種，葉綠體基因組組裝出來，大小是 170Kb，說實話，相對于 600Mb 的基因組來說，確實很小。SNP Calling 也會很快。
但為了加速流程，我選擇了最簡單粗暴的方式，使用 samtools + bcftools，具體可參考流程（發(fā)現(xiàn)這個流程時，我還是有點激動，畢竟是官方文檔）。其實我剛接觸生信的之后（2015年），那會不知道是否有 bcftools 或者 SNP calling 的功能，我是直接手動解 samtools 的 mpileup 輸出。
http://samtools.sourceforge.net/mpileup.shtml
讀段回帖、排序與去重復
安裝必要的程序，除非環(huán)境中已經(jīng)有相關程序
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
序列比對到葉綠體基因組上
# 將雙端測序數(shù)據(jù)，比對到 葉綠體基因組
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
# 發(fā)現(xiàn) bcftools 也有 mplieup 功能，有意思
bcftools mpileup -Ou --threads 40 -f 108.ref.cpGenome.fa *.mkdup.bam|bcftools call -vm --ploidy 1 --threads 40 > direct.vcf
進行 PCA 分析
# 考慮調(diào)整 --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 網(wǎng)站工具上，基于距離矩陣，繪制進化樹
寫在最后
整體流程簡單，不做贅述。有時候想想，有些東西還是記錄下來，以免后續(xù)用到時，還得翻翻舊硬盤，找找流程。
當然，更多時候，似乎有新的軟件和工具可以使用了。
