1. #理解inferCNV结果文件
    2. rm(list=ls())
    3. options(stringsAsFactors = F)
    4. library(phylogram)
    5. library(gridExtra)
    6. library(grid)
    7. require(dendextend)
    8. require(ggthemes)
    9. library(tidyverse)
    10. library(Seurat)
    11. library(infercnv)
    12. library(miscTools)
    14. #  Import inferCNV dendrogram
    15. infercnv.dend <- read.dendrogram(file = "inferCNV_output/infercnv.observations_dendrogram.txt")
    17. # Cut tree 
    18. infercnv.labels <- cutree(infercnv.dend,
    19.                           k = 6, 
    20.                           order_clusters_as_data = FALSE)
    21. table(infercnv.labels)
    22. # Color labels
    23. the_bars <- as.data.frame(tableau_color_pal("Tableau 20")(20)[infercnv.labels])
    24. colnames(the_bars) <- "inferCNV_tree"
    25. the_bars$inferCNV_tree <- as.character(the_bars$inferCNV_tree)
    27. infercnv.dend %>% set("labels",rep("", nobs(infercnv.dend)) )  %>% plot(main="inferCNV dendrogram") %>%
    28.   colored_bars(colors = as.data.frame(the_bars), dend = infercnv.dend, sort_by_labels_order = FALSE, add = T, y_scale=10, y_shift = 0)
    30. infercnv.labels=as.data.frame(infercnv.labels)
    31. groupFiles='groupFiles.txt'   
    32. meta=read.table(groupFiles,sep = '\t')
    33. infercnv.labels$V1=rownames(infercnv.labels)
    34. meta=merge(meta,infercnv.labels,by='V1')
    35. table(meta[,2:3]) 
    37. if( ! file.exists(  "cnv_scores.csv")){
    38.   tmp=read.table("inferCNV_output/infercnv.references.txt", header=T)
    39.   down=mean(rowMeans(tmp)) - 2 * mean( apply(tmp, 1, sd))
    40.   up=mean(rowMeans(tmp)) + 2 * mean( apply(tmp, 1, sd))
    41.   oneCopy=up-down
    42.   oneCopy
    43.   a1= down- 2*oneCopy
    44.   a2= down- 1*oneCopy
    45.   down;up
    46.   a3= up +  1*oneCopy
    47.   a4= up + 2*oneCopy 
    49.   cnv_table <- read.table("inferCNV_output/infercnv.observations.txt", header=T)
    50.   # Score cells based on their CNV scores 
    51.   # Replicate the table 
    52.   cnv_score_table <- as.matrix(cnv_table)
    54.   cnv_score_mat <- as.matrix(cnv_table)
    55.   # Scoring
    56.   cnv_score_table[cnv_score_mat > 0 & cnv_score_mat < a2] <- "A" #complete loss. 2pts
    57.   cnv_score_table[cnv_score_mat >= a2 & cnv_score_mat < down] <- "B" #loss of one copy. 1pts
    58.   cnv_score_table[cnv_score_mat >= down & cnv_score_mat <  up ] <- "C" #Neutral. 0pts
    59.   cnv_score_table[cnv_score_mat >= up  & cnv_score_mat <= a3] <- "D" #addition of one copy. 1pts
    60.   cnv_score_table[cnv_score_mat > a3  & cnv_score_mat <= a4 ] <- "E" #addition of two copies. 2pts
    61.   cnv_score_table[cnv_score_mat > a4] <- "F" #addition of more than two copies. 2pts
    63.   # Check
    64.   table(cnv_score_table[,1])
    65.   # Replace with score 
    66.   cnv_score_table_pts <- cnv_table
    67.   rm(cnv_score_mat)
    68.   # 
    69.   cnv_score_table_pts[cnv_score_table == "A"] <- 2
    70.   cnv_score_table_pts[cnv_score_table == "B"] <- 1
    71.   cnv_score_table_pts[cnv_score_table == "C"] <- 0
    72.   cnv_score_table_pts[cnv_score_table == "D"] <- 1
    73.   cnv_score_table_pts[cnv_score_table == "E"] <- 2
    74.   cnv_score_table_pts[cnv_score_table == "F"] <- 2
    76.   # Scores are stored in “cnv_score_table_pts”. Use colSums to add up scores for each cell and store as vector 
    77.   cell_scores_CNV <- as.data.frame(colSums(cnv_score_table_pts))
    78.   colnames(cell_scores_CNV) <- "cnv_score"
    79.   head(cell_scores_CNV)
    80.   write.csv(x = cell_scores_CNV, file = "cnv_scores.csv")
    82. } 
    85. # 除去了reference后的走inferCNV的细胞
    86. cell_scores_CNV=read.csv('cnv_scores.csv',row.names = 1)
    87. head(cell_scores_CNV) 
    90. load(file = '../section-01-cluster/basic.sce.pbmc.Rdata') 
    91. sce=pbmc 
    92. phe=sce@meta.data
    93. phe$celltype=Idents(sce)
    94. head(rownames(phe))
    95. head(rownames(cell_scores_CNV)) 
    97. #rownames(phe)=paste0('X',rownames(phe))
    98. rownames(phe)=gsub('-','.',rownames(phe))
    100. head(rownames(phe)) 
    101. head(rownames(cell_scores_CNV))
    103. head(rownames(phe))
    104. phe=phe[rownames(phe) %in% rownames(cell_scores_CNV),]
    105. identical(rownames(phe),rownames(cell_scores_CNV))
    107. infercnv.labels <- cutree(infercnv.dend, k = 6, order_clusters_as_data = FALSE)
    108. phe$inferCNV= infercnv.labels[match(rownames(phe), names(infercnv.labels) )]
    110. phe$cnv_scores  =  cell_scores_CNV[rownames(phe),]
    112. table(phe$celltype,phe$inferCNV) 
    113. head(rownames(phe))
    114. dim(phe)
    115. library(ggpubr)
    116. p1=ggboxplot(phe,'celltype','cnv_scores', fill = "celltype") 
    117. p1= p1+ theme(axis.text.x = element_text(angle = 45, 
    118.                                     vjust = 0.5, hjust=0.5))
    119. p2=ggboxplot(phe,'inferCNV','cnv_scores', fill = "inferCNV")
    120. library(patchwork)
    121. p1+p2
    122. ggsave(filename = 'anno_CNVscore.pdf')
    123. table(phe$celltype,phe$inferCNV)
    126. #rownames(phe)=gsub('X','',rownames(phe))
    127. rownames(phe)=gsub('[.]','-',rownames(phe))
    128. head(rownames(phe))
    130. sce
    131. tail(rownames(sce@meta.data))
    132. head(rownames(phe))
    133. sce$celltype=Idents(sce)
    134. table(sce$celltype)
    135. sce=subset(sce,celltype %in% c('epithelial' )) 
    136. sce
    137. kp=rownames(sce@meta.data) %in% rownames(phe)
    138. table(kp)
    139. sce=sce[,kp]
    141. phe=phe[rownames(sce@meta.data),]
    143. sce@meta.data=phe
    144. head(phe)
    145. save(sce, file = 'epi_sce_annoCNV.Rdata')