Sử dụng tool Bacollite xác định các mẫu động vật khảo cổ tự động từ quang phổ MALDI-ToF (Phần 2)
Bài đăng này đã không được cập nhật trong 3 năm
1. Lời mở đầu
Ở bài viết trước, mình đã cùng tìm hiểu được mục đích nghiên cứu của bài báo và nắm được cách hoạt động của công cụ quang phổ MALDI-ToF, các nhà nghiên cứu đã phát triển một tool cho phép thực hiện việc đối chiếu, căn chỉnh tự động giữa mẫu lý thuyết và mẫu được MALDI-ToF phân tích. Ở bài viết lần này, chúng ta sẽ cùng cài đặt và thực hiện một căn chỉnh mô phỏng theo hướng dẫn của bài báo nhé. Bắt đầu thôi =))
Link tool mình để ở đây: Bacollite - https://github.com/bioarch-sjh/bacollite
2. Cài đặt R script
Vì công cụ này được viết bởi ngôn ngữ lập trình R script nên đầu tiên chúng ta cùng cài đặt nó đã, mở Terminal lên và gõ như sau:
sudo add-apt-repository 'deb https://cloud.r-project.org/bin/linux/ubuntu bionic-cran40/'
sudo apt-get update
sudo apt-get install r-base
sudo -i R
Bạn sẽ nhận được output tương tự như sau:
3. Giới thiệu căn chỉnh mô phỏng
Giờ mình sẽ trình bày cách tải ba bản sao của một mẫu từ tệp văn bản và căn chỉnh chúng với một số peptit từ bốn loài khác nhau bằng cách sử dụng gói bacollite. Quá trình cơ bản của việc sắp xếp các pic đối với peptit sẽ được mô tả trước tiên, sau đó là chi tiết về quá trình phân loại quang phổ dựa trên thông tin peptide sẽ được đưa ra. Tóm lại, quy trình bao gồm các bước sau:
- Tải tập dữ liệu MALDI
- Nạp peptide lý thuyết cho một số loài
- Căn chỉnh dữ liệu MALDI với các peptit
- Cho điểm cho các peptit từ mỗi loài, để suy ra ID loài
- Tạo biểu diễn đồ thị của điểm
Gói Bacollite yêu cầu có devtools package nên chúng ta cần install package đó bằng 2 câu lệnh sau:
install.packages( "devtools" )
devtools::install_github("bioarch-sjh/bacollite", force=T)
Sau khi cài đặt, package có thể được tải theo cách thông thường:
require(bacollite)
Vậy là chúng ta đã chuẩn bị xong môi trường cần thiết để thực hiện mô phỏng này. Tiếp theo sẽ là từng bước theo quy trình mình đã liệt kê ở trên.
4. Cách tải dữ liệu
4.1 Thiết lập thư mục dữ liệu
Gói Bacollite cung cấp một số tệp dữ liệu mẫu có thể được sử dụng trong hướng dẫn này. Chúng ta cần một số mã lệnh R để tìm vị trí của chúng trên hệ thống của chúng ta. Để R cho chúng ta biết thư mục chứa dữ liệu, chúng ta sử dụng lệnh sau:
fpath <- system.file("extdata",package="bacollite")
4.2 Tải tệp dữ liệu cho ba bản sao
Tệp dữ liệu MALDI mà chúng ta sẽ sử dụng:
20150326_SF_YG65_94_I11.txt
20150326_SF_YG65_94_I14.txt
20150326_SF_YG65_94_I17.txt
Bạn có thể thấy rằng phần cuối cùng của tên tệp trước “.txt” là tham chiếu tại chỗ và có ba điểm: I11, I14 và I17. Chúng ta cần thông tin này để sử dụng hàm load.sample để tải ba mẫu thành R, như thế này:
library(bacollite)
// load package
froot = sprintf("%s/20150326_SF_YG65_94_",fpath)
// tạo văn bản chung cho tệp đường dẫn của cả ba mẫu, dựa trên chuỗi văn bản fpath mà chúng ta đã tạo ở trên.
sample <- load.sample(froot = froot,spots = c("I11","I14","I17"),name="folio 42")
// tải dữ liệu từ các tệp vào một đối tượng dữ liệu R được gọi là sample.
Hãy xem cấu trúc của đối tượng sample sử dụng hàm str() của R:
Điều này cho chúng ta biết rằng sample là một list gồm 5 phần dữ liệu:
- name: tên của dữ liệu, ở đây gọi là "folio 42".
- spot: 3 spot code tương ứng với 3 mẫu.
- s1 - dữ liệu mẫu cho tệp đầu tiên (I11)
- s3 - dữ liệu mẫu cho tệp thứ hai (I14)
- s3 - dữ liệu mẫu cho tệp thứ ba (I17)
Mỗi mục trong ba mục dữ liệu mẫu trong danh sách là một khung dữ liệu chứa dữ liệu phổ MALDI mà chúng ta sẽ sử dụng để phân tích.
4.3 Cách nạp peptit
Có một số bộ trình tự peptit đi kèm với bacollite. Gói này đi kèm với một số bộ peptide mà chúng ta có thể sử dụng để phân biệt giữa các mẫu từ cừu, bò, dê và hươu.
Các peptit trong dm_sheep là:
Trong dm_cow là:
Trong dm_goat là:
Trong dm_deer là:
Ở đây mình sẽ không đề cập đến về cách tạo bộ peptit từ trình tự axit amin thô.
4.4 Chạy căn chỉnh
OK, bây giờ chúng ta đã có mẫu và một số peptit, chúng ta có thể căn chỉnh bằng cách sử dụng hàm ms_fit(). Tiến hành cho mỗi bộ điểm đánh dấu. Nếu biến doplot được đặt thành TRUE thì đồ thị sẽ được tạo cho mỗi peptit:
par(mfrow=c(4,5))
gauss <- 0.2
sheep_fit <- ms_fit(peptides = dm_sheep, sample = sample, doplot=T, force=T, gauss=0.2)
cow_fit <- ms_fit(peptides = dm_cow, sample = sample, doplot=T, force=T, gauss=0.2)
goat_fit <- ms_fit(peptides = dm_goat, sample = sample, doplot=T, force=T, gauss=0.2)
deer_fit <- ms_fit(peptides = dm_deer, sample = sample, doplot=T, force=T, gauss=0.2)
Sau khi chạy xong đống lệnh này, ui chao, nó xuất ra một đống đồ thị cũng khá bắt mắt 
Hãy nhìn vào 1 đồ thị trong số đó, mình sẽ giải thích qua như sau:
- Trục x hiển thị phạm vi khối lượng đang xem xét
- Trục y là thang điểm từ 0 đến 1
- Mỗi mẫu được vẽ bằng một vạch màu và một vạch xám. Các mẫu căn chỉnh 1,2 và 3 có màu đỏ, xanh lục và xanh lam. Tại đó các vị trí khối lượng ban đầu (chưa được dịch chuyển) cũng được biểu thị bằng đường màu xám. Điều này rất hữu ích cho các căn chỉnh không tốt
- Sự phân bố đồng vị của các peptit mục tiêu được thể hiện dưới dạng một chuỗi 5 "đỉnh nhọn" như hình trên. Chúng được chia tỷ lệ bằng giá trị tương đối của đồng vị phổ biến nhất. Đối với khối lượng thấp hơn, đây thường là đồng vị đầu tiên nhưng đối với khối lượng cao hơn nó thường là thứ hai.
Tại thời điểm này, chúng ta đã có một bộ dữ liệu tương quan cho mỗi mẫu lặp lại với các peptit cho mỗi loài.Bước tiếp theo là sử dụng dữ liệu này để thực hiện phân loại mẫu.
4.5 Thực hiện phân loại
Chúng ta sẽ thực hiện phân loại bằng cách chuyển các đối tượng dữ liệu mà chúng ta đã tạo với ms_fit vào hàm phân loại cor_id.
Đầu tiên, chúng ta cần tạo một list các căn chỉnh và một vectơ tên cho ID loài như sau:
cordata <- list()
cordata[[1]]<-sheep_fit
cordata[[2]]<-cow_fit
cordata[[3]]<-goat_fit
cordata[[4]]<-deer_fit
corlab <- c("Sheep","Cow","Goat","Deer")
Giai đoạn đầu tiên trong quá trình phân loại là xác định xem liệu mối tương quan giữa mẫu và peptit có 'hit'. Thay vì sử dụng một ngưỡng duy nhất cho mối tương quan, chúng ta cần sử dụng một loạt các ngưỡng (được lưu trữ trong biến corlim bên dưới) và xem liệu căn chỉnh có nằm trên mỗi giá trị trong phạm vi đó hay không. Đoạn mã dưới đây thực hiện quá trình này:
# Set up
corlim = seq(0,1,0.05)
scores <- vector(length = length(cordata))
scores[] <- 0
laglim <- 0.6
# Massage the raw cordata into a form we can work with:
cld <- list()
for(cc in 1:length(cordata)) {
cld[[cc]] <- corlim_data(cordata[[cc]],laglim)
}
Chúng ta có thể xem xét cấu trúc của đối tượng cld bằng hàm str(), điều này cho chúng ta thấy kết quả về những gì đoạn mã trên đạt được:
Như bạn có thể thấy, chúng ta có một danh sách gồm 4 data frame, một frame cho mỗi loài. Mỗi data frame có một bộ giá trị cho mỗi ngưỡng tương quan đang được xem xét. Biến cl là ngưỡng tương quan, nh là giá trị cung cấp số lần truy cập cho mỗi ngưỡng và giá trị sc là tổng số ion cho các đỉnh được phân loại là một lần truy cập (chúng ta sẽ không sử dụng biến sc cho phân tích này).
Phần chính của phân loại là kết hợp dữ liệu này thành điểm số cho từng ID loài ứng viên, tương ứng với phương trình (2) trong bài báo mà mình đã giới thiệu ở bài viết trước. Điểm được tích lũy cho mỗi giá trị ngưỡng tương quan. Đối với mỗi ngưỡng tương quan và số lần 'truy cập' cho ngưỡng đó đối với mỗi loài:
#initialise the scores for each sample
for(ss in 1:length(cld))
cld[[ss]]$cumscore = 0
#for each correlation threshold
for(cl in 1:length(corlim)) {
#for each species
for(ss in 1:length(cld)) {
#get the number of hits for this species
nh <- cld[[ss]]$nh[cl]
maxonh<-0
#get the max number of hits for the other candidate species
for(tt in 1:length(cld)) {
if(ss != tt) {
maxonh <- max(maxonh,cld[[tt]]$nh[cl])
}
}
#add to the score if nh is greater than maxonh
if(nh > maxonh) {
cld[[ss]]$cumscore[cl] <- (nh-maxonh)*corlim[cl]
}
}
}
Với những giá trị này, chúng ta có thể báo cáo điểm số và tích lũy điểm số cho từng loài vào một cấu trúc để giữ kết quả:
result <- data.frame("id" = corlab, "score" = 0)
for(ss in 1:length(cld)){
message(sprintf("Score for species %d (%s) = %f" ,ss,corlab[ss],sum(cld[[ss]]$cumscore)))
result$score[ss] <- sum(cld[[ss]]$cumscore)
}
#> Score for species 1 (Sheep) = 0.000000
#> Score for species 2 (Cow) = 72.700000
#> Score for species 3 (Goat) = 0.000000
#> Score for species 4 (Deer) = 0.000000
Điều này cho điểm trước mỗi loài. Trong trường hợp này, mẫu được phân loại rõ ràng là bò. Đôi khi, khi mẫu bị nhiễu hoặc bị ô nhiễm, có thể có các điểm khác 0 cho nhiều hơn một loài. Đây là hữu ích , vì nó cho chúng ta dấu hiệu về việc phân loại tự động có thực sự mạnh hay không.
4.6 Tạo biểu diễn đồ thị của điểm
Giai đoạn cuối cùng trong phân tích là tạo ra một biểu diễn đồ thị của quá trình tính điểm:
title = sprintf("Sample '%s': manual ID: '%s'; Calc ID: '%s'\nscores ",sample$name,"unknown",result$id[r
for(ss in 1:nrow(result)){
title = sprintf("%s %s = %0.3f",title,result$id[ss],result$score[ss])
}
par(mar = c(4,4,5,4))
plot(NA,xlab="Correlation Threshold",ylab = "Number of Hits",ylim=c(0,15),xlim = c(0,1), main = title)
points(x=corlim,y=cld[[1]]$nh,col="#e2ba5e")
points(x=corlim,y=cld[[2]]$nh,col="#3597c6")
points(x=corlim,y=cld[[3]]$nh,col="#e3645f")
points(x=corlim,y=cld[[4]]$nh,col="#793787")
lines(x=corlim,y=cld[[1]]$nh,col="#e2ba5e")
lines(x=corlim,y=cld[[2]]$nh,col="#3597c6")
lines(x=corlim,y=cld[[3]]$nh,col="#e3645f")
lines(x=corlim,y=cld[[4]]$nh,col="#793787")
legend("topright",legend = corlab,col=c("#e2ba5e","#3597c6","#e3645f","#793787"),lty = 1,pch=1)
Đồ thị thu được như sau:
Cùng mình phân tích đồ thị này 1 chút:
- Có năm peptit cho mỗi loài và ba lần sao chép, vì vậy số lần truy cập tối đa vào một mẫu có thể là 3 x 5 = 15.
- Cow có số lần truy cập vào 14 trong số 15 peptit lên đến ngưỡng tương quan là 0,7, vì vậy đây là một phân loại đơn giản.
- Cừu có sáu peptit chung với bò - đó là lý do tại sao số lần truy cập vẫn ở mức 5 cho đến cùng một ngưỡng phân loại.
- Ngược lại, dê và hươu không có peptit chung với bò, do đó số lần truy cập giảm dần với ngưỡng tương quan thấp hơn.
Kết luận
Vậy là chúng ta cũng đã tìm hiểu xong và code thử 1 mô phỏng để thấy được việc phân tích tự động của tool khá hữu ích. Trong khi việc phân tích thủ công cần những chuyên gia thực sự và khá tốn thời gian căn chỉnh để thu được kết quả chính xác nhất thì giờ chỉ cần thông qua công cụ Bacollite là đã có thể xác định được rồi.
Có thể bài viết khá nhàm chán và khó hiểu cho những ai không thích sinh học :v thực ra là mình cũng vậy, nhưng nhân cơ hội kì này mình được học môn Tin sinh học và được giao tìm hiểu về một bài báo nghiên cứu tin sinh nên thấy nó khá cuốn :v
Mình xin phép dừng lại bài viết ở đây vì nó quá dài rồi 
, hẹn các bạn vào những bài viết đúng chuyên môn hơn vào lần sau. Xin cảm ơn các bạn (bow)
All rights reserved
