[RAD seq classic] 關於RAD研究中的缺值。

相較於Sanger定序的資料,在RAD的資料中最常被別人詬病的不外乎幾點:缺值比例高,SNP來源不清等。這一篇就先來談談RAD研究中的缺值問題,我也挺有興趣的。

由於在RAD的研究中,我們會建構出一個Library(裡頭會有許多樣本個體,端看你所需要的coverage大小,來決定放的個體數量),一次使用次世代定序所得到的資料量是固定的 (譬如說 a HiSeq Illumina run may generate 140 million reads. )但是!即使每個放入library的樣本DNA濃度一樣,但最後分配到不同物種的資料量並非是平均,而是隨機的,如此一來缺值的情況就會發生。然而,為什麼呢?其實導致缺值的原因有很多,像是樣本的限制酶切位位置的突變,或是在後端資料處理時,像是我們是用minimal depth或變異量作為參數來決定是否要使用此loci等。但由於使用的資料量太大,我們通常很少深入一一探討。但現在使用RADseq的分子親緣關係或生物地理研究,面對缺值所造成得影響或限制,以及如何處理和解決,已成為必要議題。然而不同研究面對缺值所使用的方法不盡相同,有些研究會去除有缺值的loci (McCormack et al. 2012 ; Zellmer et al. 2012 [尚未讀過]),有些會納入使用,即使他的缺值比例甚高(像是慈鯛的例子就是,但我讀過大部分的paper都還是會納入,只是比例多少而已),但無論哪一個決定都會影響最後的結果。但目前不確定的是納入缺值所得到的結果,是否能得到精確的推論。然而2014年Huang和Knowles的模擬研究便是探討這件事情:RADseq缺值的tolerance差異,對於:1.用於分子親緣關係研究的基因組data set的性質; 2.對於單系群的推論; 3.對於分子親緣關係推論的精確度,有何影響。

最後在模擬的結果中可以推論(是的,我要跳掉中間有點難懂的模擬方法):若我們使用保守(conservatively)選擇的loci(=對於缺值低容忍度=矩陣較大=缺值較多)所得到的結果,會比randomly selected loci矩陣所建立的結果,在分子親緣關係的精確度來說會來得比較糟。在模擬結果中認為後者為優。在很多研究中也都能發現,像是慈鯛 [Wagner et al. (2013)] 和果蠅的例子 [Rubin et al. (2012)]。至於library construction和實驗材料本身的divergence history所導致的缺值影響,在這篇研究中並未有確切的結論,無法提供一個確定多少coverage的值或矩陣大小,依舊和不同物種有不同的情況而定。

這篇研究的價值在於,這是第一篇探討RADseq缺值的文章(尤其文中整理出那張可能導致缺值的圖很直觀,也很容易懂,我很喜歡用於解釋缺值成因),其中引用的研究其實不算多(可能也不夠多吧),用來模擬測試的現存研究也是,大略只能提出:the intuitive appeals about being conservative by removing loci may be misguided. 在我的研究亦有相似的推論。在我的研究中亦發現,在最後align好的SNP檔案裡,其實可以逐一打開調查不同矩陣的缺值比例,甚至一一去調查每個樣本的缺值比例,在我的結果中就可以發現,矩陣整體的缺值比例最多約為20幾%,最少有8點多%(STACKS分析),但是後者所建立的親緣關係結果就非常糟糕。以及一般而言,個體的缺值比例約為十幾%,但有少數外群的種類缺值比例會高達七十幾%(pyRAD分析),而在之後便會將它剔除,不納入矩陣中分析。然而我目前所得到的RAD的raw data亦只能透過後端參數品質篩選來挑整缺值比例而已。目前也尚未找到關於缺值要如何deal的確切標準相關研究,應該要繼續追蹤。

這篇研究中幾個名詞的概念整理:1.矩陣,在英文中常用data set, matrix。2. 對於納多/少缺值矩陣=big /small data set,較少缺值的矩陣=low tolerance for missing data (對於缺值低容忍度)=conservatively selected loci (保守),含多缺值的矩陣=randomly selected loci。

提一下supermatrix(常看到卻不懂之釐清):In the supertree method, different systematic data sets are analyzed separately, and then the trees derived from these independent analyses are used to produce a single, joint estimate of phylogeny.(先解釋supertree)The alternative supermatrix approach instead involves combining all systematic characters into a single, giant phylogenetic matrix and then analyzing all the characters simultaneously(supermatrix就是把所有的systematic characters一起考慮,在分子中就是同時計算所有串連序列特徵,在RAD的資料中勢必要使用這樣的運算策略)[詳情請看de Queiroz, Alan, and John Gatesy. “The supermatrix approach to systematics." Trends in ecology & evolution 22.1 (2007): 34-41.,有圖示易懂]

 

Huang, H., & Knowles, L. L. (2014). Unforeseen consequences of excluding missing data from next-generation sequences: simulation study of RAD sequences. Systematic Biology, syu046.
Rubin, Benjamin ER, Richard H. Ree, and Corrie S. Moreau. “Inferring phylogenies from RAD sequence data." PloS one 7.4 (2012): e33394.

[未讀]
McCormack, John E., et al. “Next-generation sequencing reveals phylogeographic structure and a species tree for recent bird divergences." Molecular Phylogenetics and Evolution 62.1 (2012): 397-406.
Zellmer, Amanda J., et al. “Deep phylogeographic structure and environmental differentiation in the carnivorous plant Sarracenia alata." Systematic Biology 61.5 (2012): 763-777.

[RAD seq classic] 輻射演化的親緣關係研究。

開始我的小檗研究之前,這篇2013年關於慈鯛的研究大概是我的精神指標,研究背景和我的材料極像。這篇研究的標題下得聳動“unprecedented resolution”,對於當時剛開始接觸演化領域的我來說,不太懂這有什麼好驚喜的。

經歷快速輻射演化的物種對於演化研究來說都是件不容易的事情,東非維多莉雅湖的慈鯛就是很經典的例子。雖然他們的外觀型態有所差異,但很常因之間的DNA變異累積不夠多,只用少數的分子標記(maker)就難以將其解開他們種間的界線或親緣關係,再加上我們的研究常面對沒有reference genome的非模式物種(Non-model organisms),這會增加這類型研究的困難。然而RAD seq補足了Sanger定序這方面的缺陷。用更有效率(包括節省金錢和實驗勞力)的方法,得到更多的分子標記,進而找到更多的SNP,便能得到更好的結果。當然之中還是有一些問題,之後有時間再來討論。

在這篇慈鯛的研究中,便使用了RAD seq 方法 (Etter et al., 2011) 搭配次世代定序 (NGS, Illumina HiSeq 2000),使用Sbf I限制脢切割,用了三個lanes定序了156個個體(使用 52和60 barcodes, 六個mer),之後再補兩個lanes來增加個體的coverage(由於此用RAD方法會受到許多因素的影響導致每個樣本最後得到的資料量會有差異,通常都會在補data來使整體研究完善,有什麼因素之後再詳細說明。)大體而言,我所使用的方法和這篇的實驗制備流程也差不多,看來輕鬆簡單,對於一個沒有實驗經驗的我來說,就卡關了好幾個月,畢竟從DNA品質開始就有所要求(當然是越新鮮越好,和過去Sanger所能使用的品質有所不同),接下來RAD library地製備更是折磨(想到就不寒而慄),當然克服之後又是海闊天空。

接下來得到raw data後,研究先使用了FastX toolkit先做了品質過濾(end-trimmed to a length of 90 bp and reads containing one or more bases with a Phred quality score below 10, or more than 5% of the positions below 30, were discarded.),接著再使用STACKS再做一次,救回barcode有一個mer定序錯誤的序列,接著做de novo assembly ,依造所設定的參數[minimum stack size of 125 reads (-m),maximum distance between stacks (-M) within a locus as 2。] 這裡可以稍微解釋一下m值和M值,這也是主要在做RAD資料品質篩選常用到的參數,無論是STACKS或pyRAD軟體的概念都一樣。當我們最後得到RAD raw data,就是大量短序列的讀長(稱為“read”,前面會包含barcode序列和平台序列,長度不一定,憑定序決定,約莫100 bp),這裡可以稍微想像一下:透過軟體參數將同一個個體中一樣的reads疊起來 (read depth) 形成一個Stack (一個棧),那要多少一樣的read疊起來呢?這個參數就是m值。那不同個體間所組成的Stacks,差異多少個base才能被當作是候選同源基因座(contig或putative loci)並且輸出以供後續使用?這個參數就是M值。接著軟體也會自動除去過分高coverage的putative loci(會被稱為lumberjack stacks),因為這可能會是來自highly repetitive regions,便有可能會是non-orthologous sequences,所以次步驟也會一併除去這些潛在影響因子。接著此研究使用Genome Analysis Tool kit得到Genotypes:包括SNP和nonvariable sites,再透過設定不同的 ‘min individuals’ 來調整缺值(missing data)的比例(此研究中設定125, 115, 110, 100, 75 and 15 ,共145個個體。)以得到不同的矩陣(這也是在RAD研究中常被拿來探討的環節)。最後這步,我在使用STACKS的時候,使用‘population’參數,把一個個體當作一個族群,亦可以達到‘min individuals’的設定需求。在以上段落的步驟,大概使用STACKS和pyRAD都能做到,尤其pyRAD就是專為這類研究的所設計,就會更加便利。(STACKS的使用研究範圍比pyRAD更廣,之後有機會再來做兩者比較。)

在最後的結果中,除了釐清了種間關係以外,另一個在分析的結果中有個比較有趣的結果:就是缺值比例較高(矩陣大)的矩陣,他最後所得到的親緣關係結果支持度和結構,都會比參數設定較嚴格(矩陣小,缺值比例較低)的矩陣結果更來得好。其中也對於RAD data提出了一些問題:這些loci來自於哪裡?“currently unclear” ,由於當我們將此方法用於沒有參考序列的非模式物種,故無法將我們擁有的序列做比對,這會影響到在建立親緣關係時,所需要選擇的modle等後續分析問題。以上兩點也是我目前感到興趣的問題,在我的研究中也有遇到,之後亦可再加以討論。

本文主要討論了RAD seq在親緣關係研究上的主要大綱,以及其資料的處理流程。還有許多細節(像是實驗操作, 分析軟體和設定)有空再一一補上。

Wagner, Catherine E., et al. “Genome‐wide RAD sequence data provide unprecedented resolution of species boundaries and relationships in the Lake Victoria cichlid adaptive radiation." Molecular ecology 22.3 (2013): 787-798.

Etter, Paul D., et al. “SNP discovery and genotyping for evolutionary genetics using RAD sequencing." Molecular methods for evolutionary genetics (2011): 157-178.

20170424 完成論文之後的雜談。

終於在1/24 完成論文口試,大概告一個段落。

到當兵之前的幾個月(3, 4, 5月)就被召回中研院整理研究,希望能轉化成paper。

但很快就知道還差了一大段,關於分析和結果的穩定性。能畢業只是代表一個階段的任務達到而已,我想還有一大段路要走才行。

逃了很久的試驗紀錄,這裡也就停滯了一段時間。想整理一些關於從零開始的研究歷程讓別人不小心google到然後參考,但遲遲沒有動筆。

想要打幾篇看完paper的紀錄,也當作我的筆記。

【關於一個階段性任務的結束。】

從18分鐘口試的突然結束到親手拿到一張薄薄的失業證書,以及把所有家當搬離台北,就這樣不小心沉澱了幾個禮拜。

遙想這個研究的開始,是2014年的那個春天,和游某在科博館的麥當勞開始,還記得我點了過於甜膩的冰炫風,那大概是接觸小檗的第一個記憶的味道。看著學長把各種小檗的葉和花一一畫出,那是我第一次開始認識這群高山的謎樣植物,當學長結束幾頁筆記紙的介紹後,拋下一句話:『你覺得你可以在這之中問什麼問題?』當時只覺得問這個問題的本身,好難。我甚至不知道DNA要怎麼從一棵植物"弄"出來,然後又要怎麼變成這個(我還沒想到的)問題的解答,我只覺得腦中一片混亂,然後眼前冰炫風就這樣融光了。

接下來不怎麼踏實的訓練就開始了,穿插好幾次高海拔稀薄空氣的沁涼、中海拔的泥濘和血漬、化學藥品的刺鼻、指令碼的乏味,甚至是不同時區的氣味,揉合了無奈和失落的巴掌,以及過於輕忽和得意忘形的嘴臉,難免銘心但享受。現在回頭看著『小孽的路,順順地走~』的字跡,絕對是個威力無比的祝福。相較這一路所相伴的幸運,我的努力就變得非常渺小,就像那晚南湖圈谷那夜的滿斗星空所包容著,無比奢侈。

向來對於階段性結束或開始的慶賀或制式化的典禮不感興趣,甚至會想刻意逃避。或許只是潛意識越來越討厭成為焦點或避免起伏太大的情緒轉變,因為日子依然繼續,那些重要的時刻和重要的人,我都會一直放在心底,並用有限的腦力將你們感激,直到阿茲海默症發起。

最後,感謝山,感謝台北。

20170217