Fungi database에는 UNITE가 가장 대표적으로 사용됩니다. 그러나 Targeted Host-associated Fungi ITS Database(THF)도 유명한 데이터베이스 중 하나입니다. 이를 dada2 버전으로 바꾼 다음에 사용해 봅시다.
1. 파일 다운로드
가장 최신 버전인 1.6.1을 다운로드하여보겠습니다.
- THFv1.6.1
- FASTA sequencesv1.6.1
위 파일은 QIIME2에 사용하기 좋은 형식으로 나누어져 있습니다. 즉 taxonomy annotation정보와 sequence정보가 따로 저장되어 있습니다.
이때 Accession 번호를 통해서, 각 서열에 매치된 데이터베이스의 정보를 파악할 수 있습니다.
그러나 DADA2 의 database의 형식은 아래와 같이, 서열의 이름에 계통 정보를 모두 포함합니다.
이는 ; 로 나누어져 있습니다.
즉 > 다음에 오는 accession 번호를 계통 정보로 바꾸어주면 손쉽게 사용 가능합니다.
이를 R에서 해볼까요?
2. Taxonomy 형식 바꾸기
저는 엑셀에서 왼쪽 형식을 오른쪽으로 먼저 바꾸어 주고, txt 파일로 저장하였습니다. 생각해 보니 이 작업도 R로 하는 것이 편했을 것 같습니다.
3. R에서 DADA2 format으로 변환하기
library(stringr)
library(readr)
library(glue)
fasta_data <- readLines( "~/Reference/dada2/THF_fungi/1.6.1/THFv1.6.1.sequence.fasta")
tax_data <- readr::read_table(~/Reference/dada2/THF_fungi/1.6.1/THFv1.6.1_for_DADA2.txt")
tax_data <- tax_data %>% data.frame()
간단하게 > 가 있는 줄에 해당하는 번호를 tax_data에서 찾아서 바꾸는 코드를 사용하였습니다.
for (i in 1:length(fasta_data)) {
if (substr(fasta_data[i], 1, 1) == ">") {
# fasta 파일의 각 헤더에서 이름 추출
fasta_name <- strsplit(fasta_data[i], ">")[[1]][2]
# tax 파일에서 일치하는 이름 찾기
matching_row <- tax_data[tax_data$Accession == fasta_name,"Taxonomy"]
# 대체할 이름이 있는 경우 이름 대체
if (!is.na(matching_row)) {
new_name <- matching_row
fasta_data[i] <- glue(">", new_name)
}
}
}
fasta_data %>% head
writeLines(fasta_data, "~/Reference/dada2/THF_fungi/1.6.1/THFv1.6.1.sequence_modi.fasta")
출처
반응형
