THỰC HÀNH BIOINFORMATIC TRONG R

Cài đặt PACKAGES

BÀI 4: LÀM VIỆC VỚI LOÀI

Gọi packages taxize

library(taxize)

Trích xuất thông tin loài “Helinthus annus ssp. jigeri”

gnr_resolve("Helinthus annus ssp. jigeri", data_source_ids=3)
## # A tibble: 1 x 5
##   user_supplied_name  submitted_name   matched_name       data_source_tit… score
## * <chr>               <chr>            <chr>              <chr>            <dbl>
## 1 Helinthus annus ss… Helinthus annus… Helianthus annuus… ITIS              0.54

Lấy thông tin loài “Helinthus annus ssp. jigeri” trên CSDL thông tin đơn vị phân loại itis

get_tsn("Helianthus annuus ssp. jaegeri")
## ══  1 queries  ═══════════════
## 
## Retrieving data for taxon 'Helianthus annuus ssp. jaegeri'
## Warning: Unknown or uninitialised column: 'itisTerms'.
## ✖  Not Found:  Helianthus annuus ssp. jaegeri
## ══  Results  ═════════════════
## 
## ● Total: 1 
## ● Found: 0 
## ● Not Found: 1
## [1] NA
## attr(,"class")
## [1] "tsn"
## attr(,"match")
## [1] "not found"
## attr(,"multiple_matches")
## [1] FALSE
## attr(,"pattern_match")
## [1] FALSE

Kiểm tra tên taxon đã được chấp nhận

itis_acceptname(525928)
## Warning: Unknown or uninitialised column: 'acceptedNames'.
##   submittedtsn acceptedname acceptedtsn author
## 1       525928           NA      525928     NA

Trích xuất thông tin theo ID trên database itis

classification(36616, db="itis")
## $`36616`
## [1] NA
## 
## attr(,"class")
## [1] "classification"
## attr(,"db")
## [1] "itis"

Truy xuất thông tin trên cơ sở dữ liệu đa dạng sinh học toàn cầu GBIF

get_gbifid("Shorea parvifolia")
## ══  1 queries  ═══════════════
## 
## Retrieving data for taxon 'Shorea parvifolia'
## ✔  Found:  Shorea parvifolia
## ══  Results  ═════════════════
## 
## ● Total: 1 
## ● Found: 1 
## ● Not Found: 0
## [1] "4096343"
## attr(,"class")
## [1] "gbifid"
## attr(,"match")
## [1] "found"
## attr(,"multiple_matches")
## [1] TRUE
## attr(,"pattern_match")
## [1] TRUE
## attr(,"uri")
## [1] "https://www.gbif.org/species/4096343"

BÀI 5: PHÂN TÍCH CÁC ĐẶC ĐIỂM HÌNH THÁI VỚI R

  • Nếu chúng ta có một ma trận các đặc điểm hình thái của một hay nhiều taxa, chúng ta có thể phân tích mối tương quan (Kiểu ma trận khoảng cách).

  • Sử dụng ma trận khoảng cách đơn giản:

Ví dụ: Làm việc với dữ liệu iris có sẵn trong R

data(iris)

Thống kê số lượng cá thể mỗi loài

table(iris$Species)
## 
##     setosa versicolor  virginica 
##         50         50         50

Quy định màu sắc theo loài

iris$color[iris$Species == 'virginica'] <- "red"
iris$color[iris$Species == 'versicolor'] <- "blue"
iris$color[iris$Species == 'setosa'] <- "green"

Phân tích mối tương quan theo các đặc điểm về chiều dài, chiều rộng của cánh hoa và đài hoa

pairs(iris[,c(-5,-6)], col = iris$color, lower.panel = NULL, pch=19)

Phân tích mối quan hệ phát sinh dựa vào ma trận khoảng cách đặc điểm hình thái

plot(hclust(dist(iris[,c(-5,-6)], method="manhattan")), labels=iris$Species)

Sử dụng ma trận khoảng cách phức tạp bằng gói packages cluster

library(cluster)
data(flower)
head(flower) # Xem các dòng đầu tiên của dữ liệu flower
##   V1 V2 V3 V4 V5 V6  V7 V8
## 1  0  1  1  4  3 15  25 15
## 2  1  0  0  2  1  3 150 50
## 3  0  1  0  3  3  1 150 50
## 4  0  0  1  4  2 16 125 50
## 5  0  1  0  5  2  2  20 15
## 6  0  1  0  4  3 12  50 40
str(flower) # Xem phân loại biến
## 'data.frame':    18 obs. of  8 variables:
##  $ V1: Factor w/ 2 levels "0","1": 1 2 1 1 1 1 1 1 2 2 ...
##  $ V2: Factor w/ 2 levels "0","1": 2 1 2 1 2 2 1 1 2 2 ...
##  $ V3: Factor w/ 2 levels "0","1": 2 1 1 2 1 1 1 2 1 1 ...
##  $ V4: Factor w/ 5 levels "1","2","3","4",..: 4 2 3 4 5 4 4 2 3 5 ...
##  $ V5: Ord.factor w/ 3 levels "1"<"2"<"3": 3 1 3 2 2 3 3 2 1 2 ...
##  $ V6: Ord.factor w/ 18 levels "1"<"2"<"3"<"4"<..: 15 3 1 16 2 12 13 7 4 14 ...
##  $ V7: num  25 150 150 125 20 50 40 100 25 100 ...
##  $ V8: num  15 50 50 50 15 40 20 15 15 60 ...
lapply(flower, class)  # Phân loại biến
## $V1
## [1] "factor"
## 
## $V2
## [1] "factor"
## 
## $V3
## [1] "factor"
## 
## $V4
## [1] "factor"
## 
## $V5
## [1] "ordered" "factor" 
## 
## $V6
## [1] "ordered" "factor" 
## 
## $V7
## [1] "numeric"
## 
## $V8
## [1] "numeric"
summary(daisy(flower))
## 153 dissimilarities, summarized :
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1418  0.3904  0.4829  0.4865  0.5865  0.8875 
## Metric :  mixed ;  Types = N, N, N, N, O, O, I, I 
## Number of objects : 18
plot(hclust(daisy(flower)))

summary(daisy(flower, type = list(asymm = c("V1", "V3"), symm= 2,
                ordratio= 7, logratio= 8)))
## 153 dissimilarities, summarized :
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1649  0.4378  0.5350  0.5288  0.6318  0.8972 
## Metric :  mixed ;  Types = A, S, A, N, O, O, T, I 
## Number of objects : 18
flower.dist <- daisy(flower, type = list(asymm = c("V1", "V3"), symm= 2,
                            ordratio= 7, logratio= 8))
plot(hclust(flower.dist))

plot(hclust(flower.dist, method="average"))

Xây dựng mối quan hệ phát sinh dựa vào ma trận các đặc điểm hình thái

Ví dụ: Thiết lặp ma trận các đặc điểm hình thái của: Cá mập (shark), Cá vàng (goldfish), cá vây Tây (coelecanth),

  • Cự đà (iguana), Chim cổ đỏ (Robin), Con người (human), Vượn cáo (lemur), Dơi (bat), Lợn (pig), Cá voi (whale), bò (cow)
  • Với ma trận các đặc điểm: Vẩy, vây, cánh, lông, máu, móng guốc, chân. đặc điểm nào có = 1, không có = 0

Các bước thực hiện như sau:

Thiết lập các ma trận đặc điểm của từng loài

shark  <- c(1, 1, 0, 0, 0, 0, 0)

Tương tự các loài khác

goldfish   <- c(1, 1, 0, 0, 0, 0, 0)
coelecanth <- c(1, 1, 0, 0, 0, 0, 0)
iguana     <- c(1, 0, 0, 0, 0, 0, 1)
robin      <- c(0, 0, 1, 0, 1, 0, 1)
human      <- c(0, 0, 0, 1, 1, 0, 1)
lemur      <- c(0, 0, 0, 1, 1, 0, 1)
bat        <- c(0, 0, 1, 1, 1, 0, 1)
pig        <- c(0, 0, 0, 1, 1, 1, 1)
whale      <- c(0, 1, 0, 0, 1, 0, 0)
cow        <- c(0, 0, 0, 1, 1, 1, 1)

Tạo dataframe mới gồm các loài có đặc điểm trên

animals <- rbind(shark, goldfish, coelecanth, iguana, 
                 robin, human, lemur, bat, pig, whale, cow)
colnames(animals) <- c("scales","fins","wings","hair",
                       "warm-blood","hooves","legs")
animals <- as.data.frame(animals)
animals
##            scales fins wings hair warm-blood hooves legs
## shark           1    1     0    0          0      0    0
## goldfish        1    1     0    0          0      0    0
## coelecanth      1    1     0    0          0      0    0
## iguana          1    0     0    0          0      0    1
## robin           0    0     1    0          1      0    1
## human           0    0     0    1          1      0    1
## lemur           0    0     0    1          1      0    1
## bat             0    0     1    1          1      0    1
## pig             0    0     0    1          1      1    1
## whale           0    1     0    0          1      0    0
## cow             0    0     0    1          1      1    1

Tìm một loài bất kỳ, vd: có đặc điểm có vẩy và có chân

row.names(animals[animals$scales == 1 & animals$legs == 1,])
## [1] "iguana"

Bài tập: Phân tích mối tương quan giữa các loài trên bằng xây dựng cây phát sinh dựa trên các đặc điểm hình thái trên?

BÀI 6: LÀM VIỆC VỚI TRÌNH TỰ DNA TRÊN CÁC CSDL

Dowload trình tự nucleotide từ các cơ sở dữ liệu online

Gọi packages “rentrez” của NCBI

library(rentrez)

Download trình tự

sea <- entrez_search(db="nuccore", term="Lithocarpus dahuoaiensis[Organism] AND matK[Gene]") 
sea
## Entrez search result with 2 hits (object contains 2 IDs and no web_history object)
##  Search term (as translated):  "Lithocarpus dahuoaiensis"[Organism] AND matK[Gene ...
sea$ids
## [1] "1240296417" "1240296411"
sp1.fasta <- entrez_fetch(db="nuccore", sea$ids, rettype="fasta")
sp1.fasta
## [1] ">LC318551.1 Lithocarpus dahuoaiensis chloroplast matK gene for maturase K, partial cds, specimen_voucher: FU<JPN_:V3194\nATTCATTCAATATTTCCCTTTTTAGAGGAAAAGTTCCCACATTTAAATTATTCGGCGGATATACTAATAC\nCCTACCCTGCCCATCTGGAAATATTGGTTCAAACCCTTCGTTACCGGGTGAAAGATGCTTCTTATTTGCA\nTTTATTGCGGTTCTTTCTTCATGAGTATTCTAATTGTAACAGTCTTATTATTACAAATAAATCTATTTCC\nATTTTTTCAAAAAGTAATCCGAGATTCTTTTTATTCCTATATAATTCTTATATATGTGAATACGAATCCA\nTCTTCCTTTTTCTCCGTAACCAATCTTCTCATTTACGATTAACATCTTCCGGAGTCCTTTTTGAACGACT\nCTGTTTATATAGAAAAATAGAACATTTTGCCGAAGTCTTTGCTAATGATTTTCCGGTCATCCCATGCTTT\nCTCAAGGATCCTTTCATGCATTATGTTAGATATCAAGGAAAATCAATTCTGGCTTCCAAAGACACACCTC\nTTCTAATGAATAAATGGAAATCTTACCTTGTCAATTTATGGCAATGTCATTTTGATGTATGGTCTCACGC\nGGCAAGTATCCGTATAAACCAATTATCCAAGCATTCCCTCGATTTTTTGAGTTACTTTTCAAGTGTTCGA\nCGAAATCCTGCAGTGGTGCGGAATCAAATGCTAGAAAGTTCATTTCTACTAAATAATGCTCCCAATAAAC\nTCGATACAATAGTTCCAATTATTCCTCTGATTGGATCATTGGCTAAAGCGAAATTTTGTAACGCAGTAGG\nGCATCCAATTAGTAAGCTGACTCGGGCCGATTTATCGGATTTTGAGATTATCAATCGATTTTTGCATATA\nTGCAGAAATCTTTCTCATTATTACAGCGGATCCTCAAAAAAAAAGAATATGTATCGAATAAAATATATAC\nTTCGAC\n\n>LC318548.1 Lithocarpus dahuoaiensis chloroplast matK gene for maturase K, partial cds, specimen_voucher: FU<JPN_:V5404\nAGAAAAGTTCCCACATTTAAATATTCGGCGGATATACTAATACCCTACCCTGCCCATCTGGAAATATTGG\nTTCAAACCCTTCGTTACCGGGTGAAGGATGCTTCTTATTTGCATTTATTGCGGTTCTTTCTTCATGAGTA\nTTCTAATTGTAACAGTCTTATTATTACAAATAAATCTATTTCCATTTTTTCAAAAAGTAATCCGAGATTC\nTTTTTATTCCTATATAATTCTTATATATGTGAATACGAATCCATCTTCCTTTTTCTCCGTAACCAATCTT\nCTCATTTACGATTAACATCTTCCGGAGTCCTTTTTGAACGACTCTGTTTATATAGAAAAATAGAACATTT\nTGCCGAAGTCTTTGCTAATGATTTTCCGGTCATCCCATGCTTTCTCAAGGATCCTTTCATGCATTATGTT\nAGATATCAAGGAAAATCAATTCTGGCTTCCAAAGACACACCTCTTCTAATGAATAAATGGAAATCTTACC\nTTGTCAATTTATGGCAATGTCATTTTGATGTATGGTCTCACGCGGCAAGTATCCGTATAAACCAATTATC\nCAAGCATTCCCTCGATTTTTTGAGTTACTTTTCAAGTGTTCGACGAAATCCTGCAGTGGTGCGGAATCAA\nATGCTAGAAAGTTCATTTCTACTAAATAATGCTCCCAATAAACTCGATACAATAGTTCCAATTATTCCTC\nTGATTGGATCATTGGCTAAAGCGAAATTTTGTAACGCAGTAGGGCATCCAATTAGTAAGCTGACTCGGGC\nCGATTTATCGGATTTTGAGATTATCAATCGATTTTTGCATATATGCAGAATCTTTCTCATA\n\n"
write(sp1.fasta, "Sp_matK_1.fasta")

Bài tập: Thực hiện load trình tự ITS, matK và rbcL của loài Lithocarpus vuquangensis

Dùng packages seqinr

library(seqinr)

Chọn cơ sở dữ liệu

choosebank("genbank")

Download trình tự

sp <- query("sp", "sp=Lithocarpus dahuoaiensis AND K=matK")
sp$req
## [[1]]
##       name     length      frame     ncbicg 
## "LC318548"      "831"        "0"       "11" 
## 
## [[2]]
##       name     length      frame     ncbicg 
## "LC318551"      "916"        "0"       "11"

Ghi thành file fasta

write.fasta(getSequence(sp), getName(sp), file.out="Sp_matK_2.fasta")

Đọc trình tự

Cài đặt và dùng packages ape install.packages(“ape”)

library(ape)
## 
## Attaching package: 'ape'
## The following objects are masked from 'package:seqinr':
## 
##     as.alignment, consensus

Đọc trình tự bằng lệnh read.dna

dna1 <- read.dna("Sp_matK_1.fasta", format="fasta")
dna1
## 2 DNA sequences in binary format stored in a list.
## 
## Mean sequence length: 873.5 
##    Shortest sequence: 831 
##     Longest sequence: 916 
## 
## Labels:
## LC318551.1 Lithocarpus dahuoaiensis chloroplast matK gene fo...
## LC318548.1 Lithocarpus dahuoaiensis chloroplast matK gene fo...
## 
## Base composition:
##     a     c     g     t 
## 0.294 0.199 0.149 0.358 
## (Total: 1.75 kb)
dna2 <- read.dna("Sp_matK_2.fasta", format="fasta")
dna2
## 2 DNA sequences in binary format stored in a list.
## 
## Mean sequence length: 873.5 
##    Shortest sequence: 831 
##     Longest sequence: 916 
## 
## Labels:
## LC318548
## LC318551
## 
## Base composition:
##     a     c     g     t 
## 0.294 0.199 0.149 0.358 
## (Total: 1.75 kb)
class(dna1)
## [1] "DNAbin"

—————————– END ——————————–

Reference: Modified from Camwebb lectures at Kyushu University, 2017