单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析5
代码 分析 解析 单细胞 转录 癌症 可及 染色质
2023-06-13 09:11:31 时间
单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析1:https://cloud.tencent.com/developer/article/2055573
单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析2:https://cloud.tencent.com/developer/article/2072069
单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析3:https://cloud.tencent.com/developer/article/2078159
单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析4:https://cloud.tencent.com/developer/article/2078348
代码解析
##这里是接的上一个教程的内容,由于代码过长,分成了两个部分
# Identify immune clusters
#######################################################
# Find immune cells by relative enrichment of ESTIMATE immune signature
##
library(psych)
test <- VlnPlot(rna,features = "immune.2")
data <- describeBy(test$data$immune.2, test$data$ident, mat = TRUE)
data.immune <- dplyr::filter(data,median > 0.1)
test <- VlnPlot(rna,features = "plasma.7")
data <- describeBy(test$data$plasma.7, test$data$ident, mat = TRUE)
data.plasma <- dplyr::filter(data,median > 0.1)
immune.clusters <- intersect(data.immune$group1,levels(Idents(rna)))
plasma.clusters <- intersect(data.plasma$group1,levels(Idents(rna)))
immune.clusters <- unique(append(immune.clusters,plasma.clusters))
for (i in 1:length(immune.clusters)){
j <- which(levels(Idents(rna)) == immune.clusters[i])
levels(Idents(rna))[j] <- paste0("immune.",immune.clusters[i])
}
rna@meta.data$postdoublet.idents <- Idents(rna)
idents <- data.frame(rownames(rna@meta.data),rna@meta.data$postdoublet.idents)
saveRDS(rna,"./rna_postdoublet_preinferCNV.rds")
#MAke inferCNV input
colnames(idents) <- c("V1","V2")
rownames(idents) <- NULL
colnames(idents) <- NULL
write.table(idents,"./sample_annotation_file_inferCNV.txt",sep = "\t",row.names = FALSE)
idents <- read.delim("./sample_annotation_file_inferCNV.txt",header = F)
gtf <- read.delim(GRCH38.annotations,header = F)
library(EnsDb.Hsapiens.v86)
convert.symbol = function(Data){
ensembls <- Data$V1
ensembls <- gsub("\\.[0-9]*$", "", ensembls)
geneIDs1 <- ensembldb::select(EnsDb.Hsapiens.v86, keys= ensembls, keytype = "GENEID", columns = "SYMBOL")
Data <- cbind.fill(Data, geneIDs1, fill = NA)
Data <- na.omit(Data)
Data$feature <- Data$SYMBOL
Data.new <- data.frame(Data$SYMBOL,Data$V2,Data$V3,Data$V4)
Data.new$Data.V2 <- paste("chr",Data.new$Data.V2,sep = "")
Data.new$Data.SYMBOL <- make.unique(Data.new$Data.SYMBOL)
return(Data.new)
}
gtf <- convert.symbol(gtf)
head(gtf)
write.table(gtf,"./Homo_sapiens.GRCh38.86.symbol.txt",sep = "\t",row.names = FALSE,col.names = FALSE)
num.immune.clusters = length(immune.clusters)
# create the infercnv object
if ( num.immune.clusters == 1) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=paste0("immune.",immune.clusters[1]))
} else if (num.immune.clusters == 2) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2])))
} else if ( num.immune.clusters == 3 ) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3])))
} else if (num.immune.clusters == 4) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3]),
paste0("immune.",immune.clusters[4])))
} else if (num.immune.clusters == 5) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3]),
paste0("immune.",immune.clusters[4]),
paste0("immune.",immune.clusters[5])))
} else if (num.immune.clusters == 6) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3]),
paste0("immune.",immune.clusters[4]),
paste0("immune.",immune.clusters[5]),
paste0("immune.",immune.clusters[6])))
}else if (num.immune.clusters == 7) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3]),
paste0("immune.",immune.clusters[4]),
paste0("immune.",immune.clusters[5]),
paste0("immune.",immune.clusters[6]),
paste0("immune.",immune.clusters[7])))
}else if (num.immune.clusters == 8) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3]),
paste0("immune.",immune.clusters[4]),
paste0("immune.",immune.clusters[5]),
paste0("immune.",immune.clusters[6]),
paste0("immune.",immune.clusters[7]),
paste0("immune.",immune.clusters[8])))
}else if (num.immune.clusters == 9) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3]),
paste0("immune.",immune.clusters[4]),
paste0("immune.",immune.clusters[5]),
paste0("immune.",immune.clusters[6]),
paste0("immune.",immune.clusters[7]),
paste0("immune.",immune.clusters[8]),
paste0("immune.",immune.clusters[9])))
}else if (num.immune.clusters == 10) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3]),
paste0("immune.",immune.clusters[4]),
paste0("immune.",immune.clusters[5]),
paste0("immune.",immune.clusters[6]),
paste0("immune.",immune.clusters[7]),
paste0("immune.",immune.clusters[8]),
paste0("immune.",immune.clusters[9]),
paste0("immune.",immune.clusters[10])))
}
# perform infercnv operations to reveal cnv signal
infercnv_obj = infercnv::run(infercnv_obj,
cutoff=0.1, # use 1 for smart-seq, 0.1 for 10x-genomics
out_dir="./output_dir_CNV_postdoublet_FailedCNVTest", # dir is auto-created for storing outputs
cluster_by_groups=T, # cluster
denoise=T,scale_data = T,
HMM=T,HMM_type = "i6",analysis_mode = "samples",min_cells_per_gene = 10,
BayesMaxPNormal = 0.4, num_threads = 8
)
regions <- read.delim("./output_dir_CNV_postdoublet_FailedCNVTest/HMM_CNV_predictions.HMMi6.hmm_mode-samples.Pnorm_0.4.pred_cnv_regions.dat")
probs <- read.delim("./output_dir_CNV_postdoublet_FailedCNVTest/BayesNetOutput.HMMi6.hmm_mode-samples/CNV_State_Probabilities.dat")
probs <- as.data.frame(t(probs[3,]))
colnames(probs) <- "Prob.Normal"
probs <- dplyr::filter(probs,Prob.Normal < 0.05)
cnvs <- rownames(probs)
cnvs <- gsub("\\.","-",cnvs)
regions <- regions[regions$cnv_name %in% cnvs, ]
cnv.groups <- sub("\\..*", "", regions$cell_group_name)
length(which(rownames(rna@reductions$pca@cell.embeddings) == rownames(rna@meta.data)))
rna$PC1.loading <- rna@reductions$pca@cell.embeddings[,1]
rna$cell.barcode <- rownames(rna@meta.data)
rna$CNV.Pos <- ifelse(as.character(rna$postdoublet.idents) %in% cnv.groups,1,0)
cnv.freq <- data.frame(table(regions$cell_group_name))
cnv.freq$Var1 <- sub("\\..*", "", cnv.freq$Var1)
rna$Total_CNVs <- ifelse(as.character(rna$postdoublet.idents) %in% cnv.freq$Var1,cnv.freq$Freq,0)
saveRDS(rna,"./rna_postdoublet_FailedCNVTest.rds")
}
}else{
all.genes <- rownames(rna)
rna <- ScaleData(rna, features = all.genes,vars.to.regress = "nCount_RNA")
rna <- FindNeighbors(rna,dims = 1:50)
rna <- FindClusters(rna,resolution = 0.7)
rna <- RunUMAP(rna,dims = 1:50)
Idents(rna) <- "RNA_snn_res.0.7"
# Proceed with CNV
library(infercnv)
library(stringr)
library(Seurat)
counts_matrix = GetAssayData(rna, slot="counts")
# Identify immune clusters
#######################################################
# Find immune cells by relative enrichment of ESTIMATE immune signature
library(psych)
test <- VlnPlot(rna,features = "immune.2")
data <- describeBy(test$data$immune.2, test$data$ident, mat = TRUE)
data.immune <- dplyr::filter(data,median > 0.1)
test <- VlnPlot(rna,features = "plasma.7")
data <- describeBy(test$data$plasma.7, test$data$ident, mat = TRUE)
data.plasma <- dplyr::filter(data,median > 0.1)
immune.clusters <- intersect(data.immune$group1,levels(Idents(rna)))
plasma.clusters <- intersect(data.plasma$group1,levels(Idents(rna)))
immune.clusters <- unique(append(immune.clusters,plasma.clusters))
for (i in 1:length(immune.clusters)){
j <- which(levels(Idents(rna)) == immune.clusters[i])
levels(Idents(rna))[j] <- paste0("immune.",immune.clusters[i])
}
rna@meta.data$postdoublet.idents <- Idents(rna)
idents <- data.frame(rownames(rna@meta.data),rna@meta.data$postdoublet.idents)
colnames(idents) <- c("V1","V2")
saveRDS(rna,"./rna_postdoublet_preinferCNV.rds")
# Make inferCNV inputs
rownames(idents) <- NULL
colnames(idents) <- NULL
write.table(idents,"./sample_annotation_file_inferCNV.txt",sep = "\t",row.names = FALSE)
idents <- read.delim("./sample_annotation_file_inferCNV.txt",header = F)
gtf <- read.delim(GRCH38.annotations,header = F)
library(EnsDb.Hsapiens.v86)
convert.symbol = function(Data){
ensembls <- Data$V1
ensembls <- gsub("\\.[0-9]*$", "", ensembls)
geneIDs1 <- ensembldb::select(EnsDb.Hsapiens.v86, keys= ensembls, keytype = "GENEID", columns = "SYMBOL")
Data <- cbind.fill(Data, geneIDs1, fill = NA)
Data <- na.omit(Data)
Data$feature <- Data$SYMBOL
Data.new <- data.frame(Data$SYMBOL,Data$V2,Data$V3,Data$V4)
Data.new$Data.V2 <- paste("chr",Data.new$Data.V2,sep = "")
Data.new$Data.SYMBOL <- make.unique(Data.new$Data.SYMBOL)
return(Data.new)
}
gtf <- convert.symbol(gtf)
head(gtf)
write.table(gtf,"./Homo_sapiens.GRCh38.86.symbol.txt",sep = "\t",row.names = FALSE,col.names = FALSE)
num.immune.clusters = length(immune.clusters)
# create the infercnv object
if ( num.immune.clusters == 1) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=paste0("immune.",immune.clusters[1]))
} else if (num.immune.clusters == 2) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2])))
} else if ( num.immune.clusters == 3 ) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3])))
} else if (num.immune.clusters == 4) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3]),
paste0("immune.",immune.clusters[4])))
} else if (num.immune.clusters == 5) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3]),
paste0("immune.",immune.clusters[4]),
paste0("immune.",immune.clusters[5])))
} else if (num.immune.clusters == 6) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3]),
paste0("immune.",immune.clusters[4]),
paste0("immune.",immune.clusters[5]),
paste0("immune.",immune.clusters[6])))
}else if (num.immune.clusters == 7) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3]),
paste0("immune.",immune.clusters[4]),
paste0("immune.",immune.clusters[5]),
paste0("immune.",immune.clusters[6]),
paste0("immune.",immune.clusters[7])))
}else if (num.immune.clusters == 8) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3]),
paste0("immune.",immune.clusters[4]),
paste0("immune.",immune.clusters[5]),
paste0("immune.",immune.clusters[6]),
paste0("immune.",immune.clusters[7]),
paste0("immune.",immune.clusters[8])))
}else if (num.immune.clusters == 9) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3]),
paste0("immune.",immune.clusters[4]),
paste0("immune.",immune.clusters[5]),
paste0("immune.",immune.clusters[6]),
paste0("immune.",immune.clusters[7]),
paste0("immune.",immune.clusters[8]),
paste0("immune.",immune.clusters[9])))
}else if (num.immune.clusters == 10) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3]),
paste0("immune.",immune.clusters[4]),
paste0("immune.",immune.clusters[5]),
paste0("immune.",immune.clusters[6]),
paste0("immune.",immune.clusters[7]),
paste0("immune.",immune.clusters[8]),
paste0("immune.",immune.clusters[9]),
paste0("immune.",immune.clusters[10])))
}
# perform infercnv operations to reveal cnv signal
infercnv_obj = infercnv::run(infercnv_obj,
cutoff=0.1, # use 1 for smart-seq, 0.1 for 10x-genomics
out_dir="./output_dir_CNV_postdoublet_FailedCorTest", # dir is auto-created for storing outputs
cluster_by_groups=T, # cluster
denoise=T,scale_data = T,
HMM=T,HMM_type = "i6",analysis_mode = "samples",min_cells_per_gene = 10,
BayesMaxPNormal = 0.4, num_threads = 8
)
regions <- read.delim("./output_dir_CNV_postdoublet_FailedCorTest/HMM_CNV_predictions.HMMi6.hmm_mode-samples.Pnorm_0.4.pred_cnv_regions.dat")
probs <- read.delim("./output_dir_CNV_postdoublet_FailedCorTest/BayesNetOutput.HMMi6.hmm_mode-samples/CNV_State_Probabilities.dat")
probs <- as.data.frame(t(probs[3,]))
colnames(probs) <- "Prob.Normal"
probs <- dplyr::filter(probs,Prob.Normal < 0.05)
cnvs <- rownames(probs)
cnvs <- gsub("\\.","-",cnvs)
regions <- regions[regions$cnv_name %in% cnvs, ]
cnv.groups <- sub("\\..*", "", regions$cell_group_name)
length(which(rownames(rna@reductions$pca@cell.embeddings) == rownames(rna@meta.data)))
rna$PC1.loading <- rna@reductions$pca@cell.embeddings[,1]
rna$cell.barcode <- rownames(rna@meta.data)
rna$CNV.Pos <- ifelse(as.character(rna$postdoublet.idents) %in% cnv.groups,1,0)
cnv.freq <- data.frame(table(regions$cell_group_name))
cnv.freq$Var1 <- sub("\\..*", "", cnv.freq$Var1)
rna$Total_CNVs <- ifelse(as.character(rna$postdoublet.idents) %in% cnv.freq$Var1,cnv.freq$Freq,0)
saveRDS(rna,"./rna_postdoublet_FailedCorTest.rds")
}
}else{
rna <- FindNeighbors(rna,dims = 1:50)
rna <- FindClusters(rna,resolution = 0.7)
rna <- RunUMAP(rna,dims = 1:50)
Idents(rna) <- "RNA_snn_res.0.7"
# Proceed with CNV
library(infercnv)
library(stringr)
library(Seurat)
counts_matrix = GetAssayData(rna, slot="counts")
# Identify immune clusters
#######################################################
# Find immune cells by relative enrichment of ESTIMATE immune signature
library(psych)
test <- VlnPlot(rna,features = "immune.2")
data <- describeBy(test$data$immune.2, test$data$ident, mat = TRUE)
data.immune <- dplyr::filter(data,median > 0.1)
test <- VlnPlot(rna,features = "plasma.7")
data <- describeBy(test$data$plasma.7, test$data$ident, mat = TRUE)
data.plasma <- dplyr::filter(data,median > 0.1)
immune.clusters <- intersect(data.immune$group1,levels(Idents(rna)))
plasma.clusters <- intersect(data.plasma$group1,levels(Idents(rna)))
immune.clusters <- unique(append(immune.clusters,plasma.clusters))
for (i in 1:length(immune.clusters)){
j <- which(levels(Idents(rna)) == immune.clusters[i])
levels(Idents(rna))[j] <- paste0("immune.",immune.clusters[i])
}
rna@meta.data$postdoublet.idents <- Idents(rna)
idents <- data.frame(rownames(rna@meta.data),rna@meta.data$postdoublet.idents)
colnames(idents) <- c("V1","V2")
saveRDS(rna,"./rna_postdoublet_preinferCNV.rds")
# Make inferCNV inputs
rownames(idents) <- NULL
colnames(idents) <- NULL
write.table(idents,"./sample_annotation_file_inferCNV.txt",sep = "\t",row.names = FALSE)
idents <- read.delim("./sample_annotation_file_inferCNV.txt",header = F)
gtf <- read.delim(GRCH38.annotations,header = F)
library(EnsDb.Hsapiens.v86)
convert.symbol = function(Data){
ensembls <- Data$V1
ensembls <- gsub("\\.[0-9]*$", "", ensembls)
geneIDs1 <- ensembldb::select(EnsDb.Hsapiens.v86, keys= ensembls, keytype = "GENEID", columns = "SYMBOL")
Data <- cbind.fill(Data, geneIDs1, fill = NA)
Data <- na.omit(Data)
Data$feature <- Data$SYMBOL
Data.new <- data.frame(Data$SYMBOL,Data$V2,Data$V3,Data$V4)
Data.new$Data.V2 <- paste("chr",Data.new$Data.V2,sep = "")
Data.new$Data.SYMBOL <- make.unique(Data.new$Data.SYMBOL)
return(Data.new)
}
gtf <- convert.symbol(gtf)
head(gtf)
write.table(gtf,"./Homo_sapiens.GRCh38.86.symbol.txt",sep = "\t",row.names = FALSE,col.names = FALSE)
num.immune.clusters = length(immune.clusters)
# create the infercnv object
# create the infercnv object
if ( num.immune.clusters == 1) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=paste0("immune.",immune.clusters[1]))
} else if (num.immune.clusters == 2) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2])))
} else if ( num.immune.clusters == 3 ) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3])))
} else if (num.immune.clusters == 4) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3]),
paste0("immune.",immune.clusters[4])))
} else if (num.immune.clusters == 5) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3]),
paste0("immune.",immune.clusters[4]),
paste0("immune.",immune.clusters[5])))
} else if (num.immune.clusters == 6) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3]),
paste0("immune.",immune.clusters[4]),
paste0("immune.",immune.clusters[5]),
paste0("immune.",immune.clusters[6])))
}else if (num.immune.clusters == 7) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3]),
paste0("immune.",immune.clusters[4]),
paste0("immune.",immune.clusters[5]),
paste0("immune.",immune.clusters[6]),
paste0("immune.",immune.clusters[7])))
}else if (num.immune.clusters == 8) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3]),
paste0("immune.",immune.clusters[4]),
paste0("immune.",immune.clusters[5]),
paste0("immune.",immune.clusters[6]),
paste0("immune.",immune.clusters[7]),
paste0("immune.",immune.clusters[8])))
}else if (num.immune.clusters == 9) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3]),
paste0("immune.",immune.clusters[4]),
paste0("immune.",immune.clusters[5]),
paste0("immune.",immune.clusters[6]),
paste0("immune.",immune.clusters[7]),
paste0("immune.",immune.clusters[8]),
paste0("immune.",immune.clusters[9])))
}else if (num.immune.clusters == 10) {
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=as.matrix(counts_matrix),
annotations_file="./sample_annotation_file_inferCNV.txt",
gene_order_file="./Homo_sapiens.GRCh38.86.symbol.txt",
ref_group_names=c(paste0("immune.",immune.clusters[1]),
paste0("immune.",immune.clusters[2]),
paste0("immune.",immune.clusters[3]),
paste0("immune.",immune.clusters[4]),
paste0("immune.",immune.clusters[5]),
paste0("immune.",immune.clusters[6]),
paste0("immune.",immune.clusters[7]),
paste0("immune.",immune.clusters[8]),
paste0("immune.",immune.clusters[9]),
paste0("immune.",immune.clusters[10])))
}
# perform infercnv operations to reveal cnv signal
infercnv_obj = infercnv::run(infercnv_obj,
cutoff=0.1, # use 1 for smart-seq, 0.1 for 10x-genomics
out_dir="./output_dir_CNV_postdoublet_SkipChecks", # dir is auto-created for storing outputs
cluster_by_groups=T, # cluster
denoise=T,scale_data = T,
HMM=T,HMM_type = "i6",analysis_mode = "samples",min_cells_per_gene = 10,
BayesMaxPNormal = 0.4, num_threads = 8
)
regions <- read.delim("./output_dir_CNV_postdoublet_SkipChecks/HMM_CNV_predictions.HMMi6.hmm_mode-samples.Pnorm_0.4.pred_cnv_regions.dat")
probs <- read.delim("./output_dir_CNV_postdoublet_SkipChecks/BayesNetOutput.HMMi6.hmm_mode-samples/CNV_State_Probabilities.dat")
probs <- as.data.frame(t(probs[3,]))
colnames(probs) <- "Prob.Normal"
probs <- dplyr::filter(probs,Prob.Normal < 0.05)
cnvs <- rownames(probs)
cnvs <- gsub("\\.","-",cnvs)
regions <- regions[regions$cnv_name %in% cnvs, ]
cnv.groups <- sub("\\..*", "", regions$cell_group_name)
length(which(rownames(rna@reductions$pca@cell.embeddings) == rownames(rna@meta.data)))
rna$PC1.loading <- rna@reductions$pca@cell.embeddings[,1]
rna$cell.barcode <- rownames(rna@meta.data)
rna$CNV.Pos <- ifelse(as.character(rna$postdoublet.idents) %in% cnv.groups,1,0)
cnv.freq <- data.frame(table(regions$cell_group_name))
cnv.freq$Var1 <- sub("\\..*", "", cnv.freq$Var1)
rna$Total_CNVs <- ifelse(as.character(rna$postdoublet.idents) %in% cnv.freq$Var1,cnv.freq$Freq,0)
boxplot.cnv <- ggplot(rna@meta.data,aes(x= postdoublet.idents,y=PC1.loading,color = as.factor(CNV.Pos)))+geom_boxplot()
boxplot.cnv+ggsave("Postdoublet_CNV_PC1_boxplot.png")
saveRDS(rna,"./rna_postdoublet_SkipChecks.rds")
}
Idents(rna)<- "RNA_snn_res.0.7"
# DEG analysis with Wilcox
##########################################################################
# Wilcox
Wilcox.markers <- FindAllMarkers(object =rna, min.pct = 0.25,only.pos = F,
test.use = "wilcox")
saveRDS(Wilcox.markers,"./wilcox_DEGs.rds")
总结
这部分的内容其实与第二个教程后面的内容是一样的,代码都是相同的,所以只有去除了一个双细胞的内容。我比较好奇的是为社么要做两遍呢,已知是多细胞会对后面的分析有影响,直接做了就可以了,还又重复了一遍内容,难道是为了进行比对两个分析的结果吗。
相关文章
- 单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析3
- 单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析6
- 单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析16
- 主题挖掘LDA和情感分析图书馆话题知乎用户问答行为数据|附代码数据
- 杭州出租车行驶轨迹数据空间时间可视化分析|附代码数据
- 主题挖掘LDA和情感分析图书馆话题知乎用户问答行为数据|附代码数据
- 数据分享|R语言用主成分PCA、 逻辑回归、决策树、随机森林分析心脏病数据并高维可视化|附代码数据
- 【视频】R语言逻辑回归(Logistic回归)模型分类预测病人冠心病风险|数据分享|附代码数据
- R语言中贝叶斯网络(BN)、动态贝叶斯网络、线性模型分析错颌畸形数据|附代码数据
- Python配对交易策略统计套利量化交易分析股票市场|附代码数据
- 【Netty】Netty 入门案例分析 ( Netty 模型解析 | Netty 服务器端代码 | Netty 客户端代码 )
- 【Android 逆向】使用 Python 代码解析 ELF 文件 ( PyCharm 中进行断点调试 | ELFFile 实例对象分析 )
- [五]类加载机制双亲委派机制 底层代码实现原理 源码分析 java类加载双亲委派机制是如何实现的详解编程语言
- 关系分析利用Neo4j挖掘代码关系(neo4j代码)
- 案例Oracle事务回滚实践代码案例分析(oracle事务回滚代码)
- PHP用SAX解析XML的实现代码与问题分析
- 不能在子类或外部类发布C#事件代码分析
- 分析python服务器拒绝服务攻击代码
- JS中实现简单Formatter函数示例代码