群體進(jìn)化與選擇信號：
●生活在世界不同區(qū)域的生物群體在歷史長河中經(jīng)歷千萬年的自然選擇秽褒、人工馴化、遷徙饱苟、遺傳漂變等事件的洗禮炊邦，有的物種滅絕编矾，有的物種分化出不同的亞種或亞群,甚至逐漸演變成新的物種。群體進(jìn)化研究就是用來追溯和揭露這個(gè)進(jìn)化過程的馁害。
●現(xiàn)代群體進(jìn)化研究窄俏，開始基于全基因組的變異信息，分析群體的遺傳背景碘菜、群體結(jié)構(gòu)凹蜈、遺傳多樣性、基因交流情況忍啸、物種的形成機(jī)制以及群體的進(jìn)化方向等生物學(xué)問題仰坦。
●群體進(jìn)化主要的應(yīng)用方向有:
1.人工馴化機(jī)制研究:通過對野生型和馴化型群體的研究，發(fā)現(xiàn)人工馴化作用下受選擇的區(qū)域计雌，進(jìn)而鑒定到與重要經(jīng)濟(jì)性狀相關(guān)的基因悄晃。
2.自然選擇機(jī)制分析:研究生存于不同地理環(huán)境下的群體，挖掘適應(yīng)性進(jìn)化過程中受到選擇的區(qū)域凿滤，進(jìn)而為育種提供環(huán)境適應(yīng)性基因資源妈橄。
3.種群歷史研究:通過分析物種的起源地及各個(gè)分布區(qū)域中群體的遺傳變異信息，探索物種的進(jìn)化過程翁脆。
●選擇信號分析是群體進(jìn)化這個(gè)大的研究領(lǐng)域中的一個(gè)研究方向/方法眷细，多用于前兩個(gè)應(yīng)用方向的研究。
●由于人工選擇或自然選擇作用造成的基因組結(jié)構(gòu)特征的變化稱為選擇信號鹃祖。選擇信號分析是生物信息分析中的“考古學(xué)”溪椎。它可以幫助我們尋找生物群體在千萬年的進(jìn)化過程中，在基因組上留下的痕跡恬口。
●這些信號與動(dòng)物的選育方向以及馴化適應(yīng)機(jī)制緊密相關(guān)校读，因此對選擇信號進(jìn)行檢測有助于挖掘與動(dòng)物經(jīng)濟(jì)性狀相關(guān)的基因，了解性狀形成的潛在的遺傳學(xué)基礎(chǔ)祖能，解析動(dòng)物的馴化適應(yīng)機(jī)制歉秫，解析群體間的差異機(jī)制。對于畜禽遺傳改良具有重要意義养铸。
選擇信號分析的特性和應(yīng)用場景
●對于物種的群體遺傳研究雁芙，本身就涉及到多個(gè)領(lǐng)域，如考古學(xué)钞螟、生態(tài)學(xué)兔甘、社會(huì)學(xué)，然后才是基因組學(xué)和生物信息學(xué)分析鳞滨。所以涉及群體遺傳的研究項(xiàng)目,需要盡量將這幾類因素綜合在一起討論洞焙，才能形成完整的故事。
1.有明確的研究方向和研究目的(如:研究豬的高原適應(yīng)性)
2.樣本材料具有相關(guān)特性且與研究目的相契合(如:藏豬具有高原適應(yīng)性)
3.對材料背景非常了解，包括起源地/棲息地澡匪、(地理位置..環(huán)境和生態(tài))熔任、 遷徙歷史、馴化歷史唁情、育成歷史等(四面環(huán)山的四川盆地成為冰川時(shí)期動(dòng)物群體的避難所)
4.當(dāng)定位到受選擇的區(qū)域之后疑苔，把受選擇區(qū)域的位點(diǎn)做基因注釋，需要了解與所研究性狀相關(guān)的基因甸鸟、蛋白夯巷、通路等(受選擇區(qū)域中發(fā)現(xiàn)'vitamin B6 binding相關(guān)基因:維生素B6可幫助血紅蛋白合成和促進(jìn)氧的綁定)
5.建議每個(gè)亞群個(gè)體數(shù)量不小于10個(gè)，且來自不同家系
遺傳背景分析
一哀墓、主成分分析(Principal Component Analysis)
●PCA是一種線性代數(shù)中的數(shù)據(jù)處理方法趁餐，它利用降維的思想，從高維度數(shù)據(jù)(如測序得到的百萬級別SNP位點(diǎn)數(shù)據(jù)) 中提取關(guān)鍵的信息篮绰，以便我們使用更少的變量(指標(biāo))就可以對樣本進(jìn)行有效區(qū)分后雷。這些被提取出的信息按照其效應(yīng)從大到小排列，我們稱之為主成分1(Principal Component1)吠各、主成分2臀突、主成分3...
●PCA分析的應(yīng)用場景:
1.檢測離群樣本
2.推斷群體分層和亞群間的遺傳距離
二、進(jìn)化樹
●又稱為系統(tǒng)發(fā)生樹贾漏，它利用樣本間的差異度將樣本進(jìn)行聚類候学，用一種類似樹狀分支的圖形來概括各物種之間的親緣關(guān)系，可用來描述物種之間的進(jìn)化關(guān)系和遺傳距離遠(yuǎn)近纵散。
●不同的構(gòu)樹方法
1.基于距離的方法:首先通過各個(gè)物種之間的比較梳码，根據(jù)一 定的假設(shè)(進(jìn)化距離模型)推導(dǎo)得出分類群之間的進(jìn)化距離，構(gòu)建一個(gè)進(jìn)化距離矩陣伍掀。 進(jìn)化樹的構(gòu)建則是基于這個(gè)矩陣中的進(jìn)化距離關(guān)系掰茶。(UPGMA, NJ)
2.基于特征的方法:不計(jì)算序列間的距離，而是將序列中有差異的位點(diǎn)作為單獨(dú)的特征蜜笤，并根據(jù)這些特征來建樹(ML MP)
●進(jìn)化樹的解讀
1.枝長:枝長累積距離越近的樣本差異越小
2.自展值:進(jìn)化樹分支可信度（藍(lán)圈濒蒋，百分比75%以上比較可信）
3.標(biāo)尺:代表序列的差異程度
三、群體結(jié)構(gòu)(structure) 分析
●先預(yù)設(shè)群體由若干亞群(k=n)構(gòu)成,通過模擬算法找出在k=n的情況下把兔，最合理的樣本分類方法沪伙。最后再根據(jù)每次模擬的最大似然值，找出最適用這群體的K值县好。
●應(yīng)用場景:
1.推算亞群劃分情況
2.推算群體基因交流程度
3.推算個(gè)體的血統(tǒng)構(gòu)成比例
●主流軟件
1.STRUCTURE
2.fastSTRUCTURE
3.Admixture
四围橡、LD衰減分析
●連鎖不平衡
●當(dāng)位于某一座位的特定等位基因與另一座位的某一等位基因同時(shí)出現(xiàn)的概率大于群體中隨機(jī)分布的兩個(gè)等位基因同時(shí)出現(xiàn)的概率時(shí)，就稱這兩個(gè)座位處于連鎖不平衡狀態(tài)。
●一般而言，兩個(gè)位點(diǎn)在基因組上離得越近，相關(guān)性就越強(qiáng)，LD系數(shù)就越大黔漂。反之，LD系數(shù)越小禀酱。也就是說炬守，隨著位點(diǎn)間的距離不斷增加，LD系數(shù)通常情況下會(huì)慢慢下降剂跟。這個(gè)規(guī)律减途，通常就會(huì)使用LD衰減圖來呈現(xiàn)。
●LD衰減圖就是利用曲線圖來呈現(xiàn)基因組上分子標(biāo)記間的平均LD系數(shù)隨著標(biāo)記間距離增加而降低的過程曹洽。
●大概的計(jì)算原理就是先統(tǒng)計(jì)基因組上兩兩標(biāo)記間的LD系數(shù)大小鳍置，再按照標(biāo)記間的距離對D系數(shù)進(jìn)行分類,最終可以計(jì)算出一定距離的分子標(biāo)記間的平均LD系數(shù)大小。
圖中橫軸為標(biāo)記間的物理距離送淆，縱軸為平均r方值税产。
群體多樣性越低，衰減速度越慢偷崩。
●LD衰減分析的應(yīng)用
1.評估群體特性和選擇強(qiáng)度:馴化選擇會(huì)導(dǎo)致群體遺傳多樣性下降辟拷，位點(diǎn)間的連鎖程度更高。所以阐斜，通常馴化程度越高選擇強(qiáng)度越大的群體,LD衰減速度也越慢衫冻。例如商品化群體比自然群體通常更大的LD衰減距離。類似的自然選擇谒出、遺傳漂變導(dǎo)致的群體遺傳多樣性下降隅俘，也會(huì)減慢LD衰減的速度。
2.檢測受選擇基因組區(qū)域:與有利突變緊密連鎖的中性位點(diǎn)會(huì)由于選擇作用在基因組上形成高頻率的核心單倍型笤喳，以其為中心向基因組兩側(cè)擴(kuò)展會(huì)形成長范圍的擴(kuò)展單倍型考赛。然而隨著與有利突變間距的增加，連鎖不平衡程度會(huì)相應(yīng)衰減莉测，在一定范圍內(nèi)各擴(kuò)展單倍型純合的總和占核心單倍型純合的比例可以被用來檢測基因組范圍內(nèi)的選擇作用颜骤。
3.GWAS分析中評估標(biāo)記密度是否足夠: GWAS分析本質(zhì)就是利用標(biāo)記和功能突變的相關(guān)性(LD關(guān)系)，來檢測與性狀相關(guān)的功能突變的位置捣卤。一般而言LD系數(shù)大于0.8就是強(qiáng)相關(guān)忍抽。如果LD系數(shù)小于0.1，則可以認(rèn)為沒有相關(guān)性董朝。如果LD衰減到0.1這么大的區(qū)間內(nèi)都沒有標(biāo)記覆蓋的話即使這個(gè)區(qū)間有一一個(gè)效應(yīng)很強(qiáng)的功能突變鸠项，也是檢測不到關(guān)聯(lián)信號的。所以通匙咏可以通過比較LD衰減(到0.1)距離和標(biāo)記間的平均距離來判斷標(biāo)記是否對全基因組有足夠的覆蓋度祟绊。(GWAS最低標(biāo)記量≈基因組大小/LD衰減距離)
選擇信號檢測方法
常用群體內(nèi)檢測指標(biāo)的計(jì)算方法大致分為三種：1.基于核苷酸多態(tài)性降低的π、θw；2.基于分離位點(diǎn)頻率的Tajima’D牧抽；3.基于連鎖不平衡增加的EHH嘉熊、iHS。以上三類指標(biāo)對應(yīng)于基因組受選擇特征的三個(gè)維度扬舒，而后才有了群體間的選擇指標(biāo)：1.由π衍生的π ratio阐肤、ROD、Fst讲坎；2.由EHH衍生的XPEHH孕惜。
https://zhuanlan.zhihu.com/p/52064863
對于單個(gè)物種，基于選擇的效應(yīng)晨炕，選擇信號檢測的方法可以被分為4大類:
1.基于等位基因頻率譜的方法
2.基于連鎖不平衡增加的方法
3.基于群體分化的方法
4.基于基因組雜合度的方法
一衫画、基于等位基因頻率譜
●基因型頻率和基因頻率的改變是選擇作用在基因組上最直接的體現(xiàn)∥屠酰基因頻譜(site-frequency spectrum)就是指某種等位基因在基因組上某個(gè)目標(biāo)區(qū)域內(nèi)出現(xiàn)的頻繁程度削罩。
●符合中性模型的群體，其群體中存在廣泛的遺傳多態(tài)遵馆，當(dāng)突變發(fā)生時(shí)總能夠維持在一個(gè)較低的頻率鲸郊，只有當(dāng)群體基因組上出現(xiàn)或存在有利突變時(shí)，選擇才會(huì)發(fā)生作用货邓，從而產(chǎn)生所謂的選擇清除或搭車效應(yīng)秆撮。
●代表性的檢測方法: Tajima's D, Fu andLi'sD, Fay and Wu'sH, CLR, Hp
●Tajima's D檢驗(yàn)的目的是區(qū)分隨機(jī)演變的DNA序列(“中性”)和在非隨機(jī)過程中演化的DNA序列，包括定向選擇或平衡選擇换况。
●Tajima's D的計(jì)算原理:多態(tài)位點(diǎn)數(shù)量和平均非匹配數(shù)量的差值职辨。
●D=0時(shí)，符合中性假設(shè)戈二，群體未受到選擇; D<0時(shí)舒裤，受到定向向選擇; D>0時(shí)，受到平衡選擇觉吭。
二腾供、基于連鎖不平衡
●基于連鎖不平衡理論，位點(diǎn)間的連鎖不平衡程度會(huì)隨標(biāo)記間距離的增加而逐漸降低鲜滩。因此伴鳖，在基因組上可以觀察到選擇作用造成的不同長度的擴(kuò)展單倍型純合(Extended Haplotype Homozygousity)。
●該方法的基本原理是:在中性條件下徙硅，基因組很難形成長范圍的連鎖不平衡的單倍型榜聂，因?yàn)樾峦蛔冃枰?jīng)歷漫長的遺傳漂變才能達(dá)到較高頻率，而在漫長的時(shí)間里會(huì)發(fā)生大量基因重組事件嗓蘑，使得這種連鎖不斷被打破须肆。而當(dāng)群體處于正向選擇作用下時(shí)匿乃，致因突變及其連鎖位點(diǎn)在正選擇的作用下，在短時(shí)間內(nèi)會(huì)達(dá)到較高頻率豌汇，形成大片段的純合單倍型幢炸。擴(kuò)展單倍型純合度檢驗(yàn)正是基于這樣的特征來篩選受選擇基因。
●代表性的檢測方法: EHH, XP-EHH, iHS, nSL, OmegaPlus
三瘤礁、基于群體分化
●同一物種不同群體之間由于環(huán)境不同或選擇目標(biāo)不同阳懂，其基因組等位基因頻率會(huì)表現(xiàn)出歧化選擇的效應(yīng)梅尤。這種現(xiàn)象在相同基因座位不同等位基因均受到選擇時(shí)表現(xiàn)尤為明顯柜思，即選擇加速群體分化。因此巷燥，基于群體分化的方法赡盘，不同群體同一等位基因頻率存在的差異程度大于兩個(gè)群體處于中性條件下的期望時(shí)，就推斷該位點(diǎn)存在選擇作用缰揪。
●代表性的檢測方法: Weir and Cockerhan's Fst, LSBL, di
●Fst的取值范圍為0-1陨享，1表示群體間完全分化的位點(diǎn)，0表示在群體間完全沒有分化的位點(diǎn)钝腺。
●基于Fst的的檢測方法多采用基因組單位點(diǎn)掃描的策略抛姑，而這樣的方式容易受到遺傳漂變等因素的影響，產(chǎn)生假陽性的顯著位點(diǎn)艳狐。為盡量減少假陽性的發(fā)生定硝，通常采用滑動(dòng)窗口的策略，降低這些干擾因素毫目，增加選擇信號檢測的準(zhǔn)確性蔬啡。
四、基于基因組雜合度
●當(dāng)基因組上特定區(qū)域受到選擇時(shí)镀虐，由于“選擇性清除”作用的存在箱蟆，該區(qū)域及其連鎖的區(qū)域表現(xiàn)為多態(tài)性降低，同時(shí)純和度增加刮便。因此對基因組的雜合度進(jìn)行檢測空猜，可以推斷出基因組中受到選擇的區(qū)域『藓担基因組上受選擇程度越高辈毯，則雜合度程度越低。
●代表性的檢測方法: θπRatio, ROH
●核苷酸多態(tài)性θπ比率越偏離1窖杀，受選擇程度越高漓摩。θπ比率的檢測公式如下:θπratio=θπA/θπB
其中，θπA和θπB分別代表A群體和B群體的θπ值入客。θπ比率大于1管毙， 反映A群體的基因組雜合度大于B群體的雜合度腿椎，則B群體相應(yīng)基因組區(qū)域受到選擇。θπ 比率小于1,則A群體的基因組雜合度低于B群體夭咬，則選擇發(fā)生在A群體對應(yīng)的基因組區(qū)域啃炸。
選擇信號的應(yīng)用策略
●由于單純一種選擇信號檢測方法容易造成假陽性選擇信號的產(chǎn)生，目前選擇信號檢測普遍采用兩種或多種檢測方法進(jìn)行組合策略卓舵，不同檢測方法相互驗(yàn)證南用，用以避免假陽性，例如基于群體分化和基于基因組雜合度檢測組合的策略掏湾、基于連鎖不平衡的檢測方法與基于基因組雜合度檢測方法組合的策略等裹虫。也有文獻(xiàn)使用3種以上的方法同時(shí)進(jìn)行選擇信號檢測，篩選2種或以上方法檢測重疊的位點(diǎn)融击。
●此外筑公，在功能基因挖掘方面，選擇信號檢測可以與其他功能基因挖掘策略進(jìn)行組合尊浪，提高基因定位的效率和準(zhǔn)確性匣屡，如選擇信號檢測與全基因組關(guān)聯(lián)分析的組合，選擇信號與拷貝數(shù)變異檢測的組合拇涤。
選擇信號的計(jì)算（腳本）
群體結(jié)構(gòu)分析樣本量差異基本都會(huì)有些影響捣作。
選擇信號分析的時(shí)候沒過濾，是因?yàn)樵谀壳暗膇bs分析中鹅士，fst的頻率抽樣券躁，受影響小些。群體結(jié)構(gòu)的時(shí)候比較擔(dān)心連鎖的影響如绸。選擇信號主要就是要找這些連鎖的地方嘱朽。
一、Fst（基于SNP位點(diǎn)的多態(tài)性的選擇信號的檢測）
Fst/pi都是基于SNP位點(diǎn)的多態(tài)性來檢測潛在的選擇信號區(qū)域怔接。Fst/pi ratio基于pi值
按照點(diǎn)來計(jì)算：
# 使用vcftools對每一個(gè)SNP變異位點(diǎn)進(jìn)行計(jì)算
nohup /home/software/vcftools/src/cpp/vcftools --vcf test.vcf --weir-fst-pop 1_population.txt --weir-fst-pop 2_population.txt --out p_1_2-single &
按照窗口區(qū)域來計(jì)算（逐個(gè)位點(diǎn)計(jì)算FST時(shí)搪泳，可能會(huì)出現(xiàn)FST值很高的假陽性信號（中性選擇導(dǎo)致），所以考慮到搭載效應(yīng)同時(shí)計(jì)算了滑窗FST）
# 計(jì)算Fst （此步驟生成的XXX.windowed.weir.fst 可用于畫曼哈頓圖）
nohup /home/software/vcftools/src/cpp/vcftools --vcf test.vcf --weir-fst-pop 1_population.txt --weir-fst-pop 2_population.txt --out p_1_2_bin --fst-window-size 500000 --fst-window-step 50000 &
# test.vcf是SNP calling 過濾后生成的vcf文件扼脐；
# p_1_2_3 生成結(jié)果的prefix
# 1_population.txt是一個(gè)文件包含同一個(gè)群體中所有個(gè)體岸军，一般每行一個(gè)個(gè)體。個(gè)體名字要和vcf的名字對應(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）即可。
# 為Fst排序
nohup sort -k 5 -g XXX.windowed.weir.fst > XXX.windowed.weir.sorted.fst &
# 窗口計(jì)數(shù)
wc -l XXX.windowed.weir.sorted.fst # 假設(shè)為10000
# 取前百分之一（需過濾空集）
nohup tail -n 100 XXX.windowed.weir.sorted.fst > XXX.windowed.weir.sorted.1%.fst &
# 計(jì)算Fst后輸出表格的格式: 1.CHROM（染色體號） 2.BIN_START （開頭）3.BIN_END（結(jié)尾） 4.N_VARIANTS（SNP數(shù)） 5.WEIGHTED_FST（加權(quán)Fst） 6.MEAN_FST（平均Fst）
# SnpEff注釋
Fst（群體間分化指數(shù)）的值在0-1之間罚攀，分化指數(shù)越大党觅，兩個(gè)群體的差異越大雌澄。Fst值為0表示兩個(gè)群體是隨機(jī)交配的，基因型完全相似杯瞻。如果Fst值為1镐牺，表示兩個(gè)群體完全隔離。
當(dāng)Fst出現(xiàn)負(fù)值時(shí)魁莉，表示位點(diǎn)在群體里無分化睬涧，一般處理為0再做后續(xù)分析。
曼哈頓圖可視化
http://www.reibang.com/p/db932369b2e8
1.制作csv表格
計(jì)算Fst后輸出表格的格式為: 1.CHROM（染色體號） 2.BIN_START （開頭）3.BIN_END（結(jié)尾） 4.N_VARIANTS（SNP數(shù)） 5.WEIGHTED_FST（加權(quán)Fst） 6.MEAN_FST（平均Fst）
awk -F '\t' '{print $1"\t"$2"\t"$5}' XXX.windowed.weir.fst > XXX.windowed.weir-125.fst
# 制作表格共四列:加表頭1.SNP 2.CHR 3.BP（BIN_START） 4.P（WEIGHTED_FST加權(quán)Fst）
# 第一列為SNP的名字旗唁，第二列CHR為所在染色體畦浓，第三列BP為染色體上所在位置。要注意如果CHR中存在X逆皮，Y宅粥，需要轉(zhuǎn)化為數(shù)字如賦予23参袱，24等电谣，其中第一列SNP的名字是可選擇的，后三列是必須提供的
# 以表格形式打開抹蚀，加上第一列SNP：snp1剿牺、snp2、snp3...
2.R可視化
install.packages('qqman')
library('qqman')
mydata2<-read.table("XXX.windowed.weir-125.fst",head=T)
color_set <- c("#8DA0CB","#E78AC3","#A6D854","#FFD92F","#E5C494","#66C2A5","#FC8D62")
par(cex=0.9)
manhattan(mydata2,col=color_set,logp = F,ylab='Fst',genomewideline = 0.6) # 在表格中查找前1%的最小值环壤，畫出閾值線晒来，本例0.6
230210
library(CMplot)
data1 <- read.table("CLR畫圖_24711_2.03.csv",header = T,sep = ',')
color_set <- c("#801e91","#344fa8","#f7cb34","#a8a8aa","#ffe4c7")
CMplot(data1, plot.type="m", LOG10=F,
chr.den.col=NULL, col = color_set,
threshold = 2.03, threshold.col = "red", threshold.lwd= 2, threshold.lty =1,
amplify = FALSE, file.output=T, height=5, width = 10,
ylab = "CLR",
pch =".", cex =4, dpi = 600, file = "jpg",memo = "CLR")
3.圖解析
y坐標(biāo)：Fst值
x坐標(biāo)：染色體號
二、Tajima's D
這個(gè)是選擇相關(guān)的一個(gè)參數(shù)郑现，大于0代表群體觀測雜合度高于預(yù)期雜合度湃崩，稀有等位基因頻率降低（群體收縮或者平衡選擇），小于0說明群體觀測雜合位點(diǎn)少于預(yù)期值接箫，稀有等位基因頻率增加（群體擴(kuò)張或者低頻選擇）攒读。 也就是說，只有0是正常的辛友，其他都是選擇發(fā)生薄扁。
nohup /home/software/vcftools/src/cpp/vcftools --vcf XXX.vcf --TajimaD 500000 --out TajimaD &
三、π废累，核苷酸多樣性邓梅，越大說明核苷酸多樣性越高，越低說明兩個(gè)座位DNA序列差異越小邑滨。(在群體內(nèi)一段一段進(jìn)行比較日缨，則雜合度越低越受選擇)
https://blog.csdn.net/yangl7/article/details/109546077
π用來分析堿基多態(tài)性，多態(tài)性越低掖看，受選擇程度越高匣距。取值時(shí)與Fst相反诈铛，需要取數(shù)據(jù)的后1%。
# A墨礁，B兩個(gè)群幢竹，僅關(guān)注A群，使用 --keep
nohup /home/software/vcftools/src/cpp/vcftools --gzvcf /home/important_data/lzxhebing/cattle_328_chr1-30_resequencing_final_vcf/cattle-328-chr1-30-snp-indel.filter.pass.vcf.gz --keep bohai.txt --window-pi 500000 --window-pi-step 50000 --out bohai_pi &
# 生成兩個(gè)文件:XXX.windowed.pi恩静、XXX.log焕毫，將性染色體刪除，PI列降序排序
nohup sort -n -k 5r bohai_pi.windowed.pi > bohai_pi.windowed.sorted.pi &
wc -l bohai_pi.windowed.sorted.pi #49770行
nohup tail -n 500 bohai_pi.windowed.sorted.pi > bohai_pi.windowed.sorted.1%.pi &
# 使用bohai_pi.windowed.sorted.1%.pi文件通過R進(jìn)行繪圖
pi（θπ） 選擇消除分析圖可視化
library(ggplot2)
data<-read.table("allsample_500K.windowed.pi",header=T)
for (i in 1:29){
result_name = paste("chr",i,".png",sep="")
png(result_name,width=1500,height=300)
chrom = subset(data,CHROM==i)
xlab = paste("Chromosome",i,"(MB)",sep="")
p <- ggplot(chrom,aes(x=BIN_START/1000000,y=PI,group=Phenotype,color=Phenotype,shape=Phenotype)) + geom_line(size=0.5) + xlab(xlab)+ ylab("Pi") + theme_bw()
print(p)
dev.off()
}
https://www.plob.org/article/21645.html
四驶乾、iHS檢測群體內(nèi)的選擇信號(基于單體型haplotype的選擇信號的檢測)
（22.10.26 這里的數(shù)據(jù)僅僅經(jīng)過geno邑飒、maf過濾，不需要ld過濾级乐，注意重測序數(shù)據(jù)需要填充８硐獭！７缈啤）
在selective sweeps選擇過程中撒轮，有些強(qiáng)烈受到選擇的位點(diǎn)variants由于LD的因素會(huì)連帶著其附近的位點(diǎn)variants一起被保留，并且不會(huì)受到重組recombination的打斷贼穆。一些低重組區(qū)域的haplotypes的長度會(huì)高于那些高重組區(qū)域的haplotypes的長度题山。因此，對比同一genomic區(qū)域在不同群體中的haplotype的長度可以用來判斷是否受到選擇故痊。例如：在一個(gè)群體內(nèi)部顶瞳，如果某一個(gè)體強(qiáng)烈受到選擇，其haplotype的長度會(huì)遠(yuǎn)長于其它個(gè)體愕秫；同理慨菱，對于兩個(gè)群體之間的比較，某一群體受到選擇戴甩，則其基因組中的受選擇區(qū)域的haplotypes會(huì)比未受到選擇群體中的haplotypes更長符喝。
例如：使用selscan軟件計(jì)算了澳洲野犬的iHS，并通過常染色體上20 kb的滑動(dòng)窗口通過規(guī)范（在selscan的軟件中）對分?jǐn)?shù)進(jìn)行歸一化等恐。如果其中30%的站點(diǎn)的iHS絕對值高于閾值(或iHS絕對值的前1%)洲劣，我們將窗口確定為候選區(qū)域。（參考Genomic regions under selection in the feralization of the dingoes）
http://www.reibang.com/p/23ab344d66f7
1.數(shù)據(jù)文件準(zhǔn)備
# 提取vcf 壓縮
nohup /home/software/bcftools/bcftools view -S id.txt test.vcf.gz -Ov > sampleid.vcf &
nohup bgzip -c -f sampleid.vcf > sampleid.vcf.gz &
# 過濾课蔬，去多等位基因位點(diǎn)囱稽，最大缺失率（max-missing）<10%即檢出率>90%的SNP位點(diǎn)，最小等位基因頻率（MAF）>0.03（哈代溫伯格平衡檢驗(yàn)的P值大于10^-6 -hwe）二跋；--min-alleles 2 --max-alleles 2 去除多等位基因和等位基因缺失
nohup /home/software/vcftools/src/cpp/vcftools --vcf sampleid.vcf --remove-indels --max-missing 0.9 --maf 0.03 --min-alleles 2 --max-alleles 2 --recode --recode-INFO-all --out sampleid.miss.snp.recode.vcf &
# 去其他染色體战惊，由于參考基因組版本問題，這里存在很多位于隨機(jī)或未知染色體上的變異位點(diǎn)扎即，所以需要先過濾掉這些位點(diǎn)
nohup /home/software/vcftools/src/cpp/vcftools --vcf sampleid.miss.snp.recode.vcf --recode --recode-INFO-all --not-chr XXX --out test_noincludeXXX &
# 轉(zhuǎn)格式
nohup /home/software/plink --vcf test_noincludeXXX.miss.snp.recode.vcf --recode --allow-extra-chr --out test_noincludeXXX.miss.snp.recode &
# 換map文件位置名稱 吞获，空格分隔
awk '{print$1,$1":"$4,$3,$4}' test_noincludeXXX.misssnp.recode.map > test_noincludeXXX.misssnp.recode.pos.map
# 轉(zhuǎn)換成vcf格式
nohup /home/software/plink --allow-extra-chr --file test_noincludeXXX.misssnp.recode.pos --make-bed --out test_noincludeXXX.misssnp.recode.pos &
# 更換染色體况凉，加頭，壓縮
2.vcf文件的phasing各拷，beagle轉(zhuǎn)換格式
由于后面計(jì)算時(shí)selscan軟件對vcf文件格式要求：前九列同標(biāo)準(zhǔn)的vcf格式一樣刁绒，后面緊接著的個(gè)體的基因型盡量是phasing后的。文件格式如下：< quality ><individual 1 genotype > ... < individual N genotype >烤黍，因此首先使用beagle軟件轉(zhuǎn)換格式
1 100 rs1 0 1 . . . . 0|1 0|0
1 200 rs2 0 1 . . . . 1|1 1|0
1 300 rs3 0 1 . . . . 0|0 1|0
1 400 rs4 0 1 . . . . 0|0 1|0
1 500 rs5 0 1 . . . . 1|0 0|1
nohup perl -pe "s/\s\.:/\t.\/.:/g" test.miss.snp.recode.vcf | bgzip -c > out.vcf.gz &（重測序數(shù)據(jù)會(huì)遇到這種問題）
nohup java -jar /home/software/beagle.25Nov19.28d.jar gt=out.vcf.gz out=test.beagle ne=44 > test.beagle.vcf.log &
#有報(bào)錯(cuò)：java.lang.OutOfMemoryError: Java heap space 解決辦法：
nohup java -jar -Xmn12G -Xms24G -Xmx48G /home/software/beagle.25Nov19.28d.jar gt=test_noincludeXXX.misssnp.recode.pos-chr-36head.vcf out=test.beagle ne=44 > test.beagle.vcf.log &
3.準(zhǔn)備分染色體的vcf文件知市，單倍型以及基于單倍型的分析都需要逐條染色體運(yùn)行，所以需要把vcf文件以及map按照染色體ID進(jìn)行拆分速蕊。
vi vcftools.sh
for i in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30;
do /home/software/vcftools/src/cpp/vcftools --vcf test.beagle.vcf --chr ${i} --recode --recode-INFO-all --out ${i}.test;
done
nohup bash vcftools.sh & # 得到${i}.test.recode.vcf文件
4.準(zhǔn)備分染色體的map文件
vi map.bash
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30;
do /home/software/vcftools/src/cpp/vcftools --vcf ${k}.test.recode.vcf --plink --out ${k}.MT;
done
nohup bash map.bash & # 得到${k}.MT.map文件
5.map計(jì)算遺傳距離
vi map2.sh
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 ;
do awk 'BEGIN{OFS=" "} {print 1,".",$4/1000000,$4}' ${k}.MT.map > ${k}.MT.map2;
done
nohup bash map2.sh & # 得到${k}.MT.map2文件
6. 計(jì)算iHS
vi iHS.sh
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30;
do /home/software/selscan/bin/linux/selscan --ihs --vcf ${k}.test.recode.vcf --map ${k}.MT.map2 --out ${k}.iHS;
done
nohup bash iHS.sh & # 得到 ${k}.iHS.ihs.out文件嫂丙，運(yùn)行時(shí)間久！
（1）# 重測序數(shù)據(jù)直接計(jì)算iHS（對每個(gè)100K窗口單倍型進(jìn)行打分规哲，標(biāo)準(zhǔn)化使其正態(tài)分布跟啤，并統(tǒng)計(jì)絕對值>2的SNP所占的比例）
# 提取結(jié)果
vi hu.sh
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 ;
do awk '{print '${k}',$2,$3,$4,$5,$6}' ${k}.iHS.ihs.out > ${k}.hu.ihs.out ;
sed -i 's/ /\t/g' ${k}.hu.ihs.out;
done
nohup bash hu.sh &
# 標(biāo)準(zhǔn)化（將norm文件第7列數(shù)據(jù)即iHS數(shù)據(jù)進(jìn)行統(tǒng)計(jì)，統(tǒng)計(jì)絕對值>2的SNP所占的比例）
vi hu2.sh
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 ;
do /home/software/selscan/bin/linux/norm --ihs --files ${k}.hu.ihs.out --bp-win --winsize 100000 ;
done
nohup bash hu2.sh &
# 提取結(jié)果
vi hu3.sh
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 ;
do awk '{print '${k}',$1,$2,$4}' ${k}.hu.ihs.out.100bins.norm.100kb.windows > ${k}.hu3.ihs.out.100kb.windows;
done
nohup bash hu3.sh &
# 得到兩個(gè)文件 XXX.hu.ihs.out.100bins.norm.100kb.windows與XXX.hu.ihs.out.100bins.norm
cat ./*.hu3.ihs.out.100kb.windows > all.hu3.ihs.out.100kb.windows
sort -k 4n,4 all.hu3.ihs.out.100kb.windows > all.hu3.ihs.out.100kb.windows.sort
# 按照第四列取前5%的位點(diǎn)唉锌，畫圖（此處也可使用XXX.hu.ihs.out.100bins.norm文件畫圖隅肥，但數(shù)據(jù)量過大，推薦使用XXX.hu.ihs.out.100bins.norm.100kb.windows）
制作表格糊秆，格式如下武福，逗號分隔保存：
X SNP CHR BP P
1 snp1 1 2400001 0.265873
2 snp2 1 2500001 0.449057
3 snp3 1 2600001 0.596059
# R(220907此處重測序數(shù)據(jù)應(yīng)該使用-log的數(shù)值畫圖，比較美觀)
install.packages('qqman')
library('qqman')
data1 <- read.table("all.iHS.csv",header = T,sep = ',')
color_set <- c("#8DA0CB","#E78AC3","#A6D854","#FFD92F","#E5C494","#66C2A5","#FC8D62")
par(cex=0.9)
manhattan(data1,logp = F,col=color_set)
# SNPEFF注釋痘番，參考下面
（2）以下芯片數(shù)據(jù)操作
# 提取結(jié)果 ( 可更改為 awk '{print '${k}',$2,$6}' )
vi awk2.sh
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30;
do awk '{print '${k}',$2,$6}' ${k}.iHS.ihs.out > ${k}._.ihs.out ;
done
nohup bash awk2.sh & # 得到{k}._.ihs.out
cat ./*._.ihs.out > all.test.ihs.out # 合并所有染色體結(jié)果
# 根據(jù)第四列（iHS值）從大到小排列：
sort -k 1n,1 -k 4n,4 all.test.ihs.out > all.test.ihs.sort.out
7.制作窗口文件（芯片數(shù)據(jù)）
設(shè)置窗口滑動(dòng)區(qū)間，計(jì)算落到各個(gè)區(qū)間的iHS總和平痰，計(jì)算平均值汞舱，提取iHS絕對值為top 1%的窗口區(qū)間，認(rèn)定為選擇信號強(qiáng)的區(qū)間宗雇，根據(jù)文獻(xiàn)昂芜，設(shè)置500kb的窗口；計(jì)算各個(gè)窗口內(nèi)的均值赔蒲，計(jì)算期望和標(biāo)準(zhǔn)差泌神，利用以下公式得到標(biāo)準(zhǔn)化的iHS：
# 使用excel表格中的函數(shù)制作50k窗口位置文件（每個(gè)染色體做一個(gè)），格式如下:
1 0 50000
1 50001 100001
1 100002 200002
...
# 由于窗口文件是從表格中粘貼的舞虱，需要注意分隔符的設(shè)置欢际，將空格鍵替換為tab鍵，并刪除最末尾的tab鍵矾兜；將每一個(gè)文件中的末尾為0损趋，即未在區(qū)間內(nèi)的行刪除，更改名稱為txt格式
sed -i 's/ /\t/g' *.window.bed
（（vi sed.sh
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30;
do sed -i 's/ /\t/g' ${k}.window.bed;
sed -i 's/[ \t]*$//g' ${k}.window.bed;
cat ${k}.window.bed | tr -s [:space:] > ${k}._window.bed椅寺；
done
得到{k}._window.bed））
8.制作區(qū)間文件（芯片數(shù)據(jù)）
# 制作位置區(qū)間bed文件
vi ihs.bed.sh
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30;
do awk '{print $1,$2,$2+1,$3}' ${k}._.ihs.out > ${k}._.ihs.bed;
done
nohup bash ihs.bed.sh &
# 取交集,需要統(tǒng)計(jì)各個(gè)SNP落在哪個(gè)窗口浑槽，用bedtools統(tǒng)計(jì)兩個(gè)區(qū)間（SNP位置蒋失，SNP+1為位置區(qū)間；窗口為第二個(gè)區(qū)間）的交集桐玻。
vi intersect.sh
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27;
do /home/software/bedtools2/bin/bedtools intersect -a ${k}.window.bed -b ${k}._.ihs.bed -wao > ${k}.window.ihs.intersect.bed;
done
nohup bash intersect.sh &
然后將每一個(gè)文件中的末尾為0篙挽，即未在區(qū)間內(nèi)的行刪除，更改名稱為txt格式
9.用ihs_window_mean_xunhuan.py腳本計(jì)算各個(gè)窗口（即第2镊靴，3列區(qū)間）內(nèi)iHS的均值（芯片數(shù)據(jù)）
vi ihs_window_mean_xunhuan.py
# -*- coding: utf-8 -*-
"""
Created on Sun Jul 11 22:16:50 2021
@author: Administrator
"""
import numpy as np
import glob
dirname ="/*.txt" #### 注意此處嚴(yán)格按照./*.txt格式
files = glob.glob(dirname)
for filename in files:
#刪除不需要的數(shù)據(jù)
data = []
fp=open(filename)
lines = fp.readlines()
for line in lines:
if line[-2] != '0':
data.append(line)
fp.close()
fp=open(filename, 'w')
for row in data:
fp.write(row)
fp.close()
#計(jì)算平均值
data = np.loadtxt(filename, dtype=float, delimiter='\t')
result = np.empty((0,8))
i = 0
count = 0
temp = data[0]
for row in data:
if temp[1] == row[1]:
count += row[6]
i += 1
else:
temp[6] = count / i
result = np.concatenate((result, np.expand_dims(temp,0)), axis=0)
temp = row
i = 1
count = row[6]
temp[6] = count / i
result = np.concatenate((result, np.expand_dims(temp,0)), axis=0)
filename = filename[:-3] + 'csv'
np.savetxt(filename, result, delimiter=',')
nohup python ihs_window_mean_xunhuan.py & （更改路徑嫉髓，注意此處嚴(yán)格按照./*.txt格式）
cat ./*.window.ihs.intersect.txt > all.window.ihs.intersect.txt
10.將all.window.ihs.intersect.txt放于電腦中按照染色體與區(qū)間位置進(jìn)行排序保存，制成all.window.ihs.intersect.xlsx邑闲，并將.window.ihs.intersect.csv放于電腦中算行，取出各個(gè)表格（.window.ihs.intersect.csv）的第七列，與all.window.ihs.intersect.xlsx對應(yīng)上（即對區(qū)間進(jìn)行去重復(fù)整理一下苫耸，單做出個(gè)工作表州邢，好久啊）褪子，將所有的按順序排好的這列制作成all.window.ihs.intersect.csv量淌，在R中計(jì)算標(biāo)準(zhǔn)iHS值，得到window_biaozhun_ihs.csv文件嫌褪，取第二列數(shù)據(jù)呀枢，即標(biāo)準(zhǔn)化的iHS值(有正有負(fù))。
setwd("")
sink("window_biaozhun_ihs.csv")
data=read.table("all_window_ihs.csv")
options(max.print=40000) ###表格有多少行就設(shè)置大小以保證所有值都能讀取到
scale(data,center=TRUE,scale=TRUE)
attr(1,"scaled:center")
attr(1,"scaled:scale")
sink()
11.然后將表格window_biaozhun_ihs3.sort_ihs.csv與all.window.ihs.intersect.xlsx中的染色體笼痛、區(qū)間位置對應(yīng)上裙秋，將這列標(biāo)準(zhǔn)化的iHS值，去掉負(fù)號取絕對值缨伊，其中|iHS|>1.96對應(yīng)的區(qū)間則為顯著性選擇區(qū)間摘刑，按染色體和位置信息排序隨后畫曼哈頓圖，其中刻坊，格式包含四列（SNP CHR BP P,分別代表snp名稱枷恕，染色體號，snp位置谭胚，相應(yīng)的fst值或其他要用曼哈頓圖展現(xiàn)的值徐块，這里為iHS值），注意此處的畫圖為Snp位點(diǎn)的位置灾而，因此需要重新對應(yīng)一下胡控，新建一個(gè)表格manhadun.csv。
install.packages('qqman')
library('qqman')
data1 <- read.table("manhadun.csv",header = T,sep = ',')
color_set <- c("#04047c","#cb8404")
par(cex=0.9)
manhattan(data1,col=color_set, logp = F,ylim=c(0,7) ,ylab='fst',suggestiveline = -log10(1e-08) )
abline(lty = 1 , h = 1.96, col = "red",lwd = 2)
##### suggestiveline 設(shè)置"suggestive"線的閾值绰疤，默認(rèn)為-log10(1e-5)铜犬，因此如果不想要這個(gè)線，可以放大這個(gè)值。
五癣猾、CLR檢測群體內(nèi)的選擇信號（基于似然法）
【群體內(nèi)選擇信號】SweeD--CLR選擇 - 知乎 (zhihu.com)
群體遺傳|基于自然選擇的CLR分析 (qq.com)
CLR敛劝，composite likelihood ratio，意即復(fù)合似然比檢驗(yàn)纷宇，主要用來探測群體內(nèi)的選擇信號夸盟，探測的原理是SFS。做群體內(nèi)選擇其實(shí)就算是一個(gè)退而求其次的做法像捶，計(jì)算簡單上陕，信號也容易找到并整合，故而可以在做分化較大的群體或者物種間的選擇信號做考慮拓春。
$ 第一步释簿，拆分群體與染色體，故而我們需要人為設(shè)置好選擇的窗口大小硼莽，即每個(gè)染色體各自劃分為多大的網(wǎng)格較為合適庶溶，所以，第一步懂鸵，應(yīng)該將vcf文件按染色體拆分偏螺，而后統(tǒng)計(jì)各自的長度，設(shè)置窗口大小匆光。
$ cd /public/home/ysq/220228_hucattleCY_hebing/20220702hu/306/CLR
$ source /public/home/sh/software/Miniconda3/bin/activate /public/home/sh/software/Miniconda3/envs/mytools
$ nohup bash clr.sh &
for i in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 ;
do vcftools --vcf sample_CY_id.vcf --recode --recode-INFO-all --stdout --chr ${i} --keep CY_id.txt > CY_id.chr${i}.vcf ;
done
$ 按照100kb窗口進(jìn)行劃分
$ vi chr_gird.txt
1,1580
2,1360
3,1200
4,1190
5,1200
6,1170
7,1100
8,1130
9,1040
套像。。终息。夺巩。。采幌。
露露安裝sweed
$ cd /home/hmy/210329test/220713CY/CLR
（實(shí)驗(yàn)一下單獨(dú)做一下每條染色體腳本
$ nohup /home/ywf/software/sweed/SweeD -name CY_id.chr1.100kb -input CY_id.chr1.vcf -grid 1580 -maf 0.05 -missing 0.1 &
$ nohup /home/ywf/software/sweed/SweeD -name CY_id.chr1.50kb -input CY_id.chr1.vcf -grid 3000 -maf 0.05 -missing 0.1 &）
循環(huán)腳本
for line in `cat chr_grid.txt`
do
chr=${line%%,*}
grid=${line##*,}
/home/ywf/software/sweed/SweeD -name ${chr}.100kb \
-input QC.JBC.filter-geno005-maf005.chr${chr}.vcf \
-grid ${grid} \
-minsnps 200 \
-maf 0.05 \
-missing 0.1;
done
$ R中進(jìn)行可視化（使用原來的就行了劲够，這個(gè)參考）
setwd("/public3/user/houhzw/test/sweed_rlt/xm")
chr2 <- read.table(file="SweeD_Report.pop1.chr2.10kb",sep = "\t", header = T)
plot(chr2[,1]/1000,chr2[,2])
去除文件前3行
$ sed -i '1,3d' SweeD_Report.*.100kb
添加染色體號
$ for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do awk '{print '${k}',$1,$2,$3,$4,$5}' SweeD_Report.${k}.100kb > SweeD_Report.chr.${k}.100kb;
sed -i 's/ /\t/g' SweeD_Report.chr.${k}.100kb;
done
$ cat ./*.txt > all.txt 對結(jié)果文件進(jìn)行處理，取前5% （49480*0.05=2474）
#其中休傍，第二列為位置，第三列為似然法計(jì)算出的值蹲姐，對第三列進(jìn)行排序磨取，表格格式為：
3 28381127.59 3.27E+01 1.13E-06 17786962 38963961
3 28330729.42 3.08E+01 1.14E-06 17897325 38681860
3 28431525.77 3.05E+01 1.14E-06 17990096 38681860
3 28481923.94 2.51E+01 1.20E-06 18532592 38444354
3 28280331.25 2.12E+01 1.27E-06 18822027 37747588
3 28532322.11 1.69E+01 1.36E-06 19707239 37329898
14 570151.1094 1.58E+01 7.10E-06 68024 2187293
$ 對位置進(jìn)行處理，取整數(shù) =ROUND(A1,0) 柴墩，替換染色體號忙厌，注釋
SnpEff注釋
http://www.reibang.com/p/b2b45d2523db
https://www.cnblogs.com/zhanmaomao/p/10964636.html
1.安裝 ，構(gòu)建數(shù)據(jù)庫
#會(huì)產(chǎn)生一個(gè)snpEff目錄 所有的程序都在這里面
下載 https://pcingola.github.io/SnpEff/
unzip snpEff_latest_core.zip
2.配置自己的基因組和注釋文件
以綿羊（sheep）參考基因組為例：打開snpEFF文件夾下的snpEff.contig江咳，在Third party databases下面增加新的物種信息：
#-------------------------------------------------------------------------------
# Third party databases （原來的）
#-------------------------------------------------------------------------------
#sheep genome, version sheepv1.0 (新加的)
sheepv1.0.genome : sheep （新加的）
# Databases for human GRCh38 (hg38) （原來的）
3.注釋文件為 gff 格式
# 在SnpEff目錄下逢净，創(chuàng)建目錄 data, data/sheepv1.0, data/genomes
mkdir data && cd data
mkdir sheepv1.0
mkdir genomes
# 將 gff3 (.gff) 文件放入sheepv1.0目錄下，并改名為 genes.gff
# 將基因組序列文件（.fasta）放入 genomes 目錄下，并改名為 sheepv1.0.fa
nohup java -jar snpEff.jar build -gff3 -v sheepv1.0 -noCheckCds -noCheckProtein & # 在 snpEff 目錄下爹土，執(zhí)行命令甥雕，若不加此參數(shù)的話-noCheckCds -noCheckProtein，需提供cds與蛋白文件Ｕ鸵稹Ｉ缏丁（220926）
4.注釋文件為 gtf 格式
# 如果注釋文件為.gtf，可參考 gff 中的步驟琼娘，要將注釋文件重新命名為 genes.gtf
nohup java -jar snpEff.jar build -gtf22 -v sheepv1.0 &
5.開始注釋 (提取的要注釋的文件中染色體應(yīng)為NC格式且文件以tab鍵分隔峭弟，提取三列即可：染色體號、起始位置脱拼、終止位置瞒瘸，若為SNP位點(diǎn)則第三列用1代替，且位于snpEff目錄下操作命令）
（1）若文件為按照窗口計(jì)算Fst后輸出的文件熄浓，則提取為bed文件格式情臭，bed格式(取染色體號，起始位置和結(jié)束位置和Fst 的值)
# 計(jì)算Fst后輸出表格的格式: 1.CHROM（染色體號） 2.BIN_START （開頭）3.BIN_END（結(jié)尾） 4.N_VARIANTS（SNP數(shù)） 5.WEIGHTED_FST（加權(quán)Fst） 6.MEAN_FST（平均Fst）
awk -F '\t' '{print $1"\t"$2"\t"$3"\t"$5}' XXX.windowed.weir.sorted.1%.fst > XXX.windowed.weir.sorted.1%.bed
#如果bed文件中染色體為1.2.3格式玉组，需替換染色體
nohup java -Xmx4g -jar ./snpEff.jar sheepv1.0 -i bed XXX.windowed.weir.sorted.1%.bed-chr.bed > XXX.windowed.weir.sorted.1%.bed-chr.anno.bed &
（2）若文件為按照位點(diǎn)計(jì)算（例如：fst按位點(diǎn)計(jì)算谎柄、重測序SNP數(shù)據(jù)、重測序INDEL數(shù)據(jù)）但重測序使用中惯雳，得到的注釋文件不理想
# 用awk提取fst中的snp 位置信息朝巫，并用vcftools對應(yīng)回vcf文件中
nohup awk '{print $1,$4,$4+1}' test.map > test-awk.bed &
# 注釋
nohup java -Xmx4g -jar ./snpEff.jar sheepv1.0 -i bed test-awk.bed > test-awk.anno.bed &
（3）若文件為vcf文件 (重測序SNP數(shù)據(jù)、重測序INDEL數(shù)據(jù))
比如只關(guān)注CDS中的注釋信息石景，不考慮上游劈猿、下游、UTR潮孽、基因間區(qū)等信息
可以選擇以下參數(shù)簡化輸出
-no-downstream
-no-upstream
-no-utr
-no-intergenic
-no-intron
nohup java -jar snpEff.jar ann -no-utr -no-downstream -no-upstream -no-intergenic sheepv1.0 test.vcf.gz > snpeff.vcf &
# 此處可將變異信息刪除做一個(gè)小的vcf文件揪荣，只保留到vcf的info信息即可，則最后出現(xiàn)的文件小一些
6.結(jié)果說明
(1)SnpEff結(jié)果解讀 - 簡書 (jianshu.com)
第一列：CHROM 發(fā)生突變的染色體ID往史。
第二列：POS：發(fā)生突變的染色體上的具體位置仗颈。
第三列：ID 可以在后面的注釋信息中找到geneID。
第四列：REF 參考基因組上的堿基或者序列椎例。
第五列：ALT 發(fā)生突變后的堿基或者序列挨决。
第六列：QUAL 得分，Phred格式的數(shù)值订歪。代表著此為點(diǎn)是純和的概率脖祈。此值越大，概率越低刷晋，代表著此為點(diǎn)是變異位點(diǎn)的可能性越大盖高。
第七列：FILTER 過濾情況 ：一般分析后的結(jié)果都為PASS慎陵，則表示該位點(diǎn)是變異位點(diǎn)。
第八列：INFO 變異位點(diǎn)的相關(guān)信息喻奥。
第九列：FORMAT 變異位點(diǎn)的格式：比如 GT:PL:ADF:ADR:AD:GP:GQ
第十列：SAMPLEs 各個(gè)樣本的值席纽，這些值對應(yīng)著第9列的各個(gè)部分，不同部分之間的值使用冒號分隔映凳。
而SnpEff結(jié)果文件中胆筒，在INFO這一列中，增添了一個(gè)字段诈豌，ANN
ANN=A|upstream_gene_variant|MODIFIER|AT1G69210|AT1G69210|transcript|AT1G69210.1|protein_coding||c.-3686C>T|||||3686|,A|upstream_gene_variant|MODIFIER|AT1G69210|AT1G69210|transcript|AT1G69210.2|protein_coding||c.-3686C>T|||||3686|,A|downstream_gene_variant|MODIFIER|SP1L2|AT1G69230|transcript|AT1G69230.1|protein_coding||c.2799C>T|||||2509|,A|downstream_gene_variant|MODIFIER|MES15|AT1G69240|transcript|AT1G69240.1|protein_coding||c.4371C>T|||||4138|,A|downstream_gene_variant|MODIFIER|SP1L2|AT1G69230|transcript|AT1G69230.2|protein_coding||c.*2799C>T|||||2568|,A|intron_variant|MODIFIER|SIK1|AT1G69220|transcript|AT1G69220.1|protein_coding|8/17|c.1323+52C>T||||||,A|intron_variant|MODIFIER|SIK1|AT1G69220|transcript|AT1G69220.2|protein_coding|8/17|c.1242+52C>T||||||
新增的字段由|進(jìn)行間隔仆救，并且這個(gè)字段中包含了突變位點(diǎn)的注釋信息，因此非常重要矫渔。
重點(diǎn)關(guān)注以下幾點(diǎn)：
Allele ANN=A 說明：突變后的堿基是A彤蔽。
Annotation upstream_gene_variant 造成的基因上游的突變。downstream_gene_variant 造成的基因下游的突變庙洼。
Annotation_Impact 突變位點(diǎn)造成的影響：一般可以劃分為四類 HIGH顿痪，MODERATE，LOW油够，MODIFILER 一般突變位點(diǎn)造成HIGH最好蚁袭，往后依次效果越低。
Gene_Name 基因名稱石咬。
Gene_ID 基因ID揩悄。
Feature_Type 想要分析的特征類型白魂，transcript, motif, miRNA 等为迈。
Feature_ID 根據(jù)Feature Type指定的特征，給出對應(yīng)的ID油坝。
Transcript_BioType 轉(zhuǎn)錄本類型, 通常采用Ensembl數(shù)據(jù)庫的轉(zhuǎn)錄本類型焕窝。
(2)snpEff_genes.txt和snpEff_summary.html這兩個(gè)文件記錄總結(jié)性信息比較簡單蹬挺。
*.ann.vcf 是一個(gè)注釋結(jié)果文件，其就在vcf的info信息新添加了anno一列信息它掂，其具體每個(gè)值含義如下：
1.Allele：突變之后的堿基巴帮，第一個(gè)突變位點(diǎn)由T堿基突變成了C堿基，對應(yīng)Allel的值為C虐秋；
2.Annotation：由sequence ontology定義的突變類型晰韵；
3.Annotation_Impact：對變異位點(diǎn)有害程度的簡單評估，取值有HIGH, MODERATE, LOW, MODIFIER 4種熟妓；
4.Gene_Name：基因名字；
5.Gene_ID：基因ID栏尚；
6.Feature_Type:想要分析的特征類型起愈，transcript, motif, miRNA 等;
7.Feature_ID:根據(jù)Feature Type指定的特征，給出對應(yīng)的ID;
8.Transcript_BioType:轉(zhuǎn)錄本類型, 通常采用Ensembl數(shù)據(jù)庫的轉(zhuǎn)錄本類型;
9.Rank:只有當(dāng)變異位點(diǎn)位于基因區(qū)域時(shí)才有值，會(huì)給出變異位點(diǎn)所處的exon/intron的編號和該基因的exon/intron的總數(shù)抬虽，比如一個(gè)突變位點(diǎn)位于基因的第3個(gè)exon上官觅，該基因一共有12個(gè)exon, 對應(yīng)的Rank的值為3/12,當(dāng)變異位點(diǎn)位于基因區(qū)域以外時(shí)，該字段的值為空;
10.HGVS.c:采用HGVS標(biāo)準(zhǔn)命名的基因水平的變異情況;
11.HGVS.p:采用HGVS標(biāo)準(zhǔn)命名的蛋白質(zhì)水平的變異情況阐污，只有當(dāng)突變位點(diǎn)位于編碼區(qū)是才會(huì)有值;
12.cDNA.pos/cDNA.length:突變位點(diǎn)在cDNA上的位置/cDNA的總長度;
13.CDS.pos/CDS.length:突變位點(diǎn)在CDS上的位置/CDS的總長度;
14.AA.pos/AA.length:突變位點(diǎn)在氨基酸序列上的位置/氨基酸序列的總長度;
15.Distance:變異位點(diǎn)與最近的特征的距離休涤，當(dāng)變異位點(diǎn)位于基因間區(qū)時(shí)，會(huì)給出與最近的基因之間的距離笛辟；當(dāng)變異位點(diǎn)位于exon區(qū)域時(shí)功氨，會(huì)給出與最近的內(nèi)含子邊界的距離，不同的情況手幢，距離的定義不同;
16.ERRORS/WARNINGS/INFO:對注釋結(jié)果的可靠程度進(jìn)行評估捷凄。
其他注釋方法
一、ensembl數(shù)據(jù)庫BioMart注釋
可觀察到所需的目標(biāo)物種位點(diǎn)注釋信息
1.ensembl——BioMart
2.Attributes #按需選擇
3.Results
二围来、ensembl vep注釋
1.通過位點(diǎn)計(jì)算Fst跺涤，再將利用Fst文件對應(yīng)出的那1%的vcf文件的前5列，用excel打開復(fù)制前5列监透，（從染色體那一行開始復(fù)制）使用NC染色體
