hpliner
commented
4 years ago Hi, I just ran through the code with the mouse kidney data used in the tutorial with results that looked fine at first glance (i.e. top 100 genes put in a go enrichment turned up kidney)...
A couple of notes
- the values are going to be very small because of the normalization - its designed more for looking at variation in activity scores among cells rather than comparing between genes.
- even in the best circumstances, there's not going to be anything near a 1:1 correspondence between expression and gene activity - while some genes seem to vary strongly with gene activity score, other don't - likely some quirk of different modes of regulation
Hope this helps
I tried to run the gene activity score pipeline using the exact code: https://cole-trapnell-lab.github.io/cicero-release/docs_m3/#cicero-gene-activity-scores
But I can't seem to get it working on 10X's 10K PBMC data set or my own B cell data set. In both cases, the activity score looks very low for some key marker genes that I know are expressed. Do you guys have some result plots of gene activity score that you can share based on some example data set?
FYI, these are the code I used to generate and plot the scores: ` hg19 <- read.table("../hg19.chrom.sizes.txt") conns <- run_cicero(cicero_cds, hg19)
gene_anno <- rtracklayer::readGFF("/somewhere/gencode.v19.annotation.gtf") gene_anno$chromosome <- gene_anno$seqid gene_anno$gene <- gene_anno$gene_id gene_anno$transcript <- gene_anno$transcript_id gene_anno$symbol <- gene_anno$gene_name pos <- subset(gene_anno, strand == "+") pos <- pos[order(pos$start),] pos <- pos[!duplicated(pos$transcript),] pos$end <- pos$start + 1 neg <- subset(gene_anno, strand == "-") neg <- neg[order(neg$start, decreasing = TRUE),] neg <- neg[!duplicated(neg$transcript),] neg$start <- neg$end - 1 gene_annotation_sub <- rbind(pos, neg)
gene_annotation_sub <- gene_annotation_sub[,c("chromosome", "start", "end", "symbol")] names(gene_annotation_sub)[4] <- "gene" input_cds <- annotate_cds_by_site(input_cds, gene_annotation_sub) unnorm_ga <- build_gene_activity_matrix(input_cds, conns) unnorm_ga <- unnorm_ga[!Matrix::rowSums(unnorm_ga) == 0, !Matrix::colSums(unnorm_ga) == 0] num_genes <- pData(input_cds)$num_genes_expressed names(num_genes) <- row.names(pData(input_cds)) cicero_gene_activities <- normalize_gene_activities(unnorm_ga, num_genes)
genes = c("CD8A", "CD4", "NKG7", "MS4A1", "PRDM1", "CD19", "IRF8")
for(gene in genes){ x1 <- cicero_gene_activities[gene,] pData(input_cds)$marker_x <- x1 cairo_pdf( paste0("gene_activityscore", gene, ".pdf"), width=4, height=4) p <- plot_cells(input_cds, show_trajectory_graph=F, label_branch_points=F, color_cells_by="marker_x") print(p) dev.off() } `