1.選擇信號分析—Fst選擇信號分析
概念
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
鏈接：Fst的計(jì)算原理與實(shí)戰(zhàn) - 簡書 (jianshu.com)
- FST(遺傳分化指數(shù)):計(jì)算使用vcftools,可視化分為箱線圖和散點(diǎn)圖，單組比較使用在染色體上的散點(diǎn)圖像樊，多組比較使用箱線圖旅敷。
- FST的原理生棍，計(jì)算方法媳谁，可視化的方法 http://www.reibang.com/p/bb0beec0ed63
http://www.reibang.com/p/c520c36fe340 - vcftools使用說明： https://www.dazhuanlan.com/2020/03/04/5e5f850672ddf/
VCFTOOLS官網(wǎng),此步驟用于計(jì)算 Weir 和 Cockerham 1984 年論文中的 Fst 估計(jì)值，這是 Fst 的首選計(jì)算方法柔纵。
兩篇少樣本數(shù)量不影響Fst文獻(xiàn)
兩篇少樣本數(shù)量不影響Fst文獻(xiàn)
FST計(jì)算方法
- Fst值的取值范圍是[0,1]锤躁，最大值為1表明兩個(gè)群體完全分化，最小值為0表明群體間無分化。
- 實(shí)際使用FST<0--0.05,表示群體分化很屑釉怠；0.05--0.15沈贝，中等程度的分化；0.15--0.25宋下，較大的分化辑莫；0.25以上，分化很大各吨。
- 其實(shí)Fst分析就是看兩個(gè)群體之間分化程度的一種方法，F(xiàn)st值越大（越接近1）表明兩個(gè)群體間分化程度越高揭蜒，親緣關(guān)系越遠(yuǎn)；Fst值越刑敫（越接近0）表明群體間分化程度越低，親緣關(guān)系越近欺冀。
- 由于負(fù)值Fst無意義萨脑，所以通常計(jì)算出的負(fù)值Fst需要將其替換為0
按照計(jì)算方法可分為：按照SNP單點(diǎn)計(jì)算和滑動窗口模式計(jì)算隐轩，一般使用的是滑動窗模式
vcftools計(jì)算Fst實(shí)戰(zhàn)
SNP單點(diǎn)計(jì)算
##對每一個(gè)SNP變異位點(diǎn)進(jìn)行計(jì)算
vcftools --vcf sample.vcf --weir-fst-pop population_1.txt --weir-fst-pop population_2.txt --out sample_1_2
? ? ? ?
滑動窗口模式計(jì)算
#按照窗口模式計(jì)算（更準(zhǔn)確）
vcftools --vcf sample.vcf --weir-fst-pop population_1.txt --weir-fst-pop population_2.txt --out P_1_2 --fst-window-size 500000 --fst-window-step 50000
# sample.vcf是SNP calling 或芯片數(shù)據(jù)過濾后生成的vcf 文件龙助；
# p_1_2_3 生成結(jié)果的prefix
# population_1.txt是一個(gè)文件包含同一個(gè)群體中所有個(gè)體砰奕，一般每行一個(gè)個(gè)體。個(gè)體名字要和vcf的名字對應(yīng)仅淑。
# population_2.txt 包含了群體二中所有個(gè)體胸哥。
##--fst-window-size # 設(shè)置計(jì)算Fst的窗口大小涯竟，根據(jù)自己的數(shù)據(jù)進(jìn)行設(shè)置，看看別人文章里怎么用的
##--fst-window-step # 設(shè)置計(jì)算Fst的步長長度银酬，根據(jù)自己的數(shù)據(jù)進(jìn)行設(shè)置
#計(jì)算的窗口是500kb筐钟，而步長是50kb （根據(jù)你的需其可以作出調(diào)整）揩瞪。我們也可以只計(jì)算每個(gè)點(diǎn)的Fst篓冲，去掉參數(shù)（--fst-window-size 500000 --fst-window-step 50000）即可。
#注：建議參考已發(fā)表文獻(xiàn)壹将，根據(jù)相關(guān)物種的在基因組單位區(qū)段面積的位點(diǎn)數(shù)目、LD平均值或在r2=0.1的時(shí)候位置來定窗口大小
#根據(jù)加權(quán)fst為Fst排序
sort -k 5 -g p1_p2_window.windowed.weir.fst > p1_p2.sorted.fst
#窗口計(jì)數(shù)
wc -l p1_p2.sorted.fst #假設(shè)為20000
#取前1%
tail -n 200 p1_p2.sorted.fst > p1_p2.sorted.1%.fst
#找出前1%中最小的加權(quán)fst值
參數(shù)說明：
--vcf 輸入vcf文件
--weir-fst-pop 輸入群體的群體ID名妇菱，該文件必須是txt格式惊畏，每個(gè)ID占一行。
--fst-window-size 500000 --fst-window-step 50000 颜启，這里窗口設(shè)置為500kb，步長設(shè)置為50kb缰盏。根據(jù)情況調(diào)整窗口大小。
(18條消息) 【群體遺傳】Fst（群體間分化指數(shù)）_陳有樸的博客-CSDN博客_fst遺傳分化指數(shù)
- population_2.txt和population_2.txt格式一樣负溪，只有一列樣品信息济炎，個(gè)體名字要和vcf的名字對應(yīng)**
例如：
Sample1
Sample2
Sample3
分別對不同結(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()
計(jì)算haploPS、XP-EHH须尚、 Fst，正向選擇分析方法尋找性狀相關(guān)的位點(diǎn)
<meta charset="utf-8">
π的計(jì)算
π耐床，核苷酸多樣性，越大說明核苷酸多樣性越高胯甩，越低說明兩個(gè)座位DNA序列差異越小。
vcftools
LD-blockde 計(jì)算
計(jì)算和可視化偎箫，參考
分析方法 參考文獻(xiàn)：中國農(nóng)業(yè)科學(xué) scientific reports
找到重測序的數(shù)據(jù)，基因分型淹办，找到單倍型，call SNP,過濾齐遵，注釋snp,可視化，分別計(jì)算haploPS,XP-EHH,Fst.求交集梗摇，公共基因進(jìn)行注釋想许，富集分析伶授。
2.計(jì)算總體PI值和TajimaD
核酸多樣性PI流纹，值越大說明核苷酸多樣性越高
TajimaD用于檢驗(yàn)DNA序列在演化國產(chǎn)過程中是否遵循中性演化模型。D>0:平衡選擇疮蹦，突然群體收縮；D<0:最近的選擇性清除愕乎，最近瓶頸后的群體擴(kuò)張壁公，與消除基因連鎖感论；D=0：群體根據(jù)突變-漂移平衡演變紊册，沒有選擇。
vcftools --vcf Filter.snp.vcf --window-pi 500000 --out total
vcftools --vcf Filter.snp.vcf --TajimaD 500000 --out total
計(jì)算每個(gè)亞群的PI值和TajimaD值
#創(chuàng)建bash腳本
vi test.sh
#寫入以下內(nèi)容：
for i in {1..3};do
#根據(jù)ID芳绩，從vcf文件中提取每個(gè)亞群的信息
vcftools --vcf Filter.snp.vcf --keep ${i}-population.txt --recode --
recode-INFO-all --out p${i}
#根據(jù)提取的信息計(jì)算每個(gè)亞群的PI值
vcftools --vcf p${i}.recode.vcf --out q${i}_pi_500kb --window-pi
500000 --window-pi-step 50000
#根據(jù)提取的信息計(jì)算每個(gè)亞群的TajimaD值
vcftools --vcf p${i}.recode.vcf --out q${i}_500kb.TajimaD --TajimaD
500000
done
#給bash腳本添加可執(zhí)行權(quán)限
chmod +x test.sh
#運(yùn)行腳本
./test.sh
根據(jù)pi值畫箱線圖和小提琴圖
#添加分組信息G1关斜、G2铺浇、G3
data=data.frame(q1_500kb_pi_windowed)
data
Group=c("G1")
data1=data.frame(data,Group)
data=data.frame(q2_500kb_pi_windowed)
data
Group=c("G2")
data2=data.frame(data,Group)
data=data.frame(q3_500kb_pi_windowed)
data
Group=c("G3")
data3=data.frame(data,Group)
#將三組數(shù)據(jù)合并
library(dplyr)
total_data<-dplyr::bind_rows(data1,data2,data3)
#畫箱線圖+小提琴圖
p <- ggplot(total_data,aes(Group,PI,fill=Group))+geom_boxplot(width=.2)+geom_violin()
mytheme <- theme(plot.title = element_text(face = "bold.italic", size = "14", colour = "brown"), axis.title = element_text(face = "bold.italic", size = "10",color = "blue"), axis.text.x = element_text(face = "bold", size = 8, angle = 45, hjust = 1, vjust = 1), axis.text.y = element_text(face = "bold",size = 8), panel.background = element_rect(fill = "white", color = "black"), panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank(), legend.text = element_text(size = 8), legend.title = element_text(size = 10, face = "bold"), panel.grid.minor.x = element_blank())
p + mytheme
根據(jù)PI值畫折線圖
#輸入數(shù)據(jù)時(shí)，每列必須是數(shù)值型數(shù)據(jù)，尤其是染色體那一列
data=data.frame(q1_500kb_pi_windowed)
Group=c("G1")
data1=data.frame(data,Group)
data=data.frame(q2_500kb_pi_windowed)
Group=c("G2")
data2=data.frame(data,Group)
data=data.frame(q3_500kb_pi_windowed)
Group=c("G3")
data3=data.frame(data,Group)
#將三組數(shù)據(jù)合并
library(dplyr)
library(ggplot2)
total_data<-dplyr::bind_rows(data1,data2,data3)
for (i in 1:7){
result_name = paste("chr",i,".png",sep="")
png(result_name,width=1500,height=300)
chrom = subset(total_data,CHROM==i)
xlab = paste("Chromosome",i,"(MB)",sep="")
p <- ggplot(chrom,aes(x=BIN_END/1000000,y=PI,color=Group,shape=Group)) + geom_line(size=0.5) + xlab(xlab)+ ylab("Pi") + theme_bw()
print(p)
dev.off()
}
鏈接：http://www.reibang.com/p/4dbde8533607
參考：http://www.reibang.com/p/c520c36fe340
http://www.reibang.com/p/b73a8d6233be
https://blog.csdn.net/u014182497/article/details/52672308
http://www.omicshare.com/forum/thread-3688-1-1.html
群體中的Fst值-學(xué)習(xí)篇 - 簡書 (jianshu.com)
http://www.reibang.com/p/4dbde8533607
