-
We have this “2N” scDNA data from the blood cells of zebrafish. Here we have 44 arms (columns) in each cells (row). 1 means this arm is present, and 0 means that this arm is absent.
Binary chromosome-arm calls for zebrafish blood cells -
The sequencing results of this 2N data are supposed to be good (all arms in all cells should be present).
-
However, we still managed to find a way to do QC (Remove bad amplicons).
- Idea ⇒ If looking at arm_1, how many cells have both arm_1 and arm_2 to be 1? how many cells have both arm_1 and arm_3 to be 1?…etc. This calculates the POSSIBILITY/FREQUENCY of appearance of other arms… when arm_1 is present in the cell.
# i,j 是 arm 的編號
# i 是要拿出來放在colnames 的 arm
# j 是跟 i 做比較，當i arm = 1時，j arm = 1 的細胞有多少個
for (i in (1:44)) {
for (j in (1:44)) {
query_arm <- paste0("arm_",i,"$")
compare_arm <- paste0("arm_",j,"$")
assign("query_arm_df",
rename_norm_2N_3amp_HighThreshold[grepl(query_arm,
names(rename_norm_2N_3amp_HighThreshold),
perl = TRUE)])
assign("compare_arm_df",
rename_norm_2N_3amp_HighThreshold[grepl(compare_arm,
names(rename_norm_2N_3amp_HighThreshold),
perl = TRUE)])
df <- cbind(query_arm_df,
compare_arm_df)
df$SUM <- rowSums(df)
assign(paste0("arm_",i,"_compare_with_arm_",j),sum(df$SUM == 2),envir=.GlobalEnv)
}
}Tidy the above result to data frame：
for (i in (1:44)) {
# Create character vector with variable names
var_names <- paste0("arm_",i,"_compare_with_arm_",1:44)
# Retrieve values of variables as a list
var_values <- mget(var_names)
# Combine list elements into a single vector
lst <- do.call(c, var_values)
# Create data frame
assign(paste0("arm_",i),data.frame(lst),envir=.GlobalEnv)
}
# cbind
# Create character vector with variable names
var_df_names <- paste0("arm_",1:44)
# Retrieve values of variables as a list
var_df <- mget(var_df_names)
arm_corr_hm_df_2N <- do.call(cbind, var_df)
rownames(arm_corr_hm_df_2N) <-c("chr1p","chr1q",
"chr2p",
"chr3q",
"chr4p","chr4q",
"chr5q",
"chr6p","chr6q",
"chr7p","chr7q",
"chr8q",
"chr9p","chr9q",
"chr10p","chr10q",
"chr11p","chr11q",
"chr12p","chr12q",
"chr13p","chr13q",
"chr14p","chr14q",
"chr15p","chr15q",
"chr16p","chr16q",
"chr17p","chr17q",
"chr18p","chr18q",
"chr19p","chr19q",
"chr20q",
"chr21p","chr21q",
"chr22p","chr22q",
"chr23p","chr23q",
"chr24p","chr24q",
"chr25p")
colnames(arm_corr_hm_df_2N) <-c("chr1p","chr1q",
"chr2p",
"chr3q",
"chr4p","chr4q",
"chr5q",
"chr6p","chr6q",
"chr7p","chr7q",
"chr8q",
"chr9p","chr9q",
"chr10p","chr10q",
"chr11p","chr11q",
"chr12p","chr12q",
"chr13p","chr13q",
"chr14p","chr14q",
"chr15p","chr15q",
"chr16p","chr16q",
"chr17p","chr17q",
"chr18p","chr18q",
"chr19p","chr19q",
"chr20q",
"chr21p","chr21q",
"chr22p","chr22q",
"chr23p","chr23q",
"chr24p","chr24q",
"chr25p")
(These values are cell numbers.)
Normalization
for (i in (1:44)) {
assign(paste0("norm_arm_corr_",i),
(arm_corr_hm_df_2N[,i]/arm_corr_hm_df_2N[,i][i])*100)
# 整合成 dataframe
var_names <- paste0("norm_arm_corr_",i)
var_values <- mget(var_names)
lst <- do.call(c, var_values)
assign(paste0("norm_arm_corr_",i),data.frame(lst),envir=.GlobalEnv)
}
# cbind
var_df_names <- paste0("norm_arm_corr_",1:44)
var_df <- mget(var_df_names)
norm_arm_corr_hm_df_2N <- do.call(cbind, var_df)
rownames(norm_arm_corr_hm_df_2N) <- c("chr1p","chr1q",
"chr2p",
"chr3q",
"chr4p","chr4q",
"chr5q",
"chr6p","chr6q",
"chr7p","chr7q",
"chr8q",
"chr9p","chr9q",
"chr10p","chr10q",
"chr11p","chr11q",
"chr12p","chr12q",
"chr13p","chr13q",
"chr14p","chr14q",
"chr15p","chr15q",
"chr16p","chr16q",
"chr17p","chr17q",
"chr18p","chr18q",
"chr19p","chr19q",
"chr20q",
"chr21p","chr21q",
"chr22p","chr22q",
"chr23p","chr23q",
"chr24p","chr24q",
"chr25p")
colnames(norm_arm_corr_hm_df_2N) <- c("chr1p","chr1q",
"chr2p",
"chr3q",
"chr4p","chr4q",
"chr5q",
"chr6p","chr6q",
"chr7p","chr7q",
"chr8q",
"chr9p","chr9q",
"chr10p","chr10q",
"chr11p","chr11q",
"chr12p","chr12q",
"chr13p","chr13q",
"chr14p","chr14q",
"chr15p","chr15q",
"chr16p","chr16q",
"chr17p","chr17q",
"chr18p","chr18q",
"chr19p","chr19q",
"chr20q",
"chr21p","chr21q",
"chr22p","chr22q",
"chr23p","chr23q",
"chr24p","chr24q",
"chr25p")Since the sequencing qualities of each arm are different, normalization is necessary. For example, if arm_1 appear in only a few cells, the above values must be very small. But in fact, arm_1 could be what we want, because it is possible that everytime it appears in a cell, all other arms appear to be in the cell.
Visualization
my_palette <- colorRampPalette(c("white","orange","red"))(n = 100)
# clustering rows
hr_2N <- hclust(Dist(norm_arm_corr_hm_df_2N,
method = "euclidean"),
method = "ward.D2")
# heatmap上，rows 由下到上的順序分別是原本的第幾rows
hr_2N$order
# 把 column 的順序更改成與 cluster後的 row order 相同
norm_arm_corr_hm_df_2N <- norm_arm_corr_hm_df_2N[,c(rev(hr_2N$order))]
pdf(file = "norm_arm_corr_hm_percent_2N_3amp_HighThreshold.pdf")
par(cex.main = 0.8)
heatmap.2(as.matrix(norm_arm_corr_hm_df_2N),
Rowv=as.dendrogram(hr_2N),
Colv=F,
scale="none",
density.info="none",
col = my_palette,
trace="none",
margins = c(4,4),
cexRow = 0.5,
cexCol = 0.5,
lhei = c(0.2,1),
xlab = "Reference arm",
ylab = "Arms to be compared with",
key.xlab= "Frequency (%)",
main="Frquency of arm-arm appearing scenario of\n each reference arm,\n 2N (4481 cells), High Threshold (3 amplicons)")
dev.off()
Recorded QC criterion: at least 98% mean co-occurrence
How to explain chr19q?
- The chr19q amplicon falls below the recorded 98% mean co-occurrence threshold and is a candidate for exclusion in this QC analysis. → the amplicons that should be filtered out:
"chr1p" "chr2p" "chr17q" "chr19q" "chr24q"