方法簡述
- 將推定惡性細(xì)胞比例調(diào)整至20%以下
- 過濾掉信息量較少的基因(默認(rèn):少于 10 個(gè)細(xì)胞表達(dá)毁葱、在 log2 尺度上的平均表達(dá)少于 0.1)
- 轉(zhuǎn)化為Z-score,將尺度限制為-3~3
- 按染色體位置對基因進(jìn)行排序椭迎,估計(jì) CNV 信號(默認(rèn) = 100 個(gè)基因)
- 用兩個(gè)參數(shù)匯總CNV信號庶溶,對惡性細(xì)胞和非惡性細(xì)胞進(jìn)行分類而咆。
- CNV 信號(MS买决,平方均值):CNV信號估計(jì)均方惯疙。
- CORR(與高 CNV 信號細(xì)胞的相關(guān)性):每個(gè)細(xì)胞的 CNV 與前 5% 細(xì)胞的平均值的相關(guān)性呕乎。
- 如果它們的 CNV 信號 (MS) > 0.02(X軸) 或 CNV 相關(guān)性 (CORR) > 0.2(Y軸)示罗,則鑒定為惡性細(xì)胞惩猫。
代碼可通過原作者Github獲取。這里以作者給出的example進(jìn)行示例蚜点。
library(ggplot2)
library(pals)
library(plyr)
library(dplyr)
library(Seurat)
library(gplots)
library(RColorBrewer)
####################################################################################################
source("R/calInferredCNA_for_CEP.R") ## calculate inferredCNV value
source("R/makingTCIDEA_for_CEP.R") ## make TCIDEA object
source("R/TCIDEA_obj_for_CEP.R") ## initiate TCIDEA object
source("R/calCNVScore_for_CEP.R") ## calculate CNV score (MS, Corr) with CEP result
###################################################################################################
## Example data ###################################################################################
#腫瘤樣本注釋
cell_annotation_with_tumor <- readRDS(file = "example/cell_info_tumor_example.Rds") # celltype = annotation cell types in transcriptome data
#腫瘤樣本log2表達(dá)矩陣
tumor_example <- readRDS(file = "example/log2TPM_tumor_example.Rds")
#正常細(xì)胞log2表達(dá)矩陣
normal_example <- readRDS(file = "example/log2TPM_normal_example.Rds")
#基因組位置信息
ref_genome_example <- readRDS(file = "example/refgenome_example.Rds")
#輸出目錄
output.dir = paste0(getwd(), "/", "result")
###################################################################################################
## PARAMETERS #####################################################################################
#設(shè)定推斷惡性細(xì)胞最大比例
EP_cutoff = 20 ## count if > 20% of EP -> add
#ep細(xì)胞亞群的名稱
target.celltypes = "EP" ## declare name of epithelial cells in metadata
label = "example"
###################################################################################################
## 1. check proportion of epithelial cells in tumor tissues ##
prop <- as.data.frame(table(cell_annotation_with_tumor$celltype))
prop$Percent = prop$Freq / nrow(cell_annotation_with_tumor) * 100
##2. Check the proportion (adding normal cells or not)
if(prop[prop$Var1 %in% target.celltypes,]$Percent > EP_cutoff){
list <- addNormalDataset(tumor.data = tumor_example, tumor.ident = cell_annotation_with_tumor, target.celltypes = target.celltypes,
normal.data = normal_example)
addnormal_example <- as.matrix(list$data); addnormal_annotation <- list$ident
runCEP(target.normalized = addnormal_example,
sample.info = addnormal_annotation, label = paste0(label,"_AddNormal"),
annotationdata = ref_genome_example, target.celltypes = target.celltypes, output.dir = output.dir,
min.cells = 10, MYwalk = 100) ## Sample list of EP proportion > EP_cutoff (20%)
}else{
runCEP(target.normalized = tumor_example,
sample.info = cell_annotation_with_tumor, label = label,
annotationdata = ref_genome_example,target.celltypes = target.celltypes, output.dir = output.dir,
min.cells = 10, MYwalk = 100) ## Sample list of EP proportion <= EP_cutoff (20%)
}
############################################
結(jié)果提取
result=readRDS(file = "result/example_AddNormal_after_calc_CNV_score.Rds")
identification=data.frame(Barcord=result$Index,cell_type=result$cell_index)
head(identification)
后續(xù)自行整合至seurat注釋進(jìn)行后續(xù)分析
參考來源:
https://github.com/SGI-LungCancer/SingleCell
參考文獻(xiàn):
Single-cell RNA sequencing demonstrates the molecular and cellular reprogramming of metastatic lung adenocarcinoma
鳴謝:
I thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.
問題交流:
Email: xuran@hrbmu.edu.cn