Background
We have targeted 22 amplicons located on certain arm of certain chromosome in each cell. And on each amplicon, we have sequenced 1 to 5 SNPs which should be heterozygous. In total, there are 39 SNPs on these amplicons. The idea is that if this cell is a good cell, the ploidy of SNPs on the same amplicon should be consistent (i.e., if there are 3 SNPs on the amplicon “chr3q1”, all of them should be either missing, 1N or 2N). Overall, 10 amplicons which have more than 2 SNPs are used to filter out bad cells. If any of these 10 amplicons show the inconsistency, the cell is omitted.
Workflow
- Pre-filter sequencing depth and allele frequency for each SNP in each cell.
- For each amplicon, check the sequencing depth of SNPs on it: If the amplicon has 3 SNPs that are all present (sequencing depth ≠ 0) → the amplicon exist. If the amplicon has 3 SNPs that are all absent (sequencing depth = 0) → the amplicon is missing. Otherwise the amplicon is inconsistent → assign NA.
- Same procedure with allele frequency. The consistent amplicons should have their SNPs all 2N or all 1N.
- Following with the above two steps: If the amplicon exist → check what’s its allele frequency (ploidy). If the amplicon is missing or NA → do nothing.
- Eventually, all values of amplicons will be converted to 0 (missing), 1 (1N), 2 (2N) and NA (inconsistent).
- For each cell, we only look at the 10 amplicons which have more than 2 SNPs. Only keep the cells that all 10 amplicons are not NA.
Data Import
AF_62_mtx <- read_csv("62_SNPs_matrix/AF_matrix_62.csv")
depth_62_mtx <- read_csv("62_SNPs_matrix/depth_matrix_62.csv")
# good SNPs
SNPs_after_filter2 <- readRDS("SNPs_after_filter2_list.rds")
AF_39_mtx <- AF_62_mtx %>% select(sample, cell_barcode, all_of(SNPs_after_filter2))
depth_39_mtx <- depth_62_mtx %>% select(sample, cell_barcode, all_of(SNPs_after_filter2))
AF_39_mtx_2N <- AF_39_mtx[AF_39_mtx$sample == "Z_BC_2N_0706-3",]
AF_39_mtx_8d <- AF_39_mtx[AF_39_mtx$sample == "SEC-8dpf_1X_10X",]
AF_39_mtx_14d <- AF_39_mtx[AF_39_mtx$sample == "SEC-14dpf_1X_10X",]
depth_39_mtx_2N <- depth_39_mtx[depth_39_mtx$sample == "Z_BC_2N_0706-3",]
depth_39_mtx_8d <- depth_39_mtx[depth_39_mtx$sample == "SEC-8dpf_1X_10X",]
depth_39_mtx_14d <- depth_39_mtx[depth_39_mtx$sample == "SEC-14dpf_1X_10X",]
#######
AF_39_mtx_2N <- AF_39_mtx_2N[,-1] %>% tibble::column_to_rownames("cell_barcode")
AF_39_mtx_8d <- AF_39_mtx_8d[,-1] %>% tibble::column_to_rownames("cell_barcode")
AF_39_mtx_14d <- AF_39_mtx_14d[,-1] %>% tibble::column_to_rownames("cell_barcode")
depth_39_mtx_2N <- depth_39_mtx_2N[,-1] %>% tibble::column_to_rownames("cell_barcode")
depth_39_mtx_8d <- depth_39_mtx_8d[,-1] %>% tibble::column_to_rownames("cell_barcode")
depth_39_mtx_14d <- depth_39_mtx_14d[,-1] %>% tibble::column_to_rownames("cell_barcode")
pre-filter the AF and depth matrix
# AF
AF_39_mtx_2N <- ifelse(AF_39_mtx_2N > 99, 1, # if >99 -> 1N
ifelse(AF_39_mtx_2N < 1, 1, 2)) # if <1 -> 1N, otherwise 2N
AF_39_mtx_8d <- ifelse(AF_39_mtx_8d > 99, 1, # if >99 -> 1N
ifelse(AF_39_mtx_8d < 1, 1, 2)) # if <1 -> 1N, otherwise 2N
AF_39_mtx_14d <- ifelse(AF_39_mtx_14d > 99, 1, # if >99 -> 1N
ifelse(AF_39_mtx_14d < 1, 1, 2)) # if <1 -> 1N, otherwise 2N
# depth # if <10 -> 0 (absent), otherwise 1 (present)
depth_39_mtx_2N <- ifelse(depth_39_mtx_2N < 10, 0, 1)
depth_39_mtx_8d <- ifelse(depth_39_mtx_8d < 10, 0, 1)
depth_39_mtx_14d <- ifelse(depth_39_mtx_14d < 10, 0, 1)
AmpliconList_for_39SNPs <- c("chr3q1",
"chr6p4",
"chr7q1",
"chr8q4",
"chr11q1",
"chr12p2","chr12p4",
"chr13p3",
"chr14p2","chr14p3","chr14q1",
"chr15q4",
"chr16q3",
"chr17p2",
"chr18p1","chr18p2",
"chr19p1",
"chr21q2","chr21q4",
"chr22p1","chr22p2","chr25q3")Summarize the depth mtx for each amplicon
for (data in c("2N","8d","14d")) {
cat("Processing the dataset:", data, "\\n")
matrix <- get(paste0("depth_39_mtx_", data))
depth_df <- data.frame(row.names = rownames(matrix))
###
for (amplicon in AmpliconList_for_39SNPs) {
cat("Processing the amplicon:", amplicon,".....\\n")
df_of_each_amplicon <-
matrix %>%
as.data.frame() %>%
select(starts_with(amplicon))
summary_for_each_amplicon <-
df_of_each_amplicon %>%
mutate(!!amplicon := if_else(rowSums(.) == ncol(df_of_each_amplicon),
1, # if rowSums == SNPs number -> all SNPs present
if_else(rowSums(.) == 0, # if rowSums != SNPs number -> all SNPs absent
0, NA_real_))) %>% # otherwise -> inconsistent
#ungroup() %>%
select(!!sym(amplicon))
depth_df <- cbind(depth_df, summary_for_each_amplicon)
}
###
assign(paste0("depth_", data),
depth_df,
envir = .GlobalEnv)
rm(depth_df)
cat("==============================\\n")
}
Summarize the AF mtx for each amplicon
# For loop to summarize the ploidy in 3 datasets
for (data in c("2N","8d","14d")) {
cat("Processing the dataset:", data, "\\n")
matrix <- get(paste0("AF_39_mtx_", data))
ploidy_df <- data.frame(row.names = rownames(matrix))
###
for (amplicon in AmpliconList_for_39SNPs) {
cat("Processing the amplicon:", amplicon,".....\\n")
df_of_each_amplicon <-
matrix %>%
as.data.frame() %>%
select(starts_with(amplicon))
summary_for_each_amplicon <-
df_of_each_amplicon %>%
rowwise() %>%
mutate(!!amplicon := if_else(length(unique(c_across(starts_with(amplicon)))) == 1, # if they are consistent (should be)
as.numeric(c_across(starts_with(amplicon))[1]), # just make it 1(1N) or 2(2N)
NA_real_)) %>% # otherwise
ungroup() %>%
select(!!sym(amplicon))
ploidy_df <- cbind(ploidy_df, summary_for_each_amplicon)
}
###
assign(paste0("ploidy_", data),
ploidy_df,
envir = .GlobalEnv)
rm(ploidy_df)
cat("==============================\\n")
}
In a cell, if such amplicon present → Replaced it with the ploidy value in AF mtx (1, 2 or NA)
for (data in c("2N","8d","14d")) {
depth_mtx <- get(paste0("depth_", data)) %>% as.matrix()
ploidy_mtx <- get(paste0("ploidy_", data)) %>% as.matrix()
# Create a logical matrix indicating where depth_mtx has the value 1
indices <- (depth_mtx == 1) & !is.na(depth_mtx)
# Replace values in df1 with values from ploidy_mtx where indices is TRUE
depth_mtx[indices] <- ploidy_mtx[indices]
assign(paste0("final_mtx_", data), depth_mtx, envir = .GlobalEnv)
rm(depth_mtx)
rm(ploidy_mtx)
}Filter out bad cells
Filter out cells with any NA (inconsistency) among the selected ten multi-SNP amplicons
for (data in c("2N","8d","14d")) {
matrix <- get(paste0("final_mtx_", data))
final_mtx_filtered <- matrix %>%
as.data.frame() %>%
mutate(consistency = case_when(
(is.na(chr11q1) | is.na(chr12p2) | is.na(chr14p2) | is.na(chr14p3) | is.na(chr15q4) | is.na(chr16q3) | is.na(chr17p2) | is.na(chr18p1) | is.na(chr21q2) | is.na(chr22p1)) ~ "FALSE",
TRUE ~ "TRUE"))
final_mtx_filtered <- final_mtx_filtered[final_mtx_filtered$consistency == TRUE,]
final_mtx_filtered <- final_mtx_filtered[,-ncol(final_mtx_filtered)]
assign(paste0("final_mtx_filtered_", data),
final_mtx_filtered,
envir = .GlobalEnv)
rm(final_mtx_filtered)
rm(matrix)
}Visualization
8d
hr_8d <- hclust(Dist(final_mtx_filtered_8d, method = "euclidean"),
method = "ward.D2")
mycl_8d <- cutree(hr_8d, h = 30/100*max(hr_8d$height))
mycolhc <- rainbow(length(unique(mycl_8d)), start=1, end=0.5)
mycolhc <- mycolhc[as.vector(mycl_8d)]
myheatcol <- c("#000000","#0080FF","#FF9933")
png(filename = "PDF/Heatmap/8d.png", width = 450, height = 250, units = "mm", res = 800)
par(cex.main = 1.8)
heatmap.2(as.matrix(final_mtx_filtered_8d),
Rowv=as.dendrogram(hr_8d),
Colv = FALSE,
col=myheatcol,
scale="none",
density.info="none",
trace="none",
RowSideColors=mycolhc,
labRow = FALSE,
margins = c(8,4),
lhei = c(0.2,1),
cexCol = 1.5,
key.xlab= "ploidy (missing, 1N, 2N)",
main="8d data, ward.D2/euclidean,\\n 326 consistent cells; cutree: 30%")
dev.off()
