生信分析16:基于PSMC估計有效群體大小

本次推送是文獻分享16的對應內容。

Fig 1

有效群體大小的估計或者種群歷史動態(tài)的估計經常出現(xiàn)在基因組/重測序相關文章中,可視化形式一般是Fig 1。

Fig 1的解讀:由于該分析是以個體來推測群體,所以作者隨機選取了4份圓果化香 (Platycarya longipes) 和4份化香樹 (Platycarya strobilacea) 的重測序數(shù)據,分別比對到各自的參考基因組上。然后基于每一份數(shù)據去估計該物種的歷史有效群體大小。圖中橫坐標是距離今天的時間,從左到右是距離今天104年-距離今天107年??v坐標是有效群體大小。圖中淺藍色線有四條,代表隨機選取的4份圓果化香的估計值,一條深藍色線代表4份圓果化香的估計值擬合成的值。黃色代表化香樹。陰影部分的豎線代表兩物種的估計分化時間。g是指該物種繁殖一代的時間,單位為年。μ是該物種的突變速率。

通過這個圖我們可以看出兩物種的分化時間,以及分化之后各自種群大小的變化,可以與各種地質事件結合。

相關概念

有效群體大小 (Effective population size, Ne)

(1)是指與實際群體具有相同基因頻率方差或相同雜合度衰減率的理想群體大小,它反映了群體平均近交系數(shù)增量的大小以及群體遺傳結構中基因的平均純合度。

(2)有效群體數(shù)量并不是根據統(tǒng)計群體數(shù)量來定義的,而是根據遺傳變異在群體中的表現(xiàn)來定義群體大小的。如果一個群體在遺傳漂變的作用下,等位基因頻率在一定時間內發(fā)生的改變很緩慢,那么我們就認為該群體的有效群體數(shù)量很大;反之,如果等位基因頻率變化很大,那么我們就認為該群體的有效群體數(shù)量很小。

(3)在一個理想群體中,在隨機遺傳漂變影響下,能夠產生相同的等位基因分布或者等量的同系繁殖的個體數(shù)量,通常小于絕對的種群大小。

以上是我整理的不同來源的定義方式,看起來還是有些模糊。我自己的理解是在一個群體中,可以交配產生后代的個體數(shù)。比如在一個10000人的群體中,有2000人不具有生殖能力,有1000人不具有生殖意愿,只有剩余的7000個可以“自由交配”產生后代。那我們根據后代的基因型只能估計出這7000個自由交配的人,而無法判斷出實際人群大小。(純屬個人的理解,可能有錯誤,歡迎討論)

物種的突變速率μ

這是估計Ne時必須用到的一個值。Fig 1中標記的突變速率為2.06×10-9,單位是每代每堿基。不同物種的突變速率是有差異的,因此不能套用同一個值。對于常見的物種擬南芥、水稻、玉米等,突變率一般是可以查到。但是大多數(shù)物種的μ是查不到的,可以根據公式μ = Ks/2T進行估計,其中Ks是兩物種間的synonymous divergence, T是兩物種間的divergence time。

PSMC原理介紹

Fig 2

估計Ne的最經典軟件是2011年李恒博士開發(fā)的PSMC軟件,盡管該軟件是針對動物模型,但是在植物中也有較多的應用,并且化香樹的文章也是用的PSMC,所以本次推送以PSMC為例進行講解。

PSMC的全稱是成對序列馬可夫共祖先分析 (Pairwise Sequentially Markovian coalescent),利用單個個體的重測序數(shù)據來推測該個體所屬的種群在歷史上各階段的有效群體大小。該流程的核心是計算最近共同祖先時間(the Time since the Most Recent Common Ancestor, TMRCA)。假設TMRCA屬于某個時期的片段比例越少,則該時期的Ne越大。

PSMC適用于二倍體物種,且該群體需要隨機交配。對于二倍體物種,基因組上是存在廣泛的純合和雜合位點的(Fig 2)。在最近共同祖先中,認為均是純合位點,隨著時間的推移,經歷了重組和突變積累了雜合位點。某區(qū)段的雜合位點越多說明經過的代數(shù)/時間越久,所以該區(qū)段最近共祖先離現(xiàn)在越遠。

基于以上原理,可以根據單個個體的重測序數(shù)據估計該個體所在群體的Ne.

需要注意的是PSMC只適用于2萬年-300萬年之間的Ne估計,更短或更長均會導致準確性降低,且重測序數(shù)據的深度不低于18×.

PSMC實操

Fig 3

軟件的安裝參考https://github.com/lh3/psmc

輸入文件(Fig 3):

1、重測序數(shù)據比對到參考基因組的bam文件

2、參考基因組

第一步?生成二倍體一致性序列(Fig 4)

Fig 4

需要注意兩個參數(shù):

-d 指定最小深度,作者建議設置成平均深度的三分之一。

-D 指定最大深度,作者建議是平均深度的兩倍。

Fq2psmcfa是psmc自帶的一個命令,用于將fq.gz轉換為fa格式。

第二步 運行psmc (Fig 5)

Fig 5

Fig 6

-N? 迭代的最大次數(shù)

-t? 2N0的最大世代時間,代表TMRCA的上限

-r? 起始θ/ρ率

以上三個參數(shù)對于大部分植物使用,可不用更改。

-p ?代表PSMC有效群體大小變化的圖隨時間變化劃分的時間階段單位數(shù)量,例如默認的-p “4+25*2+4+6”表示從古至今依次經歷了一個4單位+25個2單位+一個4單位+一個6單位的時期。數(shù)字越大需要的計算資源越多。可根據經驗調整,也可查找近緣物種的值用作參考。

(參考自生信技工)

第三步 畫圖(Fig 7)

Fig 7

Fig 8

PSMC內置了畫圖腳本psmc_plot.pl

其中兩個最關鍵的參數(shù)是-u指定該物種的突變率,-g指定該物種繁殖一代的時間,單位為年,比如人默認為25年,黃瓜、西瓜等認為是1年。

-p指定輸出pdf格式(Fig 8)。

Fig 9

當然除了PSMC軟件,近些年也開發(fā)出了更多的用來推斷種群歷史動態(tài)/有效群體大小的軟件,可能會在后續(xù)的文獻分享-生信分析中涉及。

最后,推薦一個網站(https://methodspopgen.com/methods-to-infer-populations-history/,F(xiàn)ig9),里面推薦了很多群體遺傳學軟件。

參考鏈接

[1] http://www.itdecent.cn/p/33ab4d2c81f6

[2]https://zhuanlan.zhihu.com/p/364652712

[3]https://baike.baidu.com/item/%E6%9C%89%E6%95%88%E7%BE%A4%E4%BD%93%E5%A4%A7%E5%B0%8F/5975113

[4]https://yanzhongsino.github.io/2022/10/17/bioinfo_psmc/

[5]https://cloud.tencent.com/developer/article/1911384

[6]Li, H., Durbin, R. Inference of human population history from individualwhole-genome sequences. Nature 475, 493–496 (2011).?

本文使用 文章同步助手 同步

最后編輯于
?著作權歸作者所有,轉載或內容合作請聯(lián)系作者
【社區(qū)內容提示】社區(qū)部分內容疑似由AI輔助生成,瀏覽時請結合常識與多方信息審慎甄別。
平臺聲明:文章內容(如有圖片或視頻亦包括在內)由作者上傳并發(fā)布,文章內容僅代表作者本人觀點,簡書系信息發(fā)布平臺,僅提供信息存儲服務。

相關閱讀更多精彩內容

友情鏈接更多精彩內容