PSMC軟件分析群體歷史有效群體大小流程
1 文件轉(zhuǎn)換
基因組文件格式為.bam耐朴，callsnp之后，再轉(zhuǎn)換為.fq.gz格式才能做為PSMC輸入文件
bcftools mpileup 和 bcftools call可以進(jìn)行SNP calling颜曾，
bcftools mpileup -Ou -I -f ref.fa .sorted.mardup.bam | bcftools call -c -Ov | vcfutiles vcf2fq -d 10 -D 100 | gzip > dilpoid.fq.gz
需要注意的是mpileup命令雖然也會輸出vcf格式文件褐健，但是并不直接進(jìn)行snp calling钳幅。
下面的命令可以生成vcf格式文件
bcftools mpileup -I -f ref.fa .sorted.markdup.bam > mpileup.vcf
直接 生成的vcf文件內(nèi)容如下：
里面的每一條記錄并不是一個SNP位點(diǎn)物蝙，而是染色體上每個堿基的比對情況匯總，這種信息官方稱為genotype likelihoods敢艰。
call命令才是真正執(zhí)行SNP calling的程序诬乞，基本用法如下
bcftools call mpileup.vcf -c -Ov -o variants.vcf
在進(jìn)行SNP calling 時(shí)，必須選擇一種算法钠导，有兩種calling算法可供選擇震嫉，分別是-c 和 -m 參數(shù)。-c參數(shù)對應(yīng)consensus-caller算法牡属，-m參數(shù)對應(yīng)multiallelic-caller算法票堵，后者更適合多種allel和罕見變異的calling。
2 生成psmc.fa和psmc文件
utils/fq2psmcfa -q20 diploid.fq.gz > diploid.psmcfa
psmc -N25 -t15 -r5 -p "4+25*2+4+6" -o diploid.psmc diploid.psmcfa
-p : 64個時(shí)間間隔逮栅，在計(jì)算時(shí)需要估計(jì)28個參數(shù)悴势，將64個時(shí)間間隔劃分一下窗宇，第一個劃分了4個時(shí)間間隔，在這四個時(shí)間間隔估計(jì)1個Ne的參數(shù)特纤；25*2劃分了25組2個時(shí)間間隔军俊，每兩個時(shí)間間隔估計(jì)一個和Ne有關(guān)的參數(shù)。
3 畫圖
生成psmc文件之后開始畫圖
perl utils/psmc_plot.pl -u 2e-09 -g 1 -p sample_plot sample.psmc
-u為突變率捧存，-g為世代間隔
-u 的說明是 absolute mutation rate per nucleotide [2.5e-08]
注意粪躬，這里填的是per generation per site的mutation，不是per year per site 的mutation昔穴，否則會導(dǎo)致PSMC曲線平移一個數(shù)量級镰官。
4 不同物種繪制在同一個圖中
cat sample1.psmc sample2.psmc sample3.psmc > combine.psmc
psmc_plot.pl -u mutation_rate -g generation -M 'sample1,sample2,sample3' combine_plot combine.psmc
-u 一般根據(jù)文獻(xiàn)來看用什么，不過突變率一般哺乳動物不會相差很大傻咖，-g時(shí)代間隔就是從個體出生到他的孩子出生的時(shí)間朋魔，可以根據(jù)經(jīng)驗(yàn)定
