第一步 获取对应物种的基因位置注释文件
- https://www.gencodegenes.org/human/
- wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_39/gencode.v39.chr_patch_hapl_scaff.annotation.gtf.gz
- gunzip gencode.v39.chr_patch_hapl_scaff.annotation.gtf.gz
- hg39 <- read.table('gencode.v39.chr_patch_hapl_scaff.annotation.gtf', sep = '\t', comment.char = "#")
- hg39['gene_id'] <- stringr::str_extract(hg39$V9,"(?<=gene_id ).+?(?=;)")
- hg39['gene_name'] <- stringr::str_extract(hg39$V9,"(?<=gene_name ).+?(?=;)")
- 或者直接使用TCGA的基因位置注释：https://file-cdn.limour.top/bix/22.02.18.tcga_rowRanges.csv.gz
第二步 获取分析数据的基因位置信息
f_df_order <- function(df, oN, groupN, decreasing = T){
df <- as.data.frame(df)
for(oi in levels(df[[groupN]])){
idx <- df[[groupN]] == oi
tmp <- df[idx, ]
tmp <- tmp[order(tmp[[oN]], decreasing = decreasing),]
df[idx, ] <- tmp
}
df
}
tcga_gene <- read.csv('22.02.18.tcga_rowRanges.csv', row.names = 1)
tcga_gene['seqnames'] <- factor(tcga_gene$seqnames, levels = paste0('chr', c(1:22,'X', 'Y')))
tcga_gene <- tcga_gene[order(as.numeric(tcga_gene$seqnames)),]
tcga_gene <- f_df_order(tcga_gene, 'start', 'seqnames', F)
sce <- readRDS('~/upload/zl_liu/data/pca.rds')
[email protected] <- readRDS('../fig1/E_meta.data.rds')
tcga_gene <- tcga_gene[tcga_gene$external_gene_name %in% rownames(sce), c('external_gene_name', 'seqnames', 'start', 'end')]
write.table(x = tcga_gene, file = 'gencode_downsampled.txt', sep = '\t', row.names = F, col.names = F, quote = F)
第三步 导出 raw_counts_matrix 和 annotations_file
library(Seurat)
sce <- readRDS('~/upload/zl_liu/data/pca.rds')
[email protected] <- readRDS('../fig1/E_meta.data.rds')
counts_matrix = GetAssayData(seurat_obj, slot="counts")
saveRDS(round(counts_matrix, digits=3), "sc.10x.counts.matrix")
df <- [email protected]
df <- subset(df, cell_type_fig1A %in% c('Fibroblasts', 'Luminal', 'Endothelial'))
tmp <- subset(df, patient_id == 'patient1')['cell_type_fig1A']
write.table(x = tmp, file = 'annotations_file_p1.txt', sep = '\t', row.names = T, col.names = F, quote = F)
tmp <- subset(df, patient_id == 'patient3')['cell_type_fig1A']
write.table(x = tmp, file = 'annotations_file_p3.txt', sep = '\t', row.names = T, col.names = F, quote = F)
tmp <- subset(df, patient_id == 'patient4')['cell_type_fig1A']
write.table(x = tmp, file = 'annotations_file_p4.txt', sep = '\t', row.names = T, col.names = F, quote = F)
tmp <- subset(df, patient_id == 'patient5')['cell_type_fig1A']
write.table(x = tmp, file = 'annotations_file_p5.txt', sep = '\t', row.names = T, col.names = F, quote = F)
第四步 进行分析
fuck <- read.table('gencode_downsampled.txt', sep = '\t')
fuck <- fuck[!duplicated(fuck$V1),]
write.table(x = fuck, file = 'gencode_downsampled.txt', sep = '\t', row.names = F, col.names = F, quote = F)
infercnv_obj_p1 = CreateInfercnvObject(raw_counts_matrix=readRDS("sc.10x.counts.matrix"),
annotations_file="annotations_file_p1.txt",
delim="\t",
gene_order_file="gencode_downsampled.txt",
ref_group_names=c('Fibroblasts', 'Endothelial')
)
infercnv_obj_p3 = CreateInfercnvObject(raw_counts_matrix=readRDS("sc.10x.counts.matrix"),
annotations_file="annotations_file_p3.txt",
delim="\t",
gene_order_file="gencode_downsampled.txt",
ref_group_names=c('Fibroblasts', 'Endothelial')
)
infercnv_obj_p4 = CreateInfercnvObject(raw_counts_matrix=readRDS("sc.10x.counts.matrix"),
annotations_file="annotations_file_p4.txt",
delim="\t",
gene_order_file="gencode_downsampled.txt",
ref_group_names=c('Fibroblasts', 'Endothelial')
)
infercnv_obj_p5 = CreateInfercnvObject(raw_counts_matrix=readRDS("sc.10x.counts.matrix"),
annotations_file="annotations_file_p5.txt",
delim="\t",
gene_order_file="gencode_downsampled.txt",
ref_group_names=c('Fibroblasts', 'Endothelial')
)
infercnv_obj_p1 = infercnv::run(infercnv_obj_p1,
cutoff=0.1, # use 1 for smart-seq, 0.1 for 10x-genomics
out_dir="p1", # dir is auto-created for storing outputs
cluster_by_groups=F,
analysis_mode="subclusters",
tumor_subcluster_partition_method = 'random_trees',
denoise=TRUE,
HMM=TRUE,
num_threads=12)
saveRDS(infercnv_obj_p1, 'infercnv_obj_p1.rds')
system(paste("/opt/conda/bin/python3 /home/jovyan/upload/zl_liu/wecomchan.py", "'task p1 successfully completed'"), intern = T)
infercnv_obj_p3 = infercnv::run(infercnv_obj_p3,
cutoff=0.1, # use 1 for smart-seq, 0.1 for 10x-genomics
out_dir="p3", # dir is auto-created for storing outputs
cluster_by_groups=F,
analysis_mode="subclusters",
tumor_subcluster_partition_method = 'random_trees',
denoise=TRUE,
HMM=TRUE,
num_threads=12)
saveRDS(infercnv_obj_p3, 'infercnv_obj_p3.rds')
system(paste("/opt/conda/bin/python3 /home/jovyan/upload/zl_liu/wecomchan.py", "'task p3 successfully completed'"), intern = T)
infercnv_obj_p4 = infercnv::run(infercnv_obj_p4,
cutoff=0.1, # use 1 for smart-seq, 0.1 for 10x-genomics
out_dir="p4", # dir is auto-created for storing outputs
cluster_by_groups=F,
analysis_mode="subclusters",
tumor_subcluster_partition_method = 'random_trees',
denoise=TRUE,
HMM=TRUE,
num_threads=12)
saveRDS(infercnv_obj_p4, 'infercnv_obj_p4.rds')
system(paste("/opt/conda/bin/python3 /home/jovyan/upload/zl_liu/wecomchan.py", "'task p4 successfully completed'"), intern = T)
infercnv_obj_p5 = infercnv::run(infercnv_obj_p5,
cutoff=0.1, # use 1 for smart-seq, 0.1 for 10x-genomics
out_dir="p5", # dir is auto-created for storing outputs
cluster_by_groups=F,
analysis_mode="subclusters",
tumor_subcluster_partition_method = 'random_trees',
denoise=TRUE,
HMM=TRUE,
num_threads=12)
saveRDS(infercnv_obj_p5, 'infercnv_obj_p5.rds')
system(paste("/opt/conda/bin/python3 /home/jovyan/upload/zl_liu/wecomchan.py", "'task p5 successfully completed'"), intern = T)
Comments NOTHING