從獲取下機(jī)數(shù)據(jù)做質(zhì)控所踊、比對(duì)，到calling snp獲得vcf文件這一步概荷，網(wǎng)上已經(jīng)有非常詳細(xì)的教程秕岛，有GATK4.0和全基因組數(shù)據(jù)分析實(shí)踐（上），還有Samtools+bcftools Call SNP等等乍赫。但是在群體遺傳的研究中瓣蛀，我們經(jīng)常需要比較不同群體的vcf文件，鑒定不同群體間基因組水平上發(fā)生分化的區(qū)域雷厂，并注釋其發(fā)生分化的基因惋增，網(wǎng)上對(duì)于這類群體選擇信號(hào)分析的分享還比較少。
我通過(guò)計(jì)算Fst和TajimaD這兩個(gè)經(jīng)典的指數(shù)改鲫，來(lái)進(jìn)行選擇信號(hào)的分析诈皿。話不多說(shuō)，先走一遍流程吧像棘！
1稽亏、工具和文件的準(zhǔn)備
工具：vcftools；snpEff缕题；
文件：突變型群體的vcf文件（DA.vcf）;包含了突變型和野生型兩個(gè)群體的vcf文件（all.vcf）;突變型群體信息（DA.txt）;野生型群體信息（non_DA.txt）;參考基因組的GFF注釋文件（如果GFF注釋文件中沒(méi)有對(duì)應(yīng)的GO編號(hào)信息截歉，則還需要該物種基因ID與GO編號(hào)的對(duì)應(yīng)信息文件annot.go.tab）
2、滑窗計(jì)算Fst和TajimaD值
vcftools可以通過(guò)設(shè)置窗口大小來(lái)計(jì)算染色體(或scaffolds)上指定區(qū)域的Fst和TajimaD的值烟零，但不足的是在計(jì)算TajimaD值時(shí)瘪松，不能設(shè)置步長(zhǎng)（可使用VCF-kit進(jìn)行可設(shè)置步長(zhǎng)的計(jì)算），因此锨阿，為了方便后續(xù)的分析宵睦，我在這里計(jì)算Fst和TajimaD時(shí)，只設(shè)置了窗口大惺睢（10000）壳嚎，未設(shè)置步長(zhǎng)。
vcftools --vcf all.vcf --weir-fst-pop DA.txt?--weir-fst-pop non_DA.txt?--fst-window-size 10000 --out da_nonda
生成da_nonda.windowed.weir.fst 文件
CHROM BIN_START BIN_END N_VARIANTS WEIGHTED_FST MEAN_FST
1? ? ? 1? ? ? 10000? 7? ? ? 0.0976319? ? ? 0.0971726
1? ? ? 10001? 20000? 61? ? ? 0.0511299? ? ? 0.0499982
1? ? ? 20001? 30000? 49? ? ? 0.0372642? ? ? 0.0402562
1? ? ? 30001? 40000? 49? ? ? 0.0679565? ? ? 0.073192
1? ? ? 40001? 50000? 17? ? ? 0.0299695? ? ? 0.0355619
1? ? ? 50001? 60000? 56? ? ? 0.101204? ? ? ? 0.105008
1? ? ? 60001? 70000? 180? ? 0.0951801? ? ? 0.100651
1? ? ? 70001? 80000? 21? ? ? 0.07005 0.0728101
1? ? ? 80001? 90000? 60? ? ? 0.170809? ? ? ? 0.178863
1? ? ? 90001? 100000? 21? ? ? 0.0920029? ? ? 0.10329
vcftools --vcf DA.vcf --TajimaD 10000 --out DA
生成DA.Tajima.D文件
CHROM BIN_START N_SNPS TajimaD
1? ? ? 0? ? ? 9? ? ? 2.88793
1? ? ? 10000? 40? ? ? 3.12075
1? ? ? 20000? 70? ? ? 2.90904
1? ? ? 30000? 90? ? ? 3.37423
1? ? ? 40000? 49? ? ? 3.24297
1? ? ? 50000? 88? ? ? 3.31503
1? ? ? 60000? 329? ? 3.52148
1? ? ? 70000? 26? ? ? 2.38513
1? ? ? 80000? 82? ? ? 2.30261
1? ? ? 90000? 82? ? ? 2.37531
1? ? ? 100000? 12? ? ? 2.48976
1? ? ? 110000? 86? ? ? 2.43792
1? ? ? 120000? 52? ? ? 2.95367
1? ? ? 130000? 57? ? ? 3.03159
1? ? ? 140000? 71? ? ? 3.06152
1? ? ? 150000? 142? ? 3.478
1? ? ? 160000? 0? ? ? nan
1? ? ? 170000? 2? ? ? 0.34569
1? ? ? 180000? 0? ? ? nan
3末早、對(duì)Fst和TajimaD值進(jìn)行篩選烟馅、排序
首先對(duì)生成Fst文件進(jìn)行排序，將文件按照第六列荐吉，也就是MEAN_FST這列值焙糟，從大到小進(jìn)行排序
sort -nr -k6?da_nonda.windowed.weir.fst?
1 12910001 12920000 1 2.47818e-17 2.47818e-17
11? ? ? 11750001? ? ? ? 11760000? ? ? ? 1? ? ? 2.47818e-17? ? 2.47818e-17
11? ? ? 3990001 4000000 2? ? ? 2.18381e-18? ? 2.29795e-18
10? ? ? 7260001 7270000 1? ? ? 2.10168e-17? ? 2.10168e-17
7? ? ? 16280001? ? ? ? 16290000? ? ? ? 2? ? ? 1.0842e-17? ? ? 1.23909e-17
1? ? ? 13520001? ? ? ? 13530000? ? ? ? 2? ? ? 1.0842e-17? ? ? 1.23909e-17
7? ? ? 23470001? ? ? ? 23480000? ? ? ? 1? ? ? 0.9375? 0.9375
8? ? ? 18430001? ? ? ? 18440000? ? ? ? 1? ? ? 0.931034? ? ? ? 0.931034
1? ? ? 8430001 8440000 1? ? ? 0.845361? ? ? ? 0.845361
9? ? ? 9550001 9560000 1? ? ? 0.833333? ? ? ? 0.833333
9? ? ? 5930001 5940000 1? ? ? 0.833333? ? ? ? 0.833333
9? ? ? 13450001? ? ? ? 13460000? ? ? ? 1? ? ? 0.833333? ? ? ? 0.833333
排序之后發(fā)現(xiàn)，由于有些值過(guò)小样屠，在使用科學(xué)計(jì)數(shù)法時(shí)排在了正常計(jì)數(shù)值的前面穿撮，為了去掉這個(gè)錯(cuò)誤缺脉，先把使用科學(xué)計(jì)數(shù)法的值，統(tǒng)一歸為0悦穿，再進(jìn)行排序
awk '{if($6~/e/)$6=0}1' da_nonda.windowed.weir.fst > da_nonda.window.clean.fst && sort -rn -k6 da_nonda.windowed.clean.fst? > da_nonda.windowed.sort.fst
對(duì)文件第2列做減1處理攻礼，方便后面與TajimaD匹配
awk '{$2 = $2-1;print$0}' sp1_da.window.sort.fst > sp1_da.sort.fst
統(tǒng)計(jì)da_nonda.sort.fst有多少行
wc -l?da_nonda.sort.fst
17824
取Fst值最大的前10%窗口，作為候選選擇區(qū)域
head -n 1782?da_nonda.sort.fst > da_nonda.top0.1.fst
對(duì)于TajimaD文件栗柒，先清除掉第4列為nan的行礁扮，再對(duì)其進(jìn)行排序，取前5%和后5%的值
sed -e '/nan/d' DA.Tajima.D > DA.clean.tajimaD && sort -nr -k4 DA.clean.tajimaD > DA.sort.tajimaD
wc -l DA.sort.tajimaD
19484
head -n 974?DA.sort.tajimaD > DA.top0.05.tajimaD
tail -n 974 DA.sort.tajimaD > DA.bot0.05.tajimaD
取Fst前10%區(qū)域和TajimaD前5%（后5%）區(qū)域的交集,并按照染色體(scaffolds)順序排序瞬沦，將得到的這些窗口作為受選擇區(qū)域
cat DA.top0.05.tajimaD da_nonda.top0.1.fst > DA.top && awk 'pass==1 { count[$1,$2]++ } pass==2 { if(count[$1,$2]>1) print }' pass=1 DA.top pass=2 DA.top > DA.top.site
wc -l DA.top.site
40
head -n 20 DA.top.site > DA.top.keep.site
sort -n -k1?DA.top.keep.site >?DA.top.sort.site
得到的文件如下
1 17540000 2 1.9458
1? ? ? 28520000? ? ? ? 3? ? ? 1.96716
2? ? ? 17140000? ? ? ? 5? ? ? 2.41094
2? ? ? 19520000? ? ? ? 2? ? ? 2.03336
2? ? ? 19620000? ? ? ? 2? ? ? 1.97256
2? ? ? 19630000? ? ? ? 2? ? ? 2.03336
2? ? ? 21790000? ? ? ? 2? ? ? 1.92526
2? ? ? 28710000? ? ? ? 2? ? ? 2.03336
3? ? ? 12170000? ? ? ? 12? ? ? 2.63244
3? ? ? 14190000? ? ? ? 12? ? ? 2.74771
3? ? ? 25930000? ? ? ? 4? ? ? 1.94295
3? ? ? 4700000 4? ? ? 2.25325
4? ? ? 8030000 5? ? ? 2.04905
5? ? ? 8680000 16? ? ? 2.6894
6? ? ? 13390000? ? ? ? 7? ? ? 2.06337
6? ? ? 22180000? ? ? ? 3? ? ? 2.0041
7? ? ? 23620000? ? ? ? 2? ? ? 1.92526
8? ? ? 8650000 3? ? ? 2.24976
11? ? ? 5080000 3? ? ? 2.04485
12? ? ? 1160000 6? ? ? 2.23514
4太伊、使用snpEff對(duì)DA.vcf文件進(jìn)行注釋
snpEff的教程很多，先SnpEff 配置基因組注釋文件逛钻，再snpEFF注釋vcf-筆記僚焦，最終得到包含注釋信息的DA.annotation.vcf
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT sample12 sample13 sample14 sample15 sample16 sample17 samp
le18? ? ? ? sample19? ? ? ? sample20? ? ? ? sample21
1? ? ? 3923? ? .? ? ? T? ? ? C? ? ? 98? ? ? PASS? ? BQB=0.169762;HOB=0.5;ICB=1;MQ=40;MQ0F=0.0322581;MQB=0.000633889;MQSB=0.975265;RPB=0.103031;SGB=-0.662043;VDB=0.07878
29;DP=198;DP4=60,33,30,28;AN=10;AC=5;ANN=C|intergenic_region|MODIFIER|CHR_START-VJ01Gene00001|CHR_START-VJ01Gene00001|intergenic_region|CHR_START-VJ01Gene00001|||n.3923T>C|
|||||? GT:AD:ADF:ADR:DP:PL:SP? ./.:.:.:.:.:.:. 0/1:12,9:8,6:4,3:21:55,0,255:0? 0/1:27,14:17,8:10,6:41:132,0,255:1? ? ? ./.:.:.:.:.:.:. 0/1:15,11:9,6:6,5:28:108,0,255:0
? ? ? ? 0/1:18,13:11,5:7,8:31:105,0,255:5? ? ? ./.:.:.:.:.:.:. 0/1:21,9:15,4:6,5:30:70,0,255:6 ./.:.:.:.:.:.:. ./.:.:.:.:.:.:.
5、從物種的基因ID到GO ID
得到分化區(qū)域內(nèi)SNP的注釋信息曙痘，并去掉intron區(qū)和同義突變的位點(diǎn)
awk '$3 = "vcftools --vcf DA.annotation.vcf --chr"" " $1" " "--from-bp"" "$2" " "--to-bp"" "$2+10000" " "--recode --recode-INFO-all --out"" "$1"_"$2 {print $3}' DA.top.sort.site > vcf.sh
sh vcf.sh
grep -v "int" *.recode.vcf|grep -v "synonymous_variant" > non_int_synon.vcf
提取SNP位點(diǎn)所對(duì)應(yīng)的基因ID信息
awk '{print $8}' non_int_synon.vcf|grep -o 'VJ01Gene[0-9][0-9][0-9][0-9][0-9]' > geneID.txt
去重復(fù)
sort -n geneID.txt | uniq > clean_geneID.txt
將物種的基因ID與GO數(shù)據(jù)庫(kù)的ID對(duì)應(yīng)
awk 'ARGIND==1{a[$1]=$0} ARGIND==2 && ($1 in a) {print $0}' clean_geneID.txt annot.go.tab > geneID_goID.txt
提取GO ID
awk '{print $2}' geneID_goID.txt > goID.txt
less goID.txt
GO:0004725
GO:0006570
GO:0035335
GO:0004190
GO:0006508
GO:0004190
GO:0005764
GO:0006508
這樣就得到了Fst top0.1和TajimaD top 0.05兩者都存在的GO ID芳悲，即受到平衡選擇的基因，TajimaD bot 0.05篩選的是受到定向選擇的基因边坤，做法也是一樣的名扛。
最后，對(duì)我們得到的基因進(jìn)行GO富集分析茧痒，在線GO富集的工具有很多肮韧，基迪奧在線云分析平臺(tái)和遺傳所王秀杰課題組開(kāi)發(fā)的GOEAST平臺(tái)都很方便做⊥基迪奧平臺(tái)的富集分析結(jié)果如下
寫(xiě)在最后：從去年11月拿到下機(jī)的重測(cè)序數(shù)據(jù)惹苗，并簡(jiǎn)單搭建了服務(wù)器分析環(huán)境，一路摸著石頭過(guò)河耸峭，戰(zhàn)戰(zhàn)兢兢，在度娘和github兩位老師的指導(dǎo)下才逐步入門淋纲，希望和更多做群體重測(cè)序的小伙伴交流袄湍帧！QQ：657231499
?
