Chao.py

library(Seurat)
library(Signac)
library(jsonlite)

tissue_map <- c("2305121", "2305122", "2305123", "2305124", 
                "V4-A", "V4-B", "V4-C", "V4-D")
names(tissue_map) <- c(2305121, 2305122, 2305123, 2305124, seq(456, 459))
cells <- read.table("/brahms/srivastavaa/chao/visium/visium-v1.txt")$V1

coordinates <- read.table("/brahms/srivastavaa/chao/visium/visium-v1_coordinates.txt")
names(coordinates) <- c("cb", "row", "col")
head(coordinates)

{
  plots <- lapply(names(tissue_map), function(frag.id) {
    well.id <- as.character(tissue_map[frag.id])
    
    ncounts <- CountFragments(
      fragments = paste0("/brahms/srivastavaa/chao/fragments/", frag.id,".bed.gz"),
      cells = cells
    )
    
    sum(ncounts$frequency_count)
  })
  cts <- unlist(plots)
  names(cts) <- names(tissue_map)
  
  df <- data.frame(cts)
  df$exp <- rownames(df)
  
  ggplot(df, aes(x=exp, y=cts, color=exp))+
    geom_bar(stat = "identity")+
    theme_classic()
}

{
  plots <- lapply(names(tissue_map), function(frag.id) {
    well.id <- as.character(tissue_map[frag.id])
    
    ncounts <- CountFragments(
      fragments = paste0("/brahms/srivastavaa/chao/fragments/", frag.id,".bed.gz"),
      cells = cells
    )
    
    frag.counts <- ncounts[, c(1,2)]
    rownames(frag.counts) <- frag.counts$CB
    frag.counts$CB <- NULL

    frag.counts$row <- unlist(lapply(rownames(frag.counts), function(cb){
      coordinates[coordinates$cb == cb, "row"]
    }))
    
    frag.counts$col <- unlist(lapply(rownames(frag.counts), function(cb){
      coordinates[coordinates$cb == cb, "col"]
    }))
    
    ggplot(frag.counts, aes(x=row, y=col, color=frequency_count)) +
      geom_point(size=0.1) +
      scale_colour_gradient2() +
      ggtitle(frag.id)
  })
  wrap_plots(plots, nrow = 2)
}
{
  plots <- lapply(names(tissue_map), function(frag.id) {
    well.id <- as.character(tissue_map[frag.id])
    
    ncounts <- CountFragments(
      fragments = paste0("/brahms/srivastavaa/chao/fragments/", frag.id,".bed.gz"),
      cells = cells
    )
    
    frag.counts <- ncounts[, c(1,2)]
    rownames(frag.counts) <- frag.counts$CB
    frag.counts$CB <- NULL
    frag.counts$condition <- frag.id
    print(sum(frag.counts$frequency_count > 0) / 4992)
    
    frag.counts
  })
  plots <- do.call(rbind, plots)
  ggplot(plots, aes(x=condition, y=log1p(frequency_count), 
                    color=condition)) +
    geom_jitter() +
    geom_violin() +
    theme_classic()
  
  #plots$nz <- plots$frequency_count > 0
  #ggplot(plots, aes(x=condition, y=nz, 
  #                  color=condition)) +
  #  geom_histogram() +
  #  theme_classic()
}

{
  #df <- lapply(names(tissue_map)[-14], function(frag.id) {
  #df <- lapply(names(tissue_map)[16:19], function(frag.id) {
  df <- lapply(names(tissue_map), function(frag.id) {
    well.id <- as.character(tissue_map[frag.id])
    
    ncounts <- CountFragments(
      fragments = paste0("/brahms/srivastavaa/chao/fragments/", frag.id,".bed.gz"),
      cells = cells
    )
    
    tissue <- fromJSON(paste0("/brahms/srivastavaa/chao/visium/", well.id,".json"))$oligo
    tissue <- tissue[!is.na(tissue$tissue), ]
    tissue$row <- tissue$row + 1
    tissue$col <- tissue$col + 1
    
    tissue$swp <- paste0(tissue$col, "_", tissue$row)
    coordinates$swp <- paste0(coordinates$row, "_", coordinates$col)
    
    cb.in <- coordinates[coordinates$swp %in% tissue$swp, "cb"]
    cb.out <- setdiff(cells, cb.in)
    
    df.in <- ncounts[ncounts$CB %in% cb.in, ]
    df.out <- ncounts[ncounts$CB %in% cb.out, ]
    
    sum.in <- sum(df.in$frequency_count)
    sum.out <- sum(df.out$frequency_count)
    fraction <- sum.in * 100.0 / (sum.out + sum.in)
    fraction <- format(round(fraction, 2), nsmall = 2)
    
    df.in$kind <- "tissue"
    df.out$kind <- "background"
    
    df.in$experiment <- paste0(frag.id, "-", sum.out + sum.in, "-", fraction, "%")
    df.out$experiment <- paste0(frag.id, "-", sum.out + sum.in, "-", fraction, "%")
    
    rbind(df.in, df.out)
  })
  
  df.all <- do.call(rbind, df)
  ggplot(df.all, aes(y = log1p(frequency_count),
                     x = experiment, 
                     fill = kind)) + 
    geom_violin() +
    #geom_jitter() +
    #geom_bar(position="dodge")
    theme_classic()
  #(p1 | p2) / (p3 | p4)
}