寫在前面最近主要忙一些植物群體基因組數據的項目,。前面提過，趕時間,，全基因組的 SNP Calling 使用 GATK 流程,，還是需要跑上兩三天。但這個還是耗時,，并不一定能夠趕上工期,。于是我將目標轉移到葉綠體基因組。我們都很清楚,，葉綠體基因組在植物中極其保守,，尤其是同一科屬內，更不提同一物種的不同亞類,。我們組裝了關注物種的一個葉綠體基因組后,，即完全可以進行相關 SNP 鑒定。 http://samtools.sourceforge.net/mpileup.shtml讀段回帖,、排序與去重復安裝必要的程序，除非環(huán)境中已經有相關程序 conda install bwaconda install samtoolsconda install bcftools注意,，后兩者安裝后,，建議直接再 Update 一下 conda update samtoolsconda 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 # 發(fā)現 bcftools 也有 mplieup 功能,，有意思bcftools mpileup -Ou --threads 40 -f 108.ref.cpGenome.fa *.mkdup.bam|bcftools call -vm --ploidy 1 --threads 40 > direct.vcf進行 PCA 分析 # 考慮調整 --mafplink --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.zipunzip 1.44.zipchmod -R a+x VCF2Dis-1.44/./VCF2Dis-1.44/bin/VCF2Dis -InPut direct.vcf -OutPut sample.vcf.dist直接到 FastME2 網站工具上,，基于距離矩陣,，繪制進化樹 寫在最后整體流程簡單，不做贅述,。有時候想想,，有些東西還是記錄下來，以免后續(xù)用到時,，還得翻翻舊硬盤,，找找流程。 |
|
|
