[RAD seq classic] 關於RAD seq (1)。

一開始從零開始,接觸到分生研究就直接碰到了RAD seq。很多時候少了基礎累積,就很難在上面堆疊,半路出家就容易遭遇到半生不熟的尷尬。又加上時間的限制,總需要在極短時間中找到一個比較有效率地權衡,先把需要的部份盡可能補齊。我想所有事情從零開始,最困難的部分就是概念上的建立。

RAD seq這個技術第一次發表是在2008年Baird等人的研究之中,在個研究裡所使用的材料是模式物種的三刺棘背魚(threespine stickleback)和粉色麵包黴菌(這篇值得看部分:方法部分和RAD的概念)。其實RAD seq可用的研究範圍很廣:Genomics of adaptation, Inbreeding and genomic diversity, Effective population size (Ne ) ,Population structure, phylogeography and conservation units, Introgression, Phylogenomics (Andrews et al., 2016),我們所使用的部分就只是其中的一部份而已。當然也必須分辨這個方法受到質疑時,問題是出在哪一個環節,而非全盤推翻。舉個例來說,Lowry在2017年的Breaking RAD和Catchen的Unbroken的爭執和討論,主要都在Genomics of adaptation和這方面的研究上,對於我的研究來說就沒有這麼直接相關。

想當然而,從2008年RAD seq這個方法被發表到現在,此方法也針對許多需求而有不一樣得調整(RADseq family像是2bRAD, ddRAD, GBS, CRoPS等,爾後有機會再來好好介紹)。但其中有幾個概念是固定的:1. 所有的方法都從genomic DNA開始→ 2. 並且會受到一至多個限制酶所切割→ 3. 接上Adaptor和條碼(Barcoding)→ 4. 片段大小的篩選→ 5. 最後定序得到RAD資料,以及資料處理。

1. genomic DNA開始:好的DNA品質、好的DNA品質、好的DNA品質,很重要說三次。如果DNA破碎會導致後續不好的影響,這影響層面很廣,像是在後續震碎步驟到挑選片段,或是在擴增PCR時就容易會有bias。

2. 限制酶:大致來說,限制酶可以分成兩大類:common cutters和rare cutters,透過限制酶切位的片段大小,達到切割頻率的差異。在很粗淺的估計可以得知:8-cutter會每65,536 bp切一次,而6-cutter為每4,096 bp切一次。當然這還會關係到本身材料的DNA狀況,像是GC content。在研究開始之前,我好像沒有認真調查自己物種本身的DNA性質,這應該要先做才行。

3. 接上Adaptor和條碼(Barcoding):以原版的RAD seq步驟中,會有兩次接上Adaptor(分別為P1和P2)。第一次的Adaptor(P1)包含擴增primer、上機所需的定序primer和條碼,並接在限制酶切位的地方。接下來會先透過隨機的震碎(mechanical shearing),讓大部分的片段達到所需的大小後,會再接上第二種Adaptor(P2),此功能為ensure that only the target genomic DNA fragments are sequenced. 簡言之,就是同時要擁有P1和P2的Adaptor才會順利在定序機器上被定序出來。像我所使用的算是Illumina sequencing-based的RAD seq,P2 Adaptor就會被設計成Y型,以在Illumina平台專用的定序技術(PCR bridge)所使用。現在也有一些專門用於RADseq family的kit,他們的Adaptor就又會有所不同(像是NEBNext® Ultra™ DNA Library Prep Kit for Illumina®)。在原版的RAD seq步驟中,一開始把有不同條碼的P1 Adaptor分別接上不同個體的片段之後(Barcoding),便可以將所有的樣本都Pool在一起,在接下來步驟中同時操作(這裡常會用一個動詞叫作multiplex。)如此以來,便可以節省下許多時間和勞力。

4. 片段大小的篩選:在RADseq family在這一步有不同的策略,大多都是想辦法獨立出所需要的片段大小。會透過間接或直接的方法,前者像是PCR(GBS和CRoPS,尚未了解),後者像是直接使用切膠的方法(RRLs, multiplexed shotgun genotyping (MSG), ezRAD and double digest RAD (ddRAD))。但在原版的RAD seq中,就直接使用機械性震碎的方法取得一定範圍的片段大小,所以每個片段會有一端是限制酶cut site,而另一端是隨機震斷。之後便透過磁珠篩選出我們所要的目標片段大小,以符合上機標準。在此片段大小篩選步驟中有個關鍵可以討論,就是希望最後所取得的片段大小要越一致越好,若不一致會在最後定序時造成浪費,也更容易產生缺值。然而像是2bRAD此方法,便透過IIB restriction enzymes得到一模一樣大小的小片段(33–36 bp)[尚未暸解,感覺很有趣]。

5. 次世代定序得到RAD資料:依不同的定序機型會有不同的片段大小(有機會再來聊聊),現在大部分都使用Illumina,目前可得到的片段大小約為50–300bp,有人估計未來會越來越長。而大致上可以提供兩種定序服務:single-end sequencing和paired-end sequencing,前者僅定序forward端,而後者為forward read和reverse端。幾乎所有的RADseq family皆可以選擇這兩者的定序方法,但對於2bRAD來說似乎不太需要paired-end(畢竟片段太小,也無需兩端確認)。接下來得到資料後的篩選和處理,便會透過幾個基本的概念來達成:de-multiplexing and trimming of barcodes,接著大致上能分成兩個處理方向:有reference genome或de novo。已有套裝軟體,像是近年較常被使用的Stacks, pyRAD和UNEAK,Stacks所能用的分析範圍較廣(像是genotyping and calculating population genetic statistics),至於pyRAD主要設計為phylogenetic所用,主打能克服insertion–deletion的問題,以得到更精準的資料,UNEAK主要用在GBS的資料處理(目前我使用過前兩個。)由於大多我所使用的材料物種多為非模式物種,de novo組裝的策略在於藉相似性將類似或相同的read cluster起來,接著便可以找到locus,以及若藉由paired-end data可以得到更長一些的contig。這些較為細節的分析策略在之前的文章中有提到。接著就會遭遇到缺值的問題,有一些研究就專注於這部分討論,但目前尚無定論。

這篇筆記就先寫到這裡,主要參考Andrews等人2016年的文章的前面部分,這篇很值得一看,很適合用於大方向觀念的建立,亦有一些Box等資訊補充一些詞彙的定義。但還是有大部分的細節需要去深入瞭解後才容易看懂或有感。另外這篇review也整理了很多很新的RADseq family策略比較,可以瞭解自己所需的目的,挑選適合自己的方法。另外這篇文章接下來後半段的部分,會集中在使用RAD seq會遭遇到的問題或bias,預告給下一篇筆記吧。

 

Baird, Nathan A., et al. “Rapid SNP discovery and genetic mapping using sequenced RAD markers." PloS one 3.10 (2008): e3376.

Andrews, Kimberly R., et al. “Harnessing the power of RADseq for ecological and evolutionary genomics." Nature Reviews Genetics 17.2 (2016): 81-92.

Lowry, David B., et al. “Breaking RAD: an evaluation of the utility of restriction site‐associated DNA sequencing for genome scans of adaptation." Molecular ecology resources 17.2 (2017): 142-152.

[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 

【Experimental Log】2016/09/19

3 Samples:CP144 (Berberis chingshuiensis), CPL19 (Berberis alpicola), CPL20 (Berberis alpicola)

DNA extraction [CPL19, CPL20]:CTab (未秤重,目視約一片正常大小葉子,加入液態氮磨碎。) 總體積用50 ug去溶。

純化kit:Monarch™ PCR & DNA Cleanup Kit (5 μg) (NEB)

DNA定量:Qbit

Qbit結果:CP144 (Berberis chingshuiensis):7.48 ng/uL, CPL19 (Berberis alpicola):34.6 ng/uL, CPL20 (Berberis alpicola):38.8 ng/uL

跑膠結果

20160919跑膠.jpg

結論:壓力來了就煩躁。

【Experimental Log】2016/09/08-10

14 Samples:CP247 (Berberis kawakamii), CP248 (Berberis kawakamii), CP249 (Berberis nantoensis ?), CP250 (Berberis nantoensis ?), CP251 (Berberis kawakamii), CP252 (Berberis kawakamii), CP253 (Berberis kawakamii), CP254 (Berberis kawakamii), CP255 (Berberis aristatoserrulata), CP256 (Berberis kawakamii), CP257 (Berberis kawakamii), CP258 (Berberis kawakamii), CP259 (Berberis nantoensis ?), CP260 (Berberis nantoensis ?)

DNA extraction:CTab (未秤重,目視約一片正常大小葉子,加入液態氮磨碎。) 總體積用50 ug去溶。

純化kit:Monarch™ PCR & DNA Cleanup Kit (5 μg) (NEB)

DNA定量:Qbit

Qbit結果:CP247 (Berberis kawakamii):95.4 ng/uL, CP248 (Berberis kawakamii):86.8 ng/uL, CP249 (Berberis nantoensis ?):86.2 ng/uL, CP250 (Berberis nantoensis ?):too high ng/uL, CP251 (Berberis kawakamii):82.6 ng/uL, CP252 (Berberis kawakamii):52.2 ng/uL, CP253 (Berberis kawakamii):too high ng/uL, CP254 (Berberis kawakamii):59.6 ng/uL, CP255 (Berberis aristatoserrulata):49.4 ng/uL, CP256 (Berberis kawakamii):108 ng/uL, CP257 (Berberis kawakamii):112 ng/uL, CP258 (Berberis kawakamii):too high ng/uL, CP259 (Berberis nantoensis ?):too high ng/uL, CP260 (Berberis nantoensis ?):too high ng/uL

跑膠結果

20160910跑膠.jpg

結論:最近很疲倦。

【Experimental Log】2016/08/19-20

8 Samples:CP239 (Berberis deinacantha), CP240 (Berberis deinacantha), CP241 (Berberis deinacantha), CP242 (Berberis grodtmanniae), CP243 (Berberis dumicola), CP244 (Berberis sanguinea), CP245 (Berberis kawakamii), CP246 (Berberis breviesepala)

DNA extraction:CTab (為乾燥材料,均質機打碎。一片葉子一管。) 總體積用50 ug去溶。[CP243 , CP244的狀況不好,抽了兩次 (8/19、20)。]

純化kit:Monarch™ PCR & DNA Cleanup Kit (5 μg) (NEB)

DNA定量:Qbit

Qbit結果

CP239 (Berberis deinacantha):too high ng/uL, CP240 (Berberis deinacantha):110 ng/uL, CP241 (Berberis deinacantha):too high ng/uL, CP242 (Berberis grodtmanniae):110 ng/uL, CP243 (Berberis dumicola):46.4 ng/uL, CP244 (Berberis sanguinea):15.8 ng/uL *(CP224的顏色很髒,在原液時,就有很多不明析出漂浮物,純化時無法清除顏色。)

CP243 (Berberis dumicola):82.6 ng/uL, CP244 (Berberis sanguinea):36.8 ng/uL, CP245 (Berberis kawakamii):83.4 ng/uL, CP246 (Berberis breviesepala):too high ng/uL *(這次沒有放overnight,只有三個小時,葉子放的量也減少,純化後顏色皆為透明。)

跑膠結果:

P1010428.JPG

結論:感覺CP224可以純縮成一管。其他的狀況都還算OK。已經可以正確判別葉子狀況適不適合抽出高純度DNA了。不適合,就也只能以量取勝法,來得到所需濃度和品質。