想入門RAD-Seq卻不知從何開始讀的文獻推薦

這篇寫給剛要入門NGS和相關實驗、或是要做RAD-Seq (或Target-enrichment) ,卻不知道要從何找起資料的你:

先釐清幾個名詞概念:

次世代定序Next Generation Sequencing (NGS) :主要以massively parallel sequencing的概念建構的高通量技術,達到同時高速大量的核酸定序。簡單說,他就只是一種定序的概念,不同公司就發展出不同的定序策略,以達到大量、平行定序的目標,常見的公司像是illumina系列。(這部分的介紹、優缺點、和Sanger定序比較很多,網路隨便google都有,中文資料也很多。)

圖書庫 建庫 Library Construction:要餵給NGS定序機器之前的實驗步驟(可以自己做或給別人做)依造你的材料或目標,又可以分成許多不同策略:像是RAD-Seq或Target-enrichment,底下也還有細分成不同方法,如[RAD seq之外] HyRAD & HyRADX就屬於在這個階層的範疇。

好,再來就幾篇文章可以先看:

RAD-Seq的原始使用,並非用在解決Phylogeny的領域,所以這幾篇介紹文獻,如2008年Baird等人寫的第一篇、第一次發表的RAD-Seq的文章,並沒有很好理解,畢竟他主要不是在做Phylogeny,但就經典就看看吧,算是對於RAD-Seq透測了解的知識:

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

Davey, John W., and Mark L. Blaxter. “RADSeq: next-generation population genetics." Briefings in functional genomics 9.5-6 (2010): 416-423.

接下來推薦一系列的文章,或是他其實是一本書。這本裡裡面有8個章節,八個大主題,很淺顯易懂,加上主題就是環繞在我們領域相關(Plant Systematics)以及NGS,讀起來也會比較輕鬆。雖然有點太淺顯,僅適合當入門閱讀,若要深入瞭解一些細節,就需要再自己找資料。裡頭的第六章就是介紹RAD-Seq。但可能要找找看,似乎太不好拿到電子檔。

Regnum Vegetabile Volume 158 Next-Generation Sequencing in Plant Systematics

這篇Nature的Review寫得很不錯,主要介紹了幾種RAD-Seq的不同變形。(我都戲稱他們為RAD-Seq家族),可以在弄懂了RAD-Seq的方法之後,擴展橫向知識的文章:

Andrews, K. R., Good, J. M., Miller, M. R., Luikart, G., & Hohenlohe, P. A. (2016). Harnessing the power of RADseq for ecological and evolutionary genomics. Nature Reviews Genetics17(2), 81.

關於使用RAD-Seq解決 維多利亞湖 慈鯛輻射演化的例子,這就是關於Phylogeny了,加上這類群又是非模式物種,和我的例子比較接近:(現在這類的應用已經非常廣泛,隨便找都有,只是就我讀的第一篇,內容也不會太艱澀)

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.

再來就是在實驗操作RAD-Seq上(可以只看材料與方法部分即可),我是先看我的共同指導老師的文獻(中英文各一篇):

謝明修, 吳東鴻, and 陳凱儀. “使用限制酶位點標定之核酸定序法進行稉稻雜交組合之穗上發芽數量性狀基因座的遺傳定位." 作物, 環境與生物資訊, 11 (1): (2013): 11-25.

Chen, Ai-Lin, et al. “Reassessment of QTLs for late blight resistance in the tomato accession L3708 using a restriction site associated DNA (RAD) linkage map and highly aggressive isolates of Phytophthora infestans." PloS one 9.5 (2014): e96417.

 

以上幾篇就是先推薦給剛入門卻不知道要從何找起資料的你,可以先看看。這幾篇都看完之後,就會慢慢有方向,再來就可以去看看新一點的文章。RAD-Seq相關的文章,幾乎每個月都有新的文章,每次稍微找一下,都有最新的文章出爐(關於分析軟體、新的實驗流程、門檻值設定和缺值的探討、用在不同族群的應用等),最後就會讀出心得,尤其一開始對於一些專有名詞的混亂,這可能要稍微適應一下,在NGS這個領域到現在,無論是RAD-Seq或Target-enrichment,很多相同意思不同字,或是不同意思或同一字的情況很多,端看那篇文章採用哪一種,但也可以自己做一些筆記,幫助自己釐清。對了,上面有提到的Nature的Review(2016)中,旁邊的box裡有專有名詞解釋。

接下來我應該會再寫一篇關於RAD-Seq的Raw data整理和分析相關的文獻整理,先醬。

廣告

20180614 Linux 基本 筆記

[Linux] 新增修改刪除 管理 使用者/群組 指令

Linux 的 passwd 指令範例教學

主要會使用到的指令:

0.登入server

ssh barbara.pena@140.109.29.210

然後輸入密碼

 

1.指定 UID

假如我想創建使用者 test002 , UID 設為 1022

# useradd -u 1022 test002

 

2.變更使用者密碼

一般的使用者執行 passwd 即可變更自己的密碼:

passwd

而變更密碼之前,必須先輸入現行密碼:

正在變更 gtwang 的密碼。
(目前的)UNIX 密碼: 
輸入新的 UNIX 密碼: 
再次輸入新的 UNIX 密碼: 
passwd:密碼已成功地變更

如果是 root 管理者的話,可以變更任何使用者的密碼:

sudo passwd gtwang

 

3.到/home下面開啟個人資料夾

bochung@bochung:/home$ sudo mkdir /home/barbara.pena

並開啟他的權限

bochung@bochung:/home$ sudo chmod 777 home/barbara.pena

 

 

 

目前用戶:

id bochung
uid=1000(bochung) gid=1000(bochung) groups=1000(bochung),4(adm),24(cdrom),27(sudo),30(dip),46(plugdev),110(lxd),115(lpadmin),116(sambashare)

id yuhsin.tseng
uid=1001(yuhsin.tseng) gid=105(systemd-bus-proxy) groups=105(systemd-bus-proxy),27(sudo)

id LiuSH
uid=1002(LiuSH) gid=105(systemd-bus-proxy) groups=105(systemd-bus-proxy)

id champ
uid=1007(champ) gid=1007(champ) groups=1007(champ),27(sudo)

[RAD seq之外] HyRAD & HyRADX

之前妤馨學姊有分享capture,target enchantment的研究,最後他有提到關於HyRAD和HyRADX這兩種方法。我整理目前一些常見的方法,簡單和大家分享。

無論是要聊capture或是RADseq,一開始會先介紹Reduced Representation Library (RRL)的概念,以獲得Loci的方法來看,這樣的Libaray製備的概念,大概就介於sanger定序和全基因體定序之間,可能會覺得不需使用到全基因體的資料,但又覺得sanger定序的資料太少或是太沒有效率,就會偏向使用這樣的方法。RRL這樣的概念,可以用於random或target的基因體上。我們比較常見的,也是之前我和宜軒用的方法,像是RADseq就是,或是之前陳老師有做的ddRADseq,或是像林老師他們做的GBS。

這裡簡單介紹一下RADseq的概念。簡單說就是靠限制酶的切位,定序這個site的前後一兩百bp。所以理論上來說,這樣就可以獲得每一個樣本,一樣的loci進行後續的分析。但事實上,並不是所有樣本都可以這麼順利。有可能是coverage的不足,或是本身的DNA品質不好,或是根本切位已經突變,就無法被取樣到。最後還會因為資料的處理,導致有些明明有定序出來的無法使用。一來這些都是變成missing data的原因。二來這些其實這些沒有用到的序列也都會被定序。因為一個lane的資料總量都是固定的,如果我們訂了一大堆不需要也不會用到的,相對也壓縮到需要的loci定序深度。最後拿到可以用的matrix時,缺值的比例就會很高。大概都高於50%以上。這也是為什麼RADseq會讓人詬病的原因之一。

另外再提另一個也屬RRL的方法,就是妤馨學姊做的capture,target enchantment的方法。上次他已經詳細介紹過了。簡單說他就是透過設計好的probe label目標DNA,最後用beads抓取這些被label的DNA,把不要的洗掉,最後就可以單純拿到target的序列。

這裡我整理一下,RADseq有幾個缺點:DNA品質限制,最後訂到的序列來路不明,大量缺值的問題。target enchantment的方法能彌補這幾個缺點,但是target enchantment的方法本身也有一些缺點,像是要研究從零開始的非模式物種,就需要花多一些時間和金錢來設計這些probe。所以就出現了介於這兩種方法的折衷。

像是Rapture,但這個方法其實沒有什麼特別,就只是在普通的RADseq library的製備流程中,在adaptor上加上biotin,可以更精準地抓到目標,減少雜訊,最後也可以提高locus的coverage。像是REAL baits,他其實和等一下要介紹的HyRAD很像(RADseq-based capture),但這是GBS-based capture。或是像這是上禮拜茂綸學長有在我們研究室社團分享的EecSeq,這就有點像等等要介紹的HyRAD-X,但是在probe的library產生的部分,他不是使用限制酶去標定的,不會因為這樣限制probe的數量。反正就是有新方法就各取名字的時代,差別可能只有幾個步驟的方法組合差異,像是RADseq的變形方法,或是用什麼來設計probe。但其實概念都差不了多少。

先來看hyRAD,首先會先備製兩個library,簡單說就是一個是要被抓的,一個是抓別人的。要被抓的就隨便做就好,這裡就是用shotgun library,一些品質比較差的樣本就會放在這邊。另一邊是要設計用來抓別人,所以本身的DNA品質就要好一些,這裡使用ddRADseq的方法,最後得到的reads在自行biotin label,自製probe就完成了。最後就兩邊混合capture,整體步驟就和target enchantment概念一樣。

在這篇研究中,分析的時候測試了三種reference,RAD-ref是用要抓別人的那個probe的library,那去定序,理論上這個作為reference應該是最合情合理的。另一個把最後抓到的序列定序後denovo assamble拿到的contig當作reference(像在RADseq後續分析一樣),還有一個是同時參考兩個定序後的結果,將RAD-ref做延伸,使用PriceTI設定比較嚴格組出的contig作為reference。最後結果可以發現單單用denovo assamble的reference所得到的contig是最多的。由於在抓的時候,專一性若不夠高,就很有可能會抓到其他的,尤其又是自製的probe,比起大公司合成(像是MYbaits和IDT),效果會更不好。接下來以single map的結果來看,用較嚴格的RAD-ref-ext作為reference的表現最好。再從最後得到的SNP總數來看,以map RAD-ref最好。

另外這研究在製備shotgun library,也就是被抓的的那個,他也有測試了要不要sonicate,結果發現舊一點的材料不需要會比較好,但過於老舊的(58歲那種,好像有沒有sonicate影響不大了),至於新鮮材料好像結果也沒影響太大。這裡看到fresh的材料獲得的SNP數量最少,感覺很怪,但裡頭有解釋說因為這裡設計用來抓人的reference sample和被抓的新鮮sample都是新鮮材料,兩種的genetic distance本來就比較小,找到的SNP就比較少。接著看到的是Matrix fullness(其實就是missing data比例反過來說而已)可以看到其實結果並沒有差很多,都有超過50%,也沒有哪一個condition的設定比較好,在這個研究中他用最後得到的結果反推,認為RAD-ref, non-sonicated的分析結果比較符合真實情況。原則上我覺得這個方法並沒有大破大立的進展,但就:缺值比例下降和平均的coverage提高而已。

接下來看到hyRAD-X,原則和概念和流程都和剛剛的hyRAD一樣,只是在自行合成probe的地方,他是用RNA做biotin labeling而已。但這個研究比較有趣的是用subfossil作為材料,就是沉積在湖底的銀冷杉松針。時間就拉到距今五千多到七千多年前。再抽DNA時也很酷,還要先照UV光。這裡主要測試了三種方法,和hyRAD做比較。這兩個長條圖為有沒有移掉PCR duplicates,map到reference的read數量所佔比例。在去除PCR duplicates之前和之後,不同方法所map到的結果也不太一樣。但仔細看,這些read所佔的百分比其實超少啊。hyRAD-X主要表現比較好的部分大概也是coverage提高。讀到最後好像有點失望,但回頭想想他可是從七千多年前的沈積物中拿到資訊來分析,並且拿到500 SNPs,缺值比例少於1/3,來其實也就很厲害了。

 

Suchan, T., Pitteloud, C., Gerasimova, N. S., Kostikova, A., Schmid, S., Arrigo, N., … & Alvarez, N. (2016). Hybridization capture using RAD probes (hyRAD), a new tool for performing genomic analyses on collection specimens. PloS one11(3), e0151651.

Schmid, S., Genevest, R., Gobet, E., Suchan, T., Sperisen, C., Tinner, W., & Alvarez, N. (2017). HyRAD‐X, a versatile method combining exome capture and RAD sequencing to extract genomic information from ancient DNA. Methods in Ecology and Evolution.

【Experimental Log】2017/09/19-28

29個紅豆樹數的Samples:

Number

Species

Champ

peng24773

Ormosia emarginata (Hook. et Arn.) Benth. 凹葉紅豆

1

peng24782

Ormosia fordiana Oliver 肥莢紅豆

2

61B

Ormosia pachycarpa Champ. ex Benth. 茸莢紅豆

3

65B

Ormosia indurata H.Y. Chen 韌莢紅豆

4

66B

Ormosia semicastrata Hance 軟莢紅豆

5

69B

Ormosia pachycarpa Champ. ex Benth. 茸莢紅豆

6

70B

Ormosia indurata H.Y. Chen 韌莢紅豆

7

71B

Ormosia emarginata (Hook. & Arn.) Benth. 凹葉紅豆

8

KFC3472

Ormosia pachycarpa Champ. ex Benth. 茸莢紅豆

9

KFC3550

Ormosa pinnata (Loureiro) Merrill 海南紅豆

10

KFC3551

Ormosia henryi Prain 花櫚木

11

KFC3552

Ormosia glaberrima Y.C. Wu 光葉紅豆

12

KFC3554

Ormosia pachycarpa Champion ex Bentham 茸莢紅豆

13

KFC3555

Ormosia merrilliana L.Chen 雲開紅豆

14

KFC3556

Ormosia balansae Drake 長臍紅豆

15

KFC3564

Ormosia

16

KFC3567

Ormosia glaberrima Y.C. Wu 光葉紅豆

17

KFC3568

Ormosia semicastrata Hance 軟莢紅豆

18

KFC3569

Ormosia fordiana Oliver 肥莢紅豆

19

KFC3571

Ormosia fordiana Oliver 肥莢紅豆

20

KFC2991

Ormosia hosiei Hemsl. & E.H. Wilson 紅豆樹

21

KFC2992

Ormosia hosiei Hemsl. & E.H. Wilson 紅豆樹

22

XWB13254

厚莢紅豆

23

y2596

小葉紅豆

24

y2695

小葉紅豆

25

y3101

肥莢紅豆

26

KFC3089

Ormosia hengchuniana 恆春紅豆樹

27

Torke1573

Ormosia holerythra

28

Michelangeli2047

Ormosia costulata

29

DNA extraction :皆使用CTab (未秤重,目視約一片正常大小葉子,加入液態氮磨碎。) 最後總體積用50 uL去溶。(使用過DNeasy Plant Mini Kit效果不佳。)

跑膠結果:

 

DNA定量:nanodrop (除29以外,濃度和品質皆佳。)

 

PCR(MasterMix 12.5uL, Primer L端 1uL,Primer F端 1uL, DNA 1uL, NgCL2 1.5uL, H2O 8uL,總體積為25uL):

ITS4 – ITS5

(A)

trnL – trnF

(B)

trnK685F – trnK2R

(C)

matK4La – matK1932Ra

(D)

1

9/15/17

1

9/15/17

29

9/15/17

15

2

9/15/17

2

Double Bind  9/26 切膠

9/15/17

30

9/15/17

16

3

9/15/17

3

D 9/26 切膠

9/15/17

31

9/15/17

17

4

9/15/17

4

D 9/26 切膠

9/15/17

32

9/15/17

18

5

9/15/17

5

D 9/26 切膠

9/15/17

33

9/15/17

19

6

9/14/17

1

9/15/17

34

9/14/17

3

7

9/15/17

6

D 9/26 切膠

9/15/17

35

9/15/17

20

8

9/15/17

7

9/18/17

1

9/15/17

21

9

9/15/17

8

D 9/26 切膠

9/18/17

2

9/15/17

22

10

9/15/17

9

D 9/26 切膠

9/18/17

3

9/15/17

23

11

9/14/17

4

9/18/17

4

9/14/17

6

12

9/15/17

10

D 9/26 切膠

9/18/17

5

9/15/17

24

13

9/15/17

11

D 9/26 切膠

9/18/17

6

9/15/17

25

14

9/15/17

12

D 9/26 切膠

9/18/17

7

9/15/17

26

15

9/15/17

13

9/18/17

8

9/15/17

27

16

9/15/17

14

9/18/17

9

9/15/17

28

17

9/17/17

1

D 9/26 切膠

9/18/17

10

9/17/17

14

18

9/17/17

2

D 9/26 切膠

9/18/17

11

9/17/17

15

19

9/17/17

3

D 9/26 切膠

9/18/17

12

9/17/17

16

20

9/17/17

4

D 9/26 切膠

9/18/17

13

9/17/17

17

21

9/17/17

5

9/18/17

14

9/17/17

18

22

9/17/17

6

D 9/22 切膠

9/18/17

15

9/17/17

19

23

9/17/17

7

D 9/22 切膠

9/18/17

16

9/17/17

20

24

9/17/17

8

9/18/17

17

9/17/17

21

9/28/17

25

9/17/17

9

9/18/17

18

9/17/17

22

 9/28/17

26

9/17/17

10

D 9/26 切膠

9/18/17

19

9/17/17

23

27

9/17/17

11

9/18/17

20

9/17/17

24

28

9/17/17

12

 

9/18/17

21

9/17/17

25

 

29

9/17/17

13

9/18/17

22

9/17/17

26

純化:主要使用ExoSAP。若有切膠或單純純化,使用Viogen PCR Advanced PCR Clean Up System。

送定序 (PCR原液 1uL, 單邊primer 0.5, 水 8.5):四組primer(依順序為A, B, C, D)搭配不同樣品材料編號,為送定序代號。

目前:等定序結果。

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

上一篇用比較概念性地介紹了RAD seq,這一篇就說說RAD seq的壞話(error和bias)

這些偏誤的產生大致上可分成兩個原因,一類來自於本身次世代定序的問題:另一類就是本質上的問題(library的備製和後續的資料分析。)

1.  Allele dropout and null alleles. 等位基因遺漏和無效等位基因。導致這樣的原因,有幾個像是在樣本本身基因組的限制脢切位上發生多態性(polymorphism),這樣就會導致此片段無法被切到,而無法順利取樣到,就會此這段基因遺漏。或是SNP來自於無效對偶基因(null alleles)也會此後續分析發生問題。若是使用雙限制酶切割,會導致片段過長而在後續的片段篩選時被篩掉。所以也有一些模擬的研究認為,由於在切位上的突變所造成的Allele dropout影響,對於使用雙限制酶的RAD seq會比用原始RAD seq(單一限制酶)的影響來得嚴重。(另外亦會影響族群遺傳FST )

2. PCR duplicates and genotyping errors. 大部分的NGS Library的建置,original template DNA都會歷經一段PCR放大階段。由於隨機的PCR可能會導致原本是heterozygotes對偶基因在偏誤地放大下,在後續分析時誤解成homozygotes。在過去研究中發現PCR duplicates在RAD seq的data中發生機率高(20–60% of reads),雖然理論上PCR不會特定放大某些alleles,若真正發生在大量資料中的話,在分析上就會有問題。總之,RADseq建庫的過程中,一旦發生PCR duplicates,就會導致偏誤。然而這樣的問題可以在後續的資訊處理時解決,以增進genotyping的精確度。[看Harnessing the power of RADseq for ecological and evolutionary genomics 圖2,易懂。]

3. Variance in depth of coverage among loci. 由於PCR duplicates和allele dropout皆會導致genotyping errors,兩者都是在定序時,僅將特定(也就是最後留下來,preferential amplification,優先放大)序列定出,導致數量的偏誤。但這裡和前者不同的是,有時雖然都會被定序出來,但不同的reads在最後定序的數量(sufficient depth?)結果不一定會一致,在最後分析深度篩選時,便容易會漏掉某些深度較低的基因座。這樣的問題通常會發生在RADseq建庫時的PCR步驟,過去通常會說preferential amplification的發生和GC content有關;以及發生在較短的reads上(和長reads相比)。所以若以長度影響為影響因素的話,這樣的情形比較容易發生在由兩個限制酶切割的RADseq方法(reads長度不一致),至於2bRAD或原始RAD就不會有影響(因為他們的reads長度是一致的)。除了PCR步驟會影響以外,另一個影響會是在機械震碎步驟(mechanical shearing,在原始的RAD中的步驟),但是目前使用原始RADseq的研究中並未發現這有特別的影響,因為大部分使用rare cutters所裁減的片段,多>10 kb。但是但是,其實這在真正研究操作時,還是會遇到這樣的情況。在陳凱儀老師那就發生過了,主要還咬考慮到一個因素:DNA本身的quality,若品質很糟的話(像是片段太破碎),會非常在震碎那一步,會比新鮮材料的DNA來得容易震短。[分析深度參差原因→PCR步驟(GC content)&短的reads(mechanical shearing or low DNA quality)]

 

Andrews, K. R., Good, J. M., Miller, M. R., Luikart, G., & Hohenlohe, P. A. (2016). Harnessing the power of RADseq for ecological and evolutionary genomics. Nature Reviews. Genetics17(2), 81.

 

之後來聊聊設計RADseq研究的前提和準備吧。

[RAD seq] 親緣關係研究分析(phylogeny & species tree)

關於RAD seq的data有幾個特性:1. 相較短的讀長(Relatively short sequence reads):2. 全基因組分佈 (Wide genomic distribution)3. 難以分辨旁系同源基因或直系同源基因(Lack of distinct, at the outset of project, between paralogous and orthologous sequence)4. 缺值(Dropout of loci and alleles)。若要我描述RAD seq data,白話來說:資料量大,來路不明,缺值。之前也提過了缺值,也是RAD seq最惱人(Vexing,學一個新單字)的。

由於以上特性(或問題),就導致在之後的分析會有某些限制,例如因為短讀長和缺值的特性,便很難使用過去用於多基因(multilocus)來建構完整gene tree的方法(包含建構species tree和concordance analysis。)所以這也是為什麼很多RAD seq 研究是使用串聯矩陣(concatenated matrix)來推論出親緣關係樹。截至目前為止,有很多的研究使用這樣的方法得到很好解析度的結果(和Sanger定序相比)但受限於串連後的資料量非常龐大,故僅有少部分分子親緣關係的分析軟體能夠運算。現在最常使用的,以及適用於這樣尺度的運算,便是RAxML,但目前亦沒有研究對此有所討論(此外大家使用此軟體時,最常用的substitution model也只是預設的GTR+gamma model。)也有一些研究將所有無變異的位點排除,這樣有可能會導致一些bias(畢竟那些程式本來設計就非屬:所有皆為有變異位點),若可以就使用所有的基因座,便可消除這類的疑慮。(以上是官方說法啦)但是在實際操作上,完整的data set(包含變異及非變異點位)的資料量會更加龐大,在運算時也會造成問題。(真是麻煩呢)所以說呢,使用串聯方法(concatenation approach)既簡單又方便,但很明顯就是將這些問題掃到地毯下而已(sweep the error under the rug),那些有可能會發生在loci間的變異:譜系學(如隨機合併 stochastic coalescence 或基因滲入 introgression);substitution
parameters 和 evolutionary rate等。然而這導致出(大家也最關注的)statistical inconsistency的主因(gene tree不等於species tree。)在更多實際案例和實驗案例做比較之前,我們對於串聯法所建構的RADseq親緣關係樹,應更謹慎並保持懷疑。至於關於ancestral demographic的評估推論,就更不適合用此方法。在 Next-Generation Sequencing in Plant Systematics. “Chapter 6:Inferring phylogeny with RADseq data”文中,並不建議使用串聯方法(concatenation approach)來推論最終的親緣關係樹,最多他僅適合用來作為啟發工具(heuristic tool),僅提供一個基礎或框架,還需要更進一步分析(像是introgression)。目前很少有更進一步對於RADseq和串聯方法的討論。

所以所以所以,到底RADseq的方法能不能推論species tree呢?其實很少有研究在討論這件事,目前也並沒有特定的軟體專門用來處理RADseq的資料,但能預估有越來越多的研究需要這樣的方法。這裡先稍微整理一下以往建構species tree的方法,主要可分成兩大類:Supermatrix (concatenation) 和 Supertree,前者為先將所有基因座串聯後,再用此串聯後的序列建構出一棵樹,方法通常就是前者所提及的maximum-parsimony, maximum-likelihood, Bayesian inference and distance methods;後者supertree的方法,也是我比較不熟悉的,是將各個基因座分別建樹之後,再使用summary methods最後合併成一棵樹。據學長姐說,後者的方法近年就越來越少使用。

目前我讀過幾個研究,所使用的軟體有:BUCKy (他最後結果的呈現,可以在樹形上,不同枝上有百分比數值,顯示此枝得到多少genome比例的支持,但他有個壞處在於:不容許有缺值,故近年來RADseq研究使用他的頻率較少。), SNAPP(此為BEASTS的外掛軟體,主要需要的資料為unlinked bialletic markers類型,像是SNP或是AFLP。最後呈現方式會像是取不同locus來跑出gene tree,最後再疊在一起。並可以再加另一個外掛程式,可計算BFD*來計算不同的species model的BF和ML值,並決定哪一個model比較好,並排出rank優先順序,來決定哪一個分種model比較合適。), SVDquartets(一次取四個物種來建樹,並且最後將所有結果其合併。算是簡化計算,可以增快計算效率。並且可以考慮到RADseq中缺值多的問題,若一次僅討論四個,比起一次討論二十個,所遇到的缺值比例會下降很多。), TreeMix(此軟體主要能呈現introgression events,也就是說一開始就可以先設定預計有幾次introgression事件,並且將其呈現在樹上。), SplitsTree(做出network。

在這裡可以很清楚了解到串聯方法所呈現的親緣關係樹並不適合用在RADseq的最終討論中,但現在所看到的RADseq親緣關係研究中,大多皆有呈現這個結果,並且再搭配其他使用supertree的方法,在一起討論。真希望可以快點看到RADseq的專用分析軟體啊。

 

Ree, R. H., & Hipp, A. L. (2015). Inferring phylogenetic history from restriction site associated DNA (RADseq). Next-Generation Sequencing in Plant Systematics’.(Eds E Hörandl, MS Appelhans), 181-204.

[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.