一、轉(zhuǎn)載https://blog.csdn.net/yangl7/article/details/109323644
vcftools --vcf AxiomGT1.calls.vcf --window-pi 1000 --window-pi-step 1000 --out GT1_pi
生成兩個(gè)文件
GT1_pi.windowed.pi
GT1_pi.log
使用GT1_pi.windowed.p文件通過R進(jìn)行繪圖
library(ggplot2)
data<-read.table("GT1_pi.windowed.pi",header=T)
pdf("pi_result.pdf",width=15,height=3)
p <- ggplot(sc3,aes(x=BIN_END/1000000,y=PI,color='N_VARIANTS')) + geom_line(size=0.5) + xlab("Chromosome1 (Mb)")+ ylab("Pi")
p + theme_bw()
dev.off()
數(shù)據(jù)的‘N_VARIANTS’應(yīng)該是有兩個(gè)因子（1,2），所以實(shí)際應(yīng)該對(duì)應(yīng)兩條折線，我也使用過color=‘factor(N_VARIANTS)’,還是只畫出了一條線，如果有大神能指點(diǎn)一二我將不勝感激。
二、轉(zhuǎn)載https://zhuanlan.zhihu.com/p/52064863
利用VCFtools計(jì)算幾個(gè)群體遺傳學(xué)里的參數(shù)，然后利用R將結(jié)果呈現(xiàn)出來。
Tajima's D計(jì)算，
這個(gè)是選擇相關(guān)的一個(gè)參數(shù)，大于0代表群體觀測雜合度高于預(yù)期雜合度，稀有等位基因頻率降低（群體收縮或者平衡選擇），小于0說明群體觀測雜合位點(diǎn)少于預(yù)期值，稀有等位基因頻率增加（群體擴(kuò)張或者低頻選擇）。
也就是說，只有0是正常的，其他都是選擇發(fā)生。
vcftools --vcf 20F.vcf --TajimaD 500000 --out TajimaD
2. π，核苷酸多樣性，越大說明核苷酸多樣性越高，越低說明兩個(gè)座位DNA序列差異越小。
vcftools --vcf 20F.vcf --window-pi 500000 --out pi
3. Fst, 分化系數(shù)，從0到1說明親緣關(guān)系越來越遠(yuǎn)。接近于0說明兩個(gè)個(gè)體親緣關(guān)系近，接近1說明親緣關(guān)系遠(yuǎn)。10number.txt和20number.txt是想要分組的個(gè)體文件，每行一個(gè)個(gè)體，跟vcf中的個(gè)體名稱一致。
vcftools
--vcf 30results/30results_SNP.vcf --weir-fst-pop 20data/10number.txt
--weir-fst-pop 20data/20number.txt --out 20data/P_VS_Y --fst-window-size
500000
4. Hardy-Weinberg平衡檢測，這個(gè)主要是檢測基因型頻率是否等于基因頻率乘積。比如A:0.3,a:0.7那么Aa的頻率是否為0.42
vcftools --vcf 20data/20F.vcf --hardy --out 20data/HW
library(ggplot2)
mydata<-read.table("TajimaD500K.Tajima.D",header = T) #其他文件只需要換一下名字就好
mydata1<-na.omit(mydata)
ggplot(mydata1, aes(x=BIN_START/1000000,y=TajimaD,group=factor(CHROM),colour=CHROM))+geom_line()+facet_wrap(mydata1$CHROM)+xlab("Physical distance(Mb)")+ylab("Tajima's D")+theme(legend.position = "none")
a<-mydata1[sample(nrow(mydata1), 1000), ] #隨機(jī)選取其中的1000行數(shù)據(jù)，這個(gè)Hardy-Weinberg平衡文件需要的，因?yàn)榻Y(jié)果太多，無法利用ggplot繪出全部結(jié)果
三、轉(zhuǎn)載http://www.itdecent.cn/p/b73a8d6233be
Fst的計(jì)算原理與實(shí)戰(zhàn)
概念回顧
Fst：群體間遺傳分化指數(shù)，是種群分化和遺傳距離的一種衡量方法，分化指數(shù)越大，差異越大。適用于亞群體間多樣性的比較。
用于衡量種群分化程度，取值從0到1，為0則認(rèn)為兩個(gè)種群間是隨機(jī)交配的，基因型完全相似；為1則表示是完全隔離的，完全不相似。它往往從基因的多樣性來估計(jì)，比如SNP或者microsatellites(串聯(lián)重復(fù)序列一種，長度小于等于10bp)。是一種以哈溫平衡為前提的種群遺傳學(xué)統(tǒng)計(jì)方法。
兩個(gè)種群之間遺傳差異的基本測量是統(tǒng)計(jì)量FST。在遺傳學(xué)中，F(xiàn)一詞通常代表“近親繁殖”，它傾向于減少群體中的遺傳變異。遺傳變異可以用雜合度來衡量，所以F一般表示群體中雜合性的減少。 FST是與它們所屬的總?cè)后w相比，亞群體中雜合性的減少量。
具體可以下面的公式表示：
Fst= (Ht-Hs）/ Ht
Hs:亞群體中的平均雜合度
Ht:復(fù)合群體中的平均雜合度
理論上計(jì)算Fst的步驟
理論上要估算FST，需要以下步驟：
找出每個(gè)亞群的等位基因頻率。
查找復(fù)合群體的平均等位基因頻率
計(jì)算每個(gè)亞群的雜合度（2pq）
計(jì)算這些亞群雜合度的平均值，這是HS。
根據(jù)總體等位基因頻率計(jì)算雜合度，這是HT。
最后，計(jì)算FST =（HT-HS）/ HT
具體例子
基因SLC24A5是黑色素表達(dá)途徑的關(guān)鍵部分，其導(dǎo)致皮膚和毛發(fā)色素沉著。與歐洲較輕的皮膚色素密切相關(guān)的SNP是rs1426654。 SNP有兩個(gè)等位基因A和G，其中G與輕度皮膚相關(guān)，在猶他州的歐裔美國人中，頻率為100％。美洲印第安人與美國印第安人混血兒的SNP在頻率上有所不同。墨西哥的樣本有38％A和62％G;在波多黎各，頻率分別為59％A和41％G，查爾斯頓的非裔美國人樣本中有19％A和81％G.這個(gè)例子中的FST是什么？
手工計(jì)算的步驟如下：
在得到每個(gè)群體中等位基因的頻率后:
首先第一步是計(jì)算每個(gè)亞群的雜合度H，這一步簡單H=2pq （p和q是對(duì)應(yīng)群體中等位基因的頻率）
第二步計(jì)算這些亞群雜合度的平均值HS，取三個(gè)不同種群中的H值，除去平均數(shù)3，得到結(jié)果
第三步，計(jì)算群體等位基因頻率，和上述方法類似，取每個(gè)等位基因在對(duì)應(yīng)種群中的頻率疊加除去平均數(shù)。
第四步，計(jì)算所有群體中的雜合度HT，使用2pq，p和q是群體中的等位基因頻率，在第三步驟已經(jīng)算好了。
最后，根據(jù)計(jì)算公式算出Fst。
在實(shí)際中，當(dāng)然不可能用手工對(duì)每個(gè)位點(diǎn)進(jìn)行這樣的計(jì)算，怎樣在電腦中分析FST等值，會(huì)在下次內(nèi)容中介紹。
實(shí)戰(zhàn)
當(dāng)然了，實(shí)際中我們并不需要像上述例子那樣逐個(gè)計(jì)算。等位基因頻率的信息都藏在snp calling后的vcf中，再使用恰當(dāng)?shù)墓ぞ呔涂梢钥焖儆?jì)算出Fst。下面給大家介紹一種比較popular的計(jì)算Fst的方法（還有其他方法不限于一種），使用vcftools：
##對(duì)每一個(gè)SNP變異位點(diǎn)進(jìn)行計(jì)算
vcftools--vcf test.vcf--weir-fst-pop1_population.txt--weir-fst-pop2_population.txt--outp_1_2—single
##按照區(qū)域來計(jì)算vcftools--vcf test.vcf--weir-fst-pop1_population.txt--weir-fst-pop2_population.txt--outp_1_2_bin--fst-window-size500000--fst-window-step50000
# test.vcf是SNP calling 過濾后生成的vcf 文件；
# p_1_2_3 生成結(jié)果的prefix
# 1_population.txt是一個(gè)文件包含同一個(gè)群體中所有個(gè)體，一般每行一個(gè)個(gè)體。個(gè)體名字要和vcf的名字對(duì)應(yīng)。
# 2_population.txt 包含了群體二中所有個(gè)體。
#計(jì)算的窗口是500kb，而步長是50kb （根據(jù)你的需其可以作出調(diào)整）。我們也可以只計(jì)算每個(gè)點(diǎn)的Fst，去掉參數(shù)（--fst-window-size 500000 --fst-window-step 50000）即可。
分別對(duì)不同結(jié)果進(jìn)行圖形繪制
##圖1
library(ggplot2)
data<-read.table("test1.out.windowed.weir.fst",header=T)
sc3=subset(data,CHROM=="Gm01")
p<-ggplot(sc3,aes(x=BIN_END/1000000,y=WEIGHTED_FST))+geom_point(size=0.5,colour="blue")+xlab("Physical distance (Mb)")+ylab("Fst")+ylim(-1,1)
p+theme_bw()
##圖2
library(ggplot2)
data<-read.table("test.out.weir.fst",header=T)
sc3=subset(data,CHROM=="Gm01")p<-ggplot(sc3,aes(x=POS,y=WEIR_AND_COCKERHAM_FST))+geom_point(size=0.5,colour="blue")+xlab("Physical distance (Mb)")+ylab("Fst")+ylim(-1,1)
p+theme_bw()
可以看到兩幅圖還是基本對(duì)應(yīng)的，然后對(duì)一些接近1，與大部分點(diǎn)偏離比較高的點(diǎn)可以將其與功能注釋相結(jié)合，還有一些選擇壓力分析的工具相結(jié)合，尋找出其對(duì)應(yīng)的基因，觀察該基因是否是被選擇。Pi也是選擇分析中一個(gè)很要的參數(shù)，這一部分內(nèi)容有時(shí)間這一點(diǎn)以后再補(bǔ)回來。
四、轉(zhuǎn)載于https://mp.weixin.qq.com/s?__biz=MzU1NDkzOTk2MQ==&mid=2247485064&idx=1&sn=6a9d4d29d5be747a38c1d4e8a6e6113b&chksm=fbdaa5deccad2cc8602ebc2f587547828a12dc4c6df1f360ff14e6c0628b785cc725ef7e5181&scene=21#wechat_redirect
群體分化指數(shù)-Fst
1、進(jìn)行質(zhì)控，剔除高缺失率（--geno 0.05）和極低等位基因頻率( --maf 0.01 )的SNP
#!/bin/bash
plink=/software/biosoft/software/plink/plink
$plink?--vcf?60dog.vcf?--geno?0.05?--maf?0.01?--dog?--recode?vcf-iid?--out?60dog_QC
2、這里以TM和YJ為例，計(jì)算FST，因?yàn)橹饌€(gè)位點(diǎn)計(jì)算FST時(shí)，可能會(huì)出現(xiàn)FST值很高的假陽性信號(hào)（中性選擇導(dǎo)致），所以這里考慮到搭載效應(yīng)同時(shí)計(jì)算了滑窗FST，二者可以對(duì)照著看。TM_4500.txt為TM樣本ID文件，一行一個(gè)ID，YJ_800.txt同左。
這里可以來看下結(jié)果文件
單位點(diǎn)：
窗口：
3、可視化
這里使用R包qqman完成曼哈頓圖可視化，以窗口結(jié)果為例。先處理數(shù)據(jù)文件格式：
sed?"1d"?TM_4500vsYJ_800_20k_5k.windowed.weir.fst|awk?'{if($5<0)print?$1"\t"$2"\t0";else?print?$1"\t"$2"\t"$5}'?>?TM_4500vsYJ_800_20k_5k_plot.txt
以下是可視化腳本代碼
###?加載要使用的R包
library(qqman)
library(Cairo)
Fstfile<-read.table("TM_4500vsYJ_800_20k_5k_plot.txt",?header=F,?stringsAsFactors=F)
SNP<-paste(Fstfile[,1],Fstfile[,2],sep?=?":")
Fstfile=cbind(SNP,Fstfile)
colnames(Fstfile)<-c("SNP","CHR","POS","Fst")
outfile<-"TM_4500vsYJ_800_20k_5k"
filePNG<-paste(outfile,"manhattan.png",sep=".")
CairoPNG(file=filePNG,width=1500,height=500)
colorset<-c("#FF0000","#FFD700","#2E8B57","#7FFFAA","#6495ED","#0000FF","#FF00FF")
###?調(diào)用qqman,畫出全部染色體
manhattan(Fstfile),chr="CHR",bp="POS",p="Fst",snp="SNP",?col=colorset,logp=F,suggestiveline?=?F,?genomewideline?=?F,ylab="Fst",ylim=c(0,1),font.lab=4,cex.lab=1.2,main="TMvsYJ",cex=0.8)
###?只繪制單條染色體
if?(FALSE){
????manhattan(subset(Fstfile,CHR=="10"),chr="CHR",bp="POS",p="Fst",snp="SNP",?col=colorset,logp=F,suggestiveline?=?F,?genomewideline?=?F,ylab="Fst",ylim=c(0,1),font.lab=4,cex.lab=1.2,"TMvsYJ_chr10",cex=0.8)
}
dev.off()
print("==================manhattan?plot====================")
可視化結(jié)果展示
全部染色體
10號(hào)染色體
如果發(fā)現(xiàn)更好的教程，后續(xù)更新的哦！歡迎一起討論！
