만약 샘플이 다르다면, 형성된 ASV를 기준으로 합치는 것은 매우 어렵습니다. 그러나 해상도가 조금 떨어지더라고 Species 기준으로는 두 샘플을 합칠 수 있습니다.
예제데이터를 사용하여 두 개의 phyloseq을 임의로 두 개로 나눈 다음에 다시 합쳐보겠습니다!
먼저 Phyloseq을 합치기 위한 조건이 있습니다.
1. 동일한 Metadata 속성을 가질 것
- 즉 sample_data()로 나오는 데이터의 colnames가 동일해야 합니다.
2. 전체 데이터를 Species-level로 합친 다음에, tax_table의 rownames이 Species여야 합니다.
3. count가 아닌 relative abundance에서만 합치는 것이 그나마..
pivot_wider를 이용해서 Phyloseq 합치기
- 이 함수는 melt 된 데이터를 다시 otu_table로 만들어주는 기능이 있습니다.
- phyloseq object를 문서로 저장할 수 없을 때, psmelt()로 이를 하나의 문서로 녹여내서 저장 가능합니다. 이때 melt 된 otu table을 다시 복구해 주는 함수가 바로 tidyr::pivot_wider입니다.
- 위 방법은 계통수를 합쳐주지 않습니다.
1. 예제 데이터 불러오기
library(ggplot2)
library(phyloseq)
library(dplyr)
`%!in%` <- Negate(`%in%`)
ps <- readRDS("./moving-picture.rds")
# phyloseq-class experiment-level object
# otu_table() OTU Table: [ 770 taxa and 34 samples ]
# sample_data() Sample Data: [ 34 samples by 8 sample variables ]
# tax_table() Taxonomy Table: [ 770 taxa by 7 taxonomic ranks ]
# phy_tree() Phylogenetic Tree: [ 770 tips and 768 internal nodes ]
2. 두 개의 phyloseq으로 분할하기
- Gut과 Gut이 아닌 샘플로 분할해 봅시다.
- 만약 두 샘플이 시퀀싱 런도 다르고, QIIME분석도 다르게 했다면, OTU ID를 기준으로 합칠 수 없습니다. 그러므로 우리는 최선의 방법인 Species-level로 합치도록 합니다.
# Species-level로 샘플 합치기
ps.sp <- tax_glom(ps, "Species") # 351 taxa and 34 samples
# 두 개의 phyloseq으로 분할
ps.gut <- subset_samples(ps.sp, body.site %in% "gut")
ps.gut <- prune_taxa(taxa_sums(ps.gut) >0, ps.gut) # 80 taxa and 8 samples
ps.other <- subset_samples(ps.sp, body.site %!in% "gut")
ps.other <- prune_taxa(taxa_sums(ps.other) >0, ps.other) # 328 taxa and 26 samples
3. melt데이터 합치기
melt.gut <- psmelt(ps.gut)
melt.other <- psmelt(ps.other)
melt.total <- rbind(melt.gut, melt.other)
4. Tax table 합치기
melt 된 데이터에서 유니크한 taxa를 추출합니다. 거기에 Species_001등으로 순번을 부여해 봅니다.
# taxonomy 의 공백 변환하기 (데이터에 따라 다름)
melt.total$Species_2 <- gsub(pattern = " ", replacement = "_", x = melt.total$Species)
melt.total$Species_2 <- gsub(pattern = " ", replacement = "_", x = melt.total$Species_2)
# 각 taxonomic rank로 그룹화하여 unique한 taxa 식별
unique_taxa <- melt.total %>%
distinct(Kingdom, Phylum, Class, Order, Family, Genus, Species) %>%
mutate(Species_num = paste0("Species_", sprintf("%03d", row_number())))
dim(unique_taxa)
# 통합된 tax data 만들기
merge.tax <- unique_taxa
# rownames은 Species_001등으로 바꾸어 주기
rownames(merge.tax) <- merge.tax$Species_num
5. OTU table 합치기
# 부여된 Species이름을 다시 melt데이터에 추가
melt.total2 <- melt.total %>%
left_join(unique_taxa, by = c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species"))
merge.otu <- melt.total2 %>%
# 합칠 데이터셋 추출
select(Species_num, Sample, Abundance) %>%
# OTU table형식으로 복구하기
pivot_wider(names_from = Sample, values_from = Abundance2) %>%
select(Species_num, everything())
# NA값은 0 으로 변환
merge.otu[is.na(merge.otu)] <- 0
merge.otu <- merge.otu %>% tibble::column_to_rownames("Species_num")
head(merge.otu)
# L1S105 L1S140 L1S208 L1S257 L1S281
# Species_001 3075 3966 3585 2075 2500
# Species_002 1197 647 964 909 531
# Species_003 910 711 509 250 565
# Species_004 196 172 626 738 707
# Species_005 664 13 285 301 648
6. Metadata 합치기
merge.meta <- melt.total %>%
select(Sample, body.site, subject, reported.antibiotic.usage) %>%
unique
rownames(merge.meta) <- merge.meta$Sample
7. Phyloseq 복구하기
phy <- phyloseq(otu_table(as.matrix(merge.otu), taxa_are_rows = T),
tax_table(as.matrix(merge.tax)),
sample_data(merge.meta)
)
phy
# phyloseq-class experiment-level object
# otu_table() OTU Table: [ 351 taxa and 34 samples ]
# sample_data() Sample Data: [ 34 samples by 4 sample variables ]
# tax_table() Taxonomy Table: [ 351 taxa by 8 taxonomic ranks ]
자 두 개의 phyloseq을 합쳤을 때 샘플은 총 34개, taxa는 총 351개입니다.
이는 기존 phyloseq을 Species-level로 합친 결과와 동일합니다.
ps.sp <- tax_glom(ps, "Species")
ps.sp
# phyloseq-class experiment-level object
# otu_table() OTU Table: [ 351 taxa and 34 samples ]
# sample_data() Sample Data: [ 34 samples by 8 sample variables ]
# tax_table() Taxonomy Table: [ 351 taxa by 7 taxonomic ranks ]
# phy_tree() Phylogenetic Tree: [ 351 tips and 350 internal nodes ]
만약 두 개의 phyloseq이 서열데이터도 가지고 있다면, R 상에서 tree를 제작 후 삽입할 수 있을 것입니다.
자 이렇게 수동으로 데이터를 선택하여 하나의 object로 만들었습니다.
저는 opendataset의 16S rRNA의 증폭영역이 다른 데이터를 이러한 형식으로 합쳐서 분석하고 있습니다.
더 좋은 방법이 있다면 추천해주세요!
참고로 psmelt 된 데이터는 데이터가 커질수록 컴퓨터에서 매우 힘들어함으로 intersect()를 활용하는 것도 좋은 방법입니다.
감사합니다.
