Analyzing Diversity Data : วิเคราะห์ความหลากหลาย
Ekaphan Kraichak
August 15, 2023
วิเคราะห์ความหลายไปทำไม
ข้อมูลความหลากหลายทางชีวภาพเป็นข้อมูลที่พวกเราได้มาอย่างยากลำบาก กว่าจะเก็บตัวอย่าง กว่าจะระบุชนิดสิ่งมีชีวิตแต่ละตัวได้นั้น ล้วนใช้ทรัพยากรเงิน บุคลากร และ เวลาไปเป็นอย่างมาก ฉะนั้นแล้ว เราจะพยายามนำข้อมูลที่ได้มาได้นี้มาใช้ให้เกิดประโยชน์มากที่สุดเท่าที่จะทำได้ ผ่านการวิเคราะห์ข้อมูลในรูปแบบต่าง ๆ ที่จะสามารถนำไปใช้อธิบายรูปแบบความหลากหลายที่เราสำรวจได้อย่างเป็นระบบมากขึ้น
ใน workshop ครั้งนี้ เราจะให้ทุกท่านได้เริ่มเรียนรู้การใช้ R (สำหรับคนที่ไม่เคยใช้) ไปพร้อม ๆ กับการวิเคราะห์ความหลากหลายทางชีวภาพ ซึ่งเราจะใช้ข้อมูลพื้นฐานคือ
หน่วยการเก็บตัวอย่าง (Sampling Unit) เช่น แปลง สถานที่ ถาด บ่อ อำเภอ จังหวัด เป็นต้น
หน่วยอนุกรมวิธานที่สนใจ (Operational Taxonomic Units: OTUs) เช่น ชนิด สกุล วงศ์ เป็นต้น
จำนวนของสิ่งมีชีวิตนั้น ๆ ในแต่ละหน่วยการเก็บตัวอย่างและหน่วยอนุกรมวิธาน (Abundance) เช่น จำนวนตัว จำนวน read ร้อยละการปกคลุม เป็นต้น หากไม่ได้เก็บค่าจำนวนมา อาจใช้การพบหรือไม่พบ (absence/presence) แทนได้เช่นกัน
ถ้าเรามีข้อมูลเหล่านี้แล้ว ก็สามารถนำมาใช้วิเคราะห์ข้อมูลได้ในหลากหลายรูปแบบ ตามวัตถุประสงค์และคำถามวิจัยของแต่ละท่าน workshop นี้จะให้เครื่องมือแก่ท่านในการนำไปใช้วิเคราะห์ข้อมูลด้วยตนเองในอนาคต หรือ อย่างน้อยที่สุด ได้ทราบถึงแนวทางการวิเคราะห์ข้อมูลที่ได้จากการทำงานวิจัยด้านอนุกรมวิธาน
Platform ที่เราใช้อยู่นี้คือ RStudio ที่อยู่ออนไลน์ อยู่บน server ของบริษัท Posit ซึ่งทำให้เราสามารถใช้ R ที่ไหนก็ได้ผ่าน browser บนคอมพิวเตอร์ หรือ แม้กระทั่ง tablet แต่จะต้องมีอินเตอร์เน็ตอยู่ตลอดเวลา ถ้าท่านใดต้องการที่จะทำการวิเคราะห์จำนวนมาก อาจจะพิจารณาติดตั้งโปรแกรม RStudio ในเครื่องของท่านเองได้
0. Use the packages
เราจะเริ่มด้วยการนำเข้า packages เสริม ซึ่งก็คือ ชุดคำสั่งเพิ่มเติมจาก R พื้นฐาน (base R) ที่จะทำให้การทำงานของเราง่ายขึ้น โดยตัว โค้ดทั้งหมด จะอยู่ในแถบสีเทาด้านล่าง และ เราสามารถให้โค้ดทำงาน (Run Code) ได้โดยการกดปุ่ม Play สีเขียวที่มุมขวาบนของแต่ละกล่อง
หากมีตัวอักษรสีแดงจำนวนมากขึ้นมา ไม่ต้องตกใจนะครับ เป็นการแสดงว่า package เหล่านี้ได้รับการโหลดขึ้นมาเรียบร้อยแล้ว
library(vegan)## Loading required package: permute## Loading required package: lattice## This is vegan 2.6-4library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.2 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(picante)## Loading required package: ape
##
## Attaching package: 'ape'
##
## The following object is masked from 'package:dplyr':
##
## where
##
## Loading required package: nlme
##
## Attaching package: 'nlme'
##
## The following object is masked from 'package:dplyr':
##
## collapselibrary(broom)
library(pheatmap)
library(RColorBrewer)
library(ggfortify)1. Importing your diversity data
เราสามารถนำข้อมูลในรูปแบบไฟล์ csv (comma separated value) ซึ่งสามารถ Save
As ได้จาก Microsoft Excel หรือ โปรแกรมอื่น ๆ โดยใช้คำสั่ง read.csv
และนำผลที่อ่านได้ไปเก็บใน object ใหม่ ที่เราจะตั้งชื่อว่า
sample.sheet
จะเห็นว่าข้อมูลแบบนี้เป็นการจดแบบง่ายโดยให้คอลัมน์แรกเป็นหน่วยการเก็บข้อมูล
(plot ในที่นี้) และ คอลัมน์ถัดมาเป็น จำนวนตัว
(abundance) และตามมาด้วยชื่อชนิดในคอลัมน์ชื่อ sp
sample_sheet <- read.csv("sample_sheet.csv")
sample_sheet## plot abundance sp
## 1 site1 2 sp_a
## 2 site1 5 sp_b
## 3 site2 3 sp_b
## 4 site2 1 sp_c
## 5 site2 2 sp_d
## 6 site2 1 sp_e
## 7 site3 4 sp_e
## 8 site3 10 sp_fข้อมูลในลักษณะข้างต้น จะต้องมีการเปลี่ยนรูปแบบไปอยู่ในรูปแบบของ เมทริกซ์ (matrix)
เพื่อการคำนวณต่าง ๆ ต่อไป โดยสามารถใช้คำสั่งสำเร็จรูป sample2matrix
ได้ดังนี้
sample_mat <- sample2matrix(sample_sheet)
sample_mat## sp_a sp_b sp_c sp_d sp_e sp_f
## site1 2 5 0 0 0 0
## site2 0 3 1 2 1 0
## site3 0 0 0 0 4 10วิธีการนี้จะช่วยลดความผิดพลาดในการบันทึกข้อมูลในรูปแบบ matrix ตั้งแต่แรกที่ต้องมาเติม 0 ในช่องที่ไม่พบสิ่งมีชีวิตแบบเอง ยิ่งถ้าข้อมูลยิ่งเยอะ จะมีโอกาสผิดพลาดได้เยอะขึ้น
แต่ถ้าเรามีข้อมูลในรูปแบบ matrix อยู่แล้ว ก็สามารถนำเข้าในโปรแกรมได้เลยดังนี้ ข้อมูลด้านล่างนี้เป็นตัวอย่างข้อมูลสังคมพืชจากพื้นที่สันทราย (sand dune) แห่งหนึ่ง ที่มีการวางแปลงย่อย 20 แปลง และพบพืชทั้งหมด 30 ชนิด แหล่งที่มาข้อมูลสังคมพืชบนสันทราย
dune_comm <- read.csv("dune_comm.csv", row.names = 1)
dune_comm## Achimill Agrostol Airaprae Alopgeni Anthodor Bellpere Bromhord Chenalbu
## 1 1 0 0 0 0 0 0 0
## 2 3 0 0 2 0 3 4 0
## 3 0 4 0 7 0 2 0 0
## 4 0 8 0 2 0 2 3 0
## 5 2 0 0 0 4 2 2 0
## 6 2 0 0 0 3 0 0 0
## 7 2 0 0 0 2 0 2 0
## 8 0 4 0 5 0 0 0 0
## 9 0 3 0 3 0 0 0 0
## 10 4 0 0 0 4 2 4 0
## 11 0 0 0 0 0 0 0 0
## 12 0 4 0 8 0 0 0 0
## 13 0 5 0 5 0 0 0 1
## 14 0 4 0 0 0 0 0 0
## 15 0 4 0 0 0 0 0 0
## 16 0 7 0 4 0 0 0 0
## 17 2 0 2 0 4 0 0 0
## 18 0 0 0 0 0 2 0 0
## 19 0 0 3 0 4 0 0 0
## 20 0 5 0 0 0 0 0 0
## Cirsarve Comapalu Eleopalu Elymrepe Empenigr Hyporadi Juncarti Juncbufo
## 1 0 0 0 4 0 0 0 0
## 2 0 0 0 4 0 0 0 0
## 3 0 0 0 4 0 0 0 0
## 4 2 0 0 4 0 0 0 0
## 5 0 0 0 4 0 0 0 0
## 6 0 0 0 0 0 0 0 0
## 7 0 0 0 0 0 0 0 2
## 8 0 0 4 0 0 0 4 0
## 9 0 0 0 6 0 0 4 4
## 10 0 0 0 0 0 0 0 0
## 11 0 0 0 0 0 2 0 0
## 12 0 0 0 0 0 0 0 4
## 13 0 0 0 0 0 0 0 3
## 14 0 2 4 0 0 0 0 0
## 15 0 2 5 0 0 0 3 0
## 16 0 0 8 0 0 0 3 0
## 17 0 0 0 0 0 2 0 0
## 18 0 0 0 0 0 0 0 0
## 19 0 0 0 0 2 5 0 0
## 20 0 0 4 0 0 0 4 0
## Lolipere Planlanc Poaprat Poatriv Ranuflam Rumeacet Sagiproc Salirepe
## 1 7 0 4 2 0 0 0 0
## 2 5 0 4 7 0 0 0 0
## 3 6 0 5 6 0 0 0 0
## 4 5 0 4 5 0 0 5 0
## 5 2 5 2 6 0 5 0 0
## 6 6 5 3 4 0 6 0 0
## 7 6 5 4 5 0 3 0 0
## 8 4 0 4 4 2 0 2 0
## 9 2 0 4 5 0 2 2 0
## 10 6 3 4 4 0 0 0 0
## 11 7 3 4 0 0 0 2 0
## 12 0 0 0 4 0 2 4 0
## 13 0 0 2 9 2 0 2 0
## 14 0 0 0 0 2 0 0 0
## 15 0 0 0 0 2 0 0 0
## 16 0 0 0 2 2 0 0 0
## 17 0 2 1 0 0 0 0 0
## 18 2 3 3 0 0 0 0 3
## 19 0 0 0 0 0 0 3 3
## 20 0 0 0 0 4 0 0 5
## Scorautu Trifprat Trifrepe Vicilath Bracruta Callcusp
## 1 0 0 0 0 0 0
## 2 5 0 5 0 0 0
## 3 2 0 2 0 2 0
## 4 2 0 1 0 2 0
## 5 3 2 2 0 2 0
## 6 3 5 5 0 6 0
## 7 3 2 2 0 2 0
## 8 3 0 2 0 2 0
## 9 2 0 3 0 2 0
## 10 3 0 6 1 2 0
## 11 5 0 3 2 4 0
## 12 2 0 3 0 4 0
## 13 2 0 2 0 0 0
## 14 2 0 6 0 0 4
## 15 2 0 1 0 4 0
## 16 0 0 0 0 4 3
## 17 2 0 0 0 0 0
## 18 5 0 2 1 6 0
## 19 6 0 2 0 3 0
## 20 2 0 0 0 4 32. Data At-a-Glance
ก่อนที่จะเริ่มวิเคราะห์ใด ๆ เราจะเริ่มแสดงภาพความหลากหลายของเราเสียก่อนเพื่อให้เราเข้าใจข้อมูลของเรามากขึ้น โดยจะใช้การแสดงผลที่เรียกว่า heatmap ซึ่งจะมาในรูปแบบของตารางที่ในแต่ละช่องจะแสดงจำนวน (abundance) ของสิ่งมีชีวิตผ่านสีต่าง ๆ โดยสามารถวาดได้หลายรูปแบบดังนี้
1) Layout 1: look like your data matrix
ในแบบนี้จะเรียงลำดับตามข้อมูลที่เราระบุใน matrix ของเรา
pheatmap(sample_mat, cluster_rows = FALSE, cluster_cols = FALSE)
2) Layout 2: Reorder with their similarities + Cluster Diagram
ในแบบที่เป็น default ของ function จะมีการเรียงลำดับ plot และ species ตามความเหมือน (similarity) เพื่อจะทำให้เราเห็น pattern ภายในข้อมูลมากขึ้น
pheatmap(sample_mat)
เราสามารถปรับสีได้ดังนี้
my_color <- colorRampPalette(c("white", "brown", "red"))(50)
pheatmap(sample_mat, color = my_color)
หรือจะใช้ palette ที่มีอยู่แล้วโดยเลือกได้จาก palette ด้านล่างนี้ได้เลย
my_color <- colorRampPalette(brewer.pal(n = 9, name = "Greens"))(50)
pheatmap(sample_mat, color = my_color)
Challenge Set 1
ในแบบฝึกหัดนี้ เราจะลองนำเข้าข้อมูลใหม่ ซึ่งเป็นข้อมูลของนกในไต้หวันที่มีการเก็บข้อมูลจำนวน 50 แปลงตามระดับความสูงจากระดับน้ำทะเลที่แตกต่างกัน (https://www.davidzeleny.net/anadat-r/doku.php/en:data:ybirds) ให้แต่ละท่านลองดำเนินการต่อไปนี้
1.1 อ่านไฟล์ bird_comm.csv เข้าโปรแกรม จะพบนกทั้งหมด 59 ชนิด
จาก 50 แปลง.
1.2 วาด heatmap เพื่อแสดงภาพความหลากหลายในข้อมูลนี้
3. Types of Diversity
ในการพิจารณาความหลากหลายทางชีวภาพเราสามารถพิจารณาความหลากหลายได้ 3 ระดับดังต่อไปนี้
\(\alpha\) diversity (alpha diversity) = ความหลากหลายในแต่ละหน่วยเก็บข้อมูล บางครั้งเรียกว่า local diversity ซึ่งขอบเขตของหน่วยเก็บข้อมูลส่วนมากจะเป็นแปลงศึกษา หรือ พรมแดนธรรมชาติที่ผู้วิจัยกำหนด
\(\beta\) diversity (beta diversity) = ความแตกต่างขององค์ประกอบชนิดระหว่างหน่วยเก็บข้อมูล ซึ่งมีสูตรคำนวณความแตกต่างนี้จำนวนมาก ขึ้นกับวัตถุประสงค์ของผู้วิจัย
\(\gamma\) diversity (gamma diversity) = ความหลากหลายของทุกหน่วยเก็บข้อมูลร่วมกัน บางครั้งอาจจะเรียกว่า regional diversity หรือ ความหลากหลายระดับภูมิภาค ซึ่งจะขึ้นอยู่กับผู้วิจัยว่าจะกำหนดให้ “ภูมิภาค” ครอบคลุมแค่ไหน แต่ส่วนมากจะหมายถึงความหลากหลายโดยรวมของการศึกษาในครั้งนั้น ๆ
ใน workshop วันนี้ เราจะพาทุกท่านไปสำรวจความหลากหลายที่ระดับต่าง ๆ จากข้อมูล
dune_comm กันนะครับ
4. Gamma diversity
4.1) Total Number of Species
เราสามารถนับจำนวนชนิดของทุกแปลงรวมกัน โดยพิจารณาชนิดที่ซ้ำกันโดยเริ่มจาก บวกรวมค่า
abundance ของแต่ละชนิด (ที่อยู่ในแต่ละคอลัมน์กันก่อน ด้วยคำสั่ง
colSums
colSums(dune_comm)## Achimill Agrostol Airaprae Alopgeni Anthodor Bellpere Bromhord Chenalbu
## 16 48 5 36 21 13 15 1
## Cirsarve Comapalu Eleopalu Elymrepe Empenigr Hyporadi Juncarti Juncbufo
## 2 4 25 26 2 9 18 13
## Lolipere Planlanc Poaprat Poatriv Ranuflam Rumeacet Sagiproc Salirepe
## 58 26 48 63 14 18 20 11
## Scorautu Trifprat Trifrepe Vicilath Bracruta Callcusp
## 54 9 47 4 49 10จากข้อมูลนี้ เราก็จะนับชนิดที่มีมากกว่า 0 เป็นจำนวนชนิดที่พบ ด้วยคำสั่ง
length ก็จะพบว่า \(\gamma\)
diversity ของทั้งชุดข้อมูลนี้คือ 30 ชนิดนั่นเอง
length(colSums(dune_comm) > 0)## [1] 304.2) Species accumulation/Collecting Effort curve
Species Accumulation Curve
จะแสดงให้เราเห็นว่าการเก็บตัวอย่างของเราทั้งหมดนั้น
ใกล้เคียงจำนวนตัวอย่างที่ควรจะพบในพื้นที่ทั้งหมดนั้นหรือยัง
โดยโปรแกรมจะทำการวาดกราฟแสดงจำนวนชนิด สะสม
ที่เพิ่มขึ้นเมื่อเพิ่มจำนวนหน่วยเก็บข้อมูล (sampling unit)
และจะทำการสุ่มให้ด้วยว่าที่จำนวนแปลงเท่ากันนี้ หากเลือกคำนวณชนิดสะสมจากคนละแปลงกัน
จะได้จำนวนชนิดสะสมเท่าไหร่ ทั้งหมดนี้ดำเนินการด้วยคำสั่ง specaccum
(vegan package)
dune_specaccum <- specaccum(dune_comm)
dune_specaccum## Species Accumulation Curve
## Accumulation method: exact
## Call: specaccum(comm = dune_comm)
##
##
## Sites 1.000000 2.000000 3.000000 4.000000 5.000000 6.000000 7.00000
## Richness 9.850000 15.110526 18.510526 20.937461 22.754321 24.149587 25.23956
## sd 2.351064 1.876385 1.572271 1.446958 1.390159 1.353035 1.31648
##
## Sites 8.000000 9.000000 10.000000 11.000000 12.000000 13.000000 14.00000
## Richness 26.103548 26.798244 27.364962 27.833984 28.227546 28.562036 28.84964
## sd 1.274903 1.228201 1.176341 1.119344 1.056454 0.987409 0.91160
##
## Sites 15.000000 16.000000 17.00000 18.000000 19.000000 20
## Richness 29.099587 29.319092 29.51403 29.689474 29.850000 30
## sd 0.828689 0.738092 0.63339 0.513971 0.357071 0ซึ่งคำสั่ง specaccum
นี้จะให้ค่าเฉลี่ยและส่วนเบี่ยงเบนมาตรฐานของจำนวนชนิดสะสมในแต่ละแปลง ซึ่งหากเราให้คำสั่ง
plot ก็จะสามารถวิเคราะห์ข้อมูลนี้ได้ดีขึ้น
plot(dune_specaccum)
เส้นโค้งคือ accumulation curve ที่แสดงแนวโน้มจำนวนชนิดสะสม และ เส้นแนวตั้งคือจำนวนชนิดสะสมเฉลี่ย ± ส่วนเบี่ยงเบนมาตรฐาน (standard deviation)
เราสามารถปรับกราฟให้สวยงามขึ้นได้ดังนี้
plot(dune_specaccum, ci.type = "polygon", ci.col = "#57F3C7", ci.lty = 0)
4.3) Expected Number of Species (Species pool)
ในกรณีที่เราคิดว่าจำนวนชนิดที่พบยังไม่ใกล้เคียงที่จำนวนที่ควรจะพบทั้งหมดได้ เราสามารถคำนวณจำนวนที่ชนิดที่ควรจะพบ หรือ จำนวนชนิดทั้งหมด (species pool) ด้วย สูตรประมาณค่าความหลากหลาย (diversity estimators) ต่าง ๆ ดังนี้
specpool(dune_comm)## Species chao chao.se jack1 jack1.se jack2 boot boot.se n
## All 30 32.1375 3.242503 32.85 1.645448 33.84474 31.5404 1.153121 20รายละเอียดของการคำนวณแต่ละค่านั้นสามารถดูได้ที่
help("specpool")
Challenge Set 2
ในแบบฝึกหัดถัดไปเราจะใช้ข้อมูลนกจาก challenge แรก มาดู gamma diversity
2.1 วาด species accumulation curve เพื่อดูว่าการเก็บตัวอย่างไรสมบูรณ์หรือยัง
2.2 คำนวณจำนวนชนิดของนกทั้งหมดที่ควรพบในการศึกษาครั้งนี้
5. Alpha Diversity
5.1) Species Richness
วิธีการอธิบายความหลากหลายที่ง่ายที่สุดก็คือจำนวนชนิด หรือ species richness ซึ่งเราสามารถจำนวนของแต่ละแปลงได้ดังต่อไปนี้
dune_richness <- specnumber(dune_comm)
dune_richness## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 5 10 10 13 14 11 13 12 13 12 9 9 10 7 8 8 7 9 9 8ในหลายกรณี เราจะมีข้อมูลเพิ่มเติมของแต่ละแปลงด้วยว่ามีลักษณะอย่างไร เราจะเรียกข้อมูลเหล่านี้ว่า “environmental data” ซึ่งเราอ่านเข้ามาได้ดังนี้
dune_env <- read.csv("dune_env.csv", row.names = 1)
dune_env## A1 Moisture Management Use Manure
## 1 2.8 1 SF Haypastu 4
## 2 3.5 1 BF Haypastu 2
## 3 4.3 2 SF Haypastu 4
## 4 4.2 2 SF Haypastu 4
## 5 6.3 1 HF Hayfield 2
## 6 4.3 1 HF Haypastu 2
## 7 2.8 1 HF Pasture 3
## 8 4.2 5 HF Pasture 3
## 9 3.7 4 HF Hayfield 1
## 10 3.3 2 BF Hayfield 1
## 11 3.5 1 BF Pasture 1
## 12 5.8 4 SF Haypastu 2
## 13 6.0 5 SF Haypastu 3
## 14 9.3 5 NM Pasture 0
## 15 11.5 5 NM Haypastu 0
## 16 5.7 5 SF Pasture 3
## 17 4.0 2 NM Hayfield 0
## 18 4.6 1 NM Hayfield 0
## 19 3.7 5 NM Hayfield 0
## 20 3.5 5 NM Hayfield 0ในข้อมูลนี้จะแสดงข้อมูลของแต่ละแปลง ดังต่อไปนี้
A1 = ความหนาของดินชั้น A1
Moisture = ระดับความชื้นจากน้อยไปมาก (1 -> 5)
Management = การจัดการพื้นที่
BFBiological farmingHFHobby farmingNMNature Conservation Management และSFStandard Farming
Use = ระดับการใช้พื้นที่ :
Hayfield(ที่ปลูกฟาง รบกวนน้อยสุด) <Haypastu<Pasture(ทุ่งหญ้าเลี้ยงสัตว์ รบกวนมากสุด)Manure = ระดับการใส่ปุ๋ยคอกจากน้อยไปมาก (0->4)
เราสามารถคำนวณจำนวนชนิดให้แยกตามปัจจัยต่าง ๆ ที่เราเพิ่งนำเข้าได้ ดังนี้
specnumber(dune_comm, groups = dune_env$Management)## BF HF NM SF
## 16 21 21 215.2) Shannon’s Index
ในหลายครั้ง จำนวนชนิดก็ไม่ได้บอกทุกอย่างเกี่ยวกับพื้นที่ที่เราศึกษา เราลองพิจารณาภาพของ 2 พื้นที่ที่มีจำนวนชนิดเท่ากัน (richness = 4) ด้านล่าง
จะเห็นว่าด้านล่างจะดูหลากหลายกว่า ทั้ง ๆ ที่จำนวนชนิดเท่ากัน อันนี้เป็นเพราะว่า แต่ละชนิดในภาพล่าง มีจำนวนพอ ๆ กัน (even abundance) เลยทำให้ดูหลายกหลายกว่า เราสามารถนับจำนวนของแต่ละชนิดมาร่วมคำนวณเพื่อหาความหลากหลายด้วย เราจะเรียกค่านี้ว่า ดัชนีความหลากหลาย (diversity index) ซึ่งจะขอยกตัวอย่างการคำนวณ Shannon’s diversity index ซึ่งมีสูตรคำนวณดังนี้
\[H = -\sum_{i=1}^{n}(p_{i}~ln~p_{i})\] โดย \(p_{i}\) คือสัดส่วนจำนวนตัวของชนิดนั้น ๆ ในสังคมสิ่งมีชีวิตทั้งหมด เราสามารถให้โปรแกรมช่วยคำนวณได้ดังนี้
dune_shannon <- diversity(dune_comm)
dune_shannon## 1 2 3 4 5 6 7 8
## 1.440482 2.252516 2.193749 2.426779 2.544421 2.345946 2.471733 2.434898
## 9 10 11 12 13 14 15 16
## 2.493568 2.398613 2.106065 2.114495 2.099638 1.863680 1.979309 1.959795
## 17 18 19 20
## 1.876274 2.079387 2.134024 2.048270ดัชนี้นี้ไม่มีหน่วยและตัวเลขที่ได้อาจจะไม่สามารถเทียบเคียงได้กับการศึกษาอื่น ๆ โดยตรง เพราะค่าที่ได้จะขึ้นอยู่กับจำนวนชนิดทั้งหมดด้วย ดังนั้น ค่านี้จึงเหมาะกับการเปรียบเทียบภายในการศึกษาเดียวกัน
และก็สามารถคำนวณแยกกลุ่มตามปัจจัยที่สนใจได้เช่นเคย
diversity(dune_comm, groups = dune_env$Management)## BF HF NM SF
## 2.566756 2.870434 2.807265 2.6845205.3) Analyzing diversity data with other factors
หนึ่งในประโยชน์ของการวิเคราะห์ข้อมูลเหล่านี้ใน R คือ เราสามารถเชื่อมโยงข้อมูลที่คำนวณได้ใหม่กับข้อมูลอื่น ๆ โดยที่ไม่ต้อง copy-paste ใน Excel เลย ช่วยลดความผิดพลาดในการคำนวณได้ดีขึ้น
ตัวอย่างด้านล่างนี้เราจะเพิ่มคอลัมน์ในไฟล์ dune_env เดิม
เพื่อรองรับข้อมูล species richness และ Shannon’s index ดังนี้
dune_env$richness <- dune_richness
dune_env$shannon <- dune_shannon
dune_env## A1 Moisture Management Use Manure richness shannon
## 1 2.8 1 SF Haypastu 4 5 1.440482
## 2 3.5 1 BF Haypastu 2 10 2.252516
## 3 4.3 2 SF Haypastu 4 10 2.193749
## 4 4.2 2 SF Haypastu 4 13 2.426779
## 5 6.3 1 HF Hayfield 2 14 2.544421
## 6 4.3 1 HF Haypastu 2 11 2.345946
## 7 2.8 1 HF Pasture 3 13 2.471733
## 8 4.2 5 HF Pasture 3 12 2.434898
## 9 3.7 4 HF Hayfield 1 13 2.493568
## 10 3.3 2 BF Hayfield 1 12 2.398613
## 11 3.5 1 BF Pasture 1 9 2.106065
## 12 5.8 4 SF Haypastu 2 9 2.114495
## 13 6.0 5 SF Haypastu 3 10 2.099638
## 14 9.3 5 NM Pasture 0 7 1.863680
## 15 11.5 5 NM Haypastu 0 8 1.979309
## 16 5.7 5 SF Pasture 3 8 1.959795
## 17 4.0 2 NM Hayfield 0 7 1.876274
## 18 4.6 1 NM Hayfield 0 9 2.079387
## 19 3.7 5 NM Hayfield 0 9 2.134024
## 20 3.5 5 NM Hayfield 0 8 2.048270ข้อมูล dune_env
ที่อัพเดทใหม่นี้จะมีข้อมูลความหลากหลายเป็นคอลัมน์ต่อท้ายเรียบร้อยแล้ว
เราก็จะนำไปใช้ตอบคำถามผ่านการวาดกราฟหรือวิเคราะห์สถิติได้เลย ดังตัวอย่างต่อไปนี้
Example 1: ระดับปุ๋ยคอกที่ใส่เข้าไปมีผลกระทบต่อดัชนีความหลากหลายของพืชหรือไม่
dune_env %>%
ggplot(aes(x = Manure, y = shannon)) +
geom_point() +
geom_smooth(se = F)## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Example 2: จำนวนชนิดของพืชแตกต่างกันตามวิธีการจัดการพื้นที่หรือไม่
dune_env %>%
ggplot(aes(x = Management, y = richness)) +
geom_boxplot()
Challenge Set 3
จากข้อมูล bird_comm.csv ด้านบน
เราจะเพิ่มเติมข้อมูลสิ่งแวดล้อมของนกเหล่านี้ ในไฟล์ bird_env.csv
ซึ่งประกอบด้วยข้อมูลเพิ่มเติมดังต่อไปนี้
Veg_type= vegetation type (ชนิดของสังคมพืช)Elevation= ความสูงจากระดับทะเลTree_diversity= ความหลากหลายของไม้ต้น
3.1 อ่านไฟล์ bird_env.csv เข้าสู่โปรแกรม
3.2 คำนวณ species richness ของแต่ละ site
3.3 คำนวณ Shannon’s index ของแต่ละ site
3.4 เพิ่มคอลัมน์ species richness และ Shannon’s index เข้าไปใน bird_env
3.5 วาดกราฟเพื่อตรวจสอบว่า diversity ของนกเกี่ยวข้องกับ elevation หรือไม่
6. Beta Diversity
6.1) Calculating beta diversity
หัวข้อนี้เป็นหัวข้อที่อาจจะเข้าใจยากสักนิดหนึ่ง เพราะค่า beta diversity ไม่ใช่การนับเหมือน alpha หรือ gamma diversity แต่เป็นการคำนวณ ความแตกต่างขององค์ประกอบชนิด ระหว่าง sampling unit ซึ่งนักนิเวศวิทยาแต่ละท่านก็ได้นิยามไว้หลากหลายมาก (ตัวอย่างเพิ่มเติม ที่เปเปอร์ของ Koleff et al. 2003) เราจะเริ่มลองพิจารณาจากตัวอย่างง่าย ๆ ที่เราได้เริ่มไว้ด้านบนกันก่อน
sample_mat## sp_a sp_b sp_c sp_d sp_e sp_f
## site1 2 5 0 0 0 0
## site2 0 3 1 2 1 0
## site3 0 0 0 0 4 10จากข้อมูลของ sample_mat เราสามารถคำนวณ beta diversity
เพื่อหาความแตกต่างกันระหว่าง site โดยใช้วิธีของ Whittaker (คนคิดเรื่องนี้เป็นคนแรก
ใช้รหัสย่อ \(\beta_{w}\))
ซึ่งมีสูตรคำนวณคือ
\[ \beta_{w} = \frac{S}{\bar{\alpha}} - 1 \]
โดย S คือ จำนวนชนิดทั้งหมดที่พบใน 2 แปลง และ \(\bar{\alpha}\) คือจำนวนชนิดเฉลี่ยต่อแปลง เราสามารถใช้คำสั่งคำนวณได้ดังต่อไปนี้
sample_beta_w <- betadiver(sample_mat, method = "w")
sample_beta_w## site1 site2
## site2 0.6666667
## site3 1.0000000 0.6666667จากวิธีข้างต้น จะเห็นได้ว่า ความแตกต่างระหว่าง site1-site2 และ site2-site3 นั้นเท่ากัน แม้ว่าจำนวนต้นจะไม่เท่ากัน หากเราต้องการนำจำนวนต้นหรือจำนวนตัวมาใช้ประกอบการคำนวณ เราสามารถใช้วิธีการคำนวณ distance (ความแตกต่าง) มาแทนได้ ตัวอย่างเช่นเราสามารถใช้ Bray-Curtis Distance มาหาความแตกต่างระหว่าง 2 แปลงได้ โดยมีสูตรคำนวณดังนี้
\[ d_{jk} = \frac{\sum|x_{ij}-x_{ik}|}{\sum(x_{ij}+x_{ik})} \]
โดยที่ \(d_{jk}\) คือความแตกต่างขององค์ประกอบชนิดระหว่างแปลง j และ k โดยด้านบนเอาบนจำนวนตัวของแต่ละชนิดมาลบกันและบวกรวมทุกชนิด และ ตัวหาร เอาจำนวนตัวของแต่ละชนิดมาบวกกันและบวกรวมกันทุกชนิด จะเห็นว่าอะไรที่คำนวณเยอะ ๆ แบบนี้เหมาะยิ่งนักกับการให้คอมพิวเตอร์คำนวณให้เรา ตามโค้ดด้านล่างนี้
sample_beta_bray <- vegdist(sample_mat)
sample_beta_bray## site1 site2
## site2 0.5714286
## site3 1.0000000 0.90476196.2) Visualizing beta diversity
Option 1: Heatmap
หลังจากที่คำนวณ beta diversity ได้แล้ว เราก็สามารถนำค่าที่ได้มาวาดกราฟต่อได้ โดยแบบแรกจะใช้ heatmap เพื่อแสดงความเหมือนหรือความแตกต่างรายคู่ โดยเราจะทำให้เห็นทั้งสองวิธีคือ
1- วิธีของ Whittaker (พิจารณาเฉพาะการพบหรือไม่พบ)
autoplot(sample_beta_w) +
labs(title = "Whittaker's Beta Diversity")## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
2- วิธีของ Bray-Curtis (พิจารณาจำนวนต้น จำนวนตัวด้วย)
autoplot(sample_beta_bray) +
labs(title = "Bray-Curtis")## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
Option 2: Ordination
อีกวิธีที่นิยมกันมากในวงการคือ การนำข้อมูลที่ได้มาทำ ordination หรือ เรียงลำดับใหม่ในมิติที่ลดรูปลงทางคณิตศาสตร์ เพื่อให้สามารถทำความเข้าใจได้ง่ายขึ้น โดยหนึ่งในวิธีที่ได้รับความนิยมมากที่สุดก็คือ principal component analysis (PCA) ซึ่งจะขออนุญาตไม่ลงหลักการใน workshop แต่ละสาธิตวิธีการคำนวณและแสดงผลดังนี้
sample_pca <- prcomp(sample_beta_bray)
autoplot(sample_pca, label = T, point = F)
สิ่งที่ควรได้จากการอ่านกราฟนี้คือ
แต่ละจุดคือ sampling unit ของเรา
ระยะห่างระหว่างจุด คือ ความแตกต่างขององค์ประกอบชนิด ระหว่าง sampling unit = ยิ่งแตกต่าง จุดยิ่งห่างกัน
% ที่อยู่ตรงแกนคือ ปริมาณ variance ภายในข้อมูลที่ได้รับการอธิบายในการคำนวณนี้ ในกรณีนี้ 2 แกนรวมกันได้ 78.9 + 21.1 = 100% ก็คือ อธิบายได้ทุกอย่างในข้อมูลนี้เลย (เพราะข้อมูลมันไม่ใหญ่มากนั่นเอง)
เรามาลองดูข้อมูลของสังคมพืชสันทรายกันบ้าง
#Calculate beta diversity
dune_beta_bray <- vegdist(dune_comm)
#Perform PCA
dune_pca <- prcomp(dune_beta_bray)
#Visualize results
autoplot(dune_pca, label = T)
PCA ของข้อมูลนี้อธิบาย Variance ได้รวมกัน 47.89+24.15 = 72.04% ซึ่งถือว่าดีเลยทีเดียว จะเห็นว่าบางจุด (เช่น 12,13) อยู่ใกล้กัน แสดงว่าองค์ประกอบชนิดคล้ายคลึงมากกว่าจุดที่ไกลออกไป
นอกจากนี้เรายังสามารถนำข้อมูล environment มาช่วย highlight ผล PCA ให้อ่านง่ายและได้ข้อมูลมากขึ้นอีกด้วย เช่นเราสามารถระบายสีแต่ละจุดตามวิธีการจัดการพื้นที่ได้ดังนี้
autoplot(dune_pca, colour = "Management", data = dune_env, frame = T)
Challenge Set 4
จากข้อมูลนก bird_comm ด้านบน ให้ทุกท่านลองศึกษา beta diversity โดยดำเนินการต่อไปนี้
3.1 คำนวณ beta diversity ด้วยวิธีที่ท่านชอบ (Whittaker หรือ Bray-Curtis)
3.2 แสดงผล beta diversity ด้วย PCA ที่ให้สีตามชนิดของสังคมพืช
(Veg_type)
Final thoughts
ใน 3.5 ชั่วโมงที่ผ่านมา เราได้เรียนรู้ไอเดียและโค้ดต่าง ๆ ค่อนข้างเยอะ ในเวลาที่ไม่มากนัก หวังว่าอย่างน้อยที่สุดผู้เรียนทุกท่านจะได้เห็นแนวทางการทำงานทางด้านความหลากหลายต่อไปในอนาคตและจะได้มีโอกาสพัฒนาทักษะเพิ่มเติมได้ต่อไปครับ หากมีข้อสงสัยเพิ่มเติม สามารถติดต่อได้ที่ ekaphan.k@ku.th