【Data process】2016/6/13

分析方向:看了一些paper,看到Stacks所設的m值:一個stack所需之最小定序深度 (read depth),其range大概在10~15左右。然而我測試是5,偏低。所以,想要開始測試高一點的值,看是否會影響樹型、分群或支持度。

48 samples (-m 10; -M 1 ; -n 1; -p 35) ‘allsample’

未命名.png

在台灣和南台灣小檗已經混在一起。清水山小檗被拆掉了。但各支序的支持度算高。

48 samples (-m 10; -M 1 ; -n 1; -p 40) ‘v1’

未命名.png

調高p值 (輸出的基因作至少須出現在多少族群中):35 → 40。發現支持度更差了。

台灣和南臺灣小檗依然為grade狀。高山和南投依然混在一起。所以我想要挑選各種下不同族群的Sample。台灣就挑合歡山、小奇萊族群;南台灣就挑北大武山族群;高山挑倫太文山、帕托魯山;南投挑合歡山小奇萊。得下面結果:

32 samples (-m 10; -M 1 ; -n 1; -p 40) ‘v2’

未命名.png

挑掉一些看似有問題的sample (不同地點的pengii和kawakamii,以及m值調高後,清水山就分掉了),便可發現各種間的分群就更好了;支持度亦高。

Next:

我認為可能是m值一下子調太高,清水山小檗本身的DNA品質就沒有其他人好,而在m值門檻調高後,沒有這麼足夠深度能組成一stack,就導致其在後面分析時的missing data較多,在樹的支序上就會誤判。以上是我的推測,所以我想接下來就逐漸調低m值 (8, 7…等)。之後再調整不同Sample的挑選,分開討論。

*分開討論:

  1. 先挑穩定且確定的各種族群,或是模式產地族群,優先分析。先測試前人的種的假說。
  2. 放入所有samples後,在其他的討論 (如:其實台灣的小檗還有很多新種?之類的,以及演化問題。)

【Data process】2016/6/10

和游學長討論,目前我的分析方向:

  1. 是否為好種?
  2. 演化關係?(外群問題)

從"外群"下手,kawakamii和pengii是否成clade的三大狀況:

  1. 在RAxML中指定wuyiensis為外群。
  2. 在RAxML中指定baradana為外群。
  3. 不加入wuyiensis,無外群。

 

結果:

↓ 在RAxML中指定wuyiensis為外群。outgroup_wu (20160610v2)

未命名

kawakamii和pengii為grade狀 (16, 21)。

↓ 在RAxML中指定baradana為外群。outgroup_bara (20160610v3)

未命名

kawakamii和pengii為grade狀 (16, 21)。

↓ 不加入wuyiensis,無根樹,無外群。no_wu (20160610v4)

未命名

kawakamii和pengii為一clade (16, 21)。

 

小結:wuyiensis似乎會影響kawakamii和pengii是否為clade或grade。到底是wuyiensis距離kawakamii和pengii的關係是太近?或太遠?在學長使用ndhF+rbcL+accD+ITS+psbA-trnH的樹中:wuyiensis距離kawakamii和pengii的關係是遠的,只是直線距離近而已。所以有可能是"長枝效應"而導致在資料上的誤判?(因為太遠,ATGC的變換已過一輪,又變回原來,兩者間在鹼基上的類似,誤判其關係太近。又尤其RAD的資料相對細緻。) 要去找解決"長枝效應"此狀況的方法 (突變率?) 又或者要如何解釋無根樹?(討論限制?有什麼是不能討論的?參考AFLP。)

在討論後續的演化前,要先弄清楚:樹本身的建構技術上的問題。

 

20160610(kawapengii)

↓ 只放所有的kawakamii和pengii,還有wuyiensis為外群。

未命名

4 (來自向陽山) 以下 (1~6),為kawakamii;以上 (7~12) 為pengii。關於他們的分布,或許有有趣的故事。pengii和kawakamii之間似乎沒有明確的species delimitation 。

 

NEXT:

  1. 放入所有pengii和kawakamii,並且先不放wuyiensis,嘗試建樹。
  2. 試跑STRUCTRUE,怎麼色塊分不開啦!!!@@ 齁 問題在哪。

【Data process】2016/6/7

20160606 和中研院越南學長討論完回來後,調整了參數 (populations的 -t 的參數),晚上重跑。

哈哈, 學習有時需要試誤呀!
不過我也剛剛也想我說得太滿, 一是你大概還不會估算需要多少記憶體? 二是
Yandi其實不時還有其他人再跑東西, 所以要給你用也是有點擔心啦!
我看到你已經開始在AgNEC1跑, 目前已經用到近70%的記憶體, 可以用"top"這個指
令看到互動畫面, 按"q"可以跳出來
如果又沒辦法登入, 大概有出事了, 再email給我吧

明後天我應該都會在農藝館318寫論文, 如果你需要問些操作上的問題, 可以來問
問, 但不能太久~

(20160606 甫錦學長的mail通知)

20160607

成功收到data,果然是"populations的 -t 的參數"問題。

*注意:Stacks出來的phylip檔的最後面,會多一行:"# Stacks v1.34; Phylip sequential; June 06, 2016″這類東西!要記得先刪掉,才能給Cipres RAxML吃!不然會出現error:顯示"最後有垃圾"而無法跑。

再接著跑Cipres RAxML-HPC BlackBox

得到RAxML_bipartitions的result檔,給FigTree吃。(我人生中的第一次,從葉子到DNA到Library到RAD data到Stacks到RAxML分析到出現樹)

未命名

自製的種階層簡圖如下:

TEST

和游學長的討論:

  1. 和原本用葉綠體跑的樹的假說做比較。原來kawakamii和pengii自成一clade,是明確的。但在RAD的樹中,kawakamii和pengii同種的同個體卻沒有成衣clade,而是有明確的分化現象;兩種之間也沒有形成一個clade。估計他們的族群內結構已出現。先使用當時學長所取樣kawakamii和pengii的地點的個體重跑樹,看是否有成為一個clade。於是出現k2p5的測試,挑選kawakamii2和pengii5,分別來自大雪山和北大武山 (為游學長當年做葉綠體的取樣)。結果出現,確實形成一clade。故將其挑出,再重跑一次。
  2. 除了尖萼clade,圓萼clade內部有了結構了。但因為一開始這部分都沒有解析度,所以可以重新開始寫我的假說。
  3. 定年問題:整合到葉綠體的樹裡來定年。
  4. 目前發現:baradana在裡面而非外面!

你要再繼續explore taroko和bara的關係,要很solid。因為這攸關 菲律賓高山植物的起源。這個關係很重要是因為 照目前你的圖看起來,菲律賓小檗是台灣傳播過去的,這樣我們等於某種程度上解開了菲律賓高地植物相的起源的一部分!是不是很興奮 百年來沒人證實過的!!

(游學長)

 

Next:

  1. 是否能跑出kawakamii和pengii同一clade。
  2. kawakamii和pengii種內、不同地點個體的結構問題。
  3. 可以先不看wuyiensis,他似乎比想像中的和kawakamii的關係還要近。似乎可以用無根樹呈現。

((((以上希望能在端午連假做好。))))

Next’s Next:

  1. 要到美國Botany 2016的poster的進度界線。(自己先想好,要和學長老師討論。)
  2. 其他分析方法的學習。(ex: structure)

【Data process】2016/6/6

【Data process】2016/5/24

大致上知道一些出錯原因和幾位學長的討論,如下。

沒有給population的PopulationMap → 最後給的phylip裏頭沒東西。

所以要分不同population,20160531時就放了PopulationMap檔案,讓一個個體一個族群,因為不確定分群 (不知道種是否為好的預先分群方式),無法自訂分群,只能先這樣跑,大概看一下誰會組成一個clade。後續再來分群跑。

Cipres不知為何出來後的檔案並不完整,並未出現最佳樹。似乎是突然結束。

我想我給Cipres的檔為fasta檔。但正常來說應當為phylip檔。突然結束或什麼,可能和這有關。

“fasta檔裡面給的是所有你過濾出來之後的allele,就是每一個"stacks"。phylip裡面輸出的是SNP,就是你比對之後所有堪用的SNP。phylip檔中的1,2,3..那些東西應該是可以直接output出population map 裡面的名字 (就是aristatoserrulata,breviesepala….)" 嘉偉學長的指教。

phylip檔:出來的東西就是已族群為單位,每一族群編號後的一列,應當就是最後align完的一整條序列。然後才餵給Cipres吃。如果想要得到一個個體的SNP,那就是把一個個體當一個族群來跑。(這算小撇步嗎?)

“population" 的program導致農藝系server當機問題。

“郁嵐, 我在上周六下午已經重開, 但目前看來, 我又無法遠端進去了,
上周六下午經過我檢視, 我想主要的問題是因為你們跑的東西把實體記憶體用光,
然後開始使用虛擬記憶體, 才會造成這樣的結果,
想請問一下, 你們最後跑的"population"程式, 是獨立的程式嗎?
如果是的話, 我建議你還是另外到Yandi去跑, Yandi有512GB的記憶體,
AgNEC1有96GB, 對一般的工作是很足夠, 但顯然你們跑的是特殊案例。
如果還需要幫忙, 再跟我說。"20160606 甫錦學長的mail通知。

populations -P ../ -b 1 -M ./PopulationMap.txt -r 0.75 “← 一個population只有一個sample沒有意義設這個" -p 5 “← 太低了!!!!至少要半數,ex:48個至少要設24或35″ -k -t 8 -f p_value –structure –phylip –vcf –fasta   20160606 越南Hùng Nguyễn 學長 (洪志銘老師研究室) 的指教。

參數門檻值設太低 → 最後篩選的資料量過大 → 才讓server當機!(目前的預計應為如此)

越南Hùng Nguyễn 學長:『如果denovo_map.pl能跑,照理說population也要能跑。因為他都能在自己的電腦裡面跑了。』『denovo_map.pl跑完後出來的三個檔案:tags, snps, allele』我應該要試著去看裏頭有什麼。『denovo_map.pl設的-m 5 -M 1 -n 1 ,他認為或許-m可以只要2, 3、M和n就設2試試看。』『因為我有一些個體,如morri或alpicola,只有單一個體,若往後要分群,其他的如kawakamii族群就有8個個體,會導致最後可用的locus數量少很多。為降低這樣的bias,就希望能盡可能把其他種的個體挑少一些,以平衡大家都差不多。反正最後也只是以族群 (我是看成種) 為單位。』『Stacks已改版,可克服pyRAD所勝出的缺失 (能否偵測到indel)、可知道不同的種的Depth coverage、BAM檔的問題 (似乎是用在Reference map所需)』

Cipers上面有不同的RAxML版本 (RAxML-HPC BlackBox、RAxML-HPC2 on XSEDE、RAxML-HPC2 Workflow on XSEDE、RAxML-HPC v.8 on XSEDE),以及之中你有下那些參數?

“我们一般用RAxML-HPC BlackBox,只用设计时间,碱基突变速率GTR+G或者+I+G,其他的都不用设置。" 刘路贤 學長的指教。

 

現在就等農藝系server復活了!Q Q