Expressed gene count / cell plot

샘플별로, cell 별로 몇개의 유전자가 발현하고 어느정도의 발현량을 갖는지 확인하기 위한 코드를 짜봤다.

bulkRNAseq data를 다루던 버릇이 남아서 그럴수도…

cellranger results에서 web-summary.html을 열어보면 mean expressed gene count가 나오긴 한다. 그래도 좀 더 자세히 알아보고 싶어서..

#Data loading
## seurat 에서 읽어오는 것 기준
rawdata <- as.data.frame(GetAssayData(SEURAT.OBJ,slot="counts")) #genes(row) X cells(col) raw counts table
rawdata <- rawdata[rowSums(rawdata) > 0,] #모든 cell에서 발현이 0인 유전자는 거르고 시작

#bulk RNA-seq data랑 비교해보기 위해 기준을 맞춤 
# 0~1/ 1~5/ 5~10/ 10~20/ 20~50/ 50~100/ 200이상

# count 최대값구하기
max(rawdata)
786
## 대충 보고 적당한 선까지 만들자

#stacked barplot으로 그리기 위해 input data frame 만들어줌
bar_input <- data.frame(colnames(rawdata))
bar_input$count.0 <- 0 #0이상 - 1미만
bar_input$count.1 <- 0 #1-5
bar_input$count.2 <- 0 #5-10
bar_input$count.3 <- 0 #10-20
bar_input$count.4 <- 0 #20-50
bar_input$count.5 <- 0 #50-100
bar_input$count.6 <- 0 #100-200
bar_input$count.7 <- 0 #200-500
bar_input$count.8 <- 0 #500-1000
bar_input$count.9 <- 0 #1000-5000
bar_input$count.10 <- 0 #5000-max

for(i in 1:ncol(rawdata)){ # cell 
    for(j in 1:nrow(rawdata)){ # gene
        counts <- rawdata[j,i]
        #print(paste0(i,':',j,'->',counts))
        if(counts >=0 & counts <1) bar_input$count.0[i] = bar_input$count.0[i]+1
        else if(counts >=1 & counts <5) bar_input$count.1[i] = bar_input$count.1[i]+1
        else if(counts >=5 & counts <10) bar_input$count.2[i] = bar_input$count.2[i]+1
        else if(counts >=10 & counts <20) bar_input$count.3[i] = bar_input$count.3[i]+1
        else if(counts >=20 & counts <50) bar_input$count.4[i] = bar_input$count.4[i]+1
        else if(counts >=50 & counts <100) bar_input$count.5[i] = bar_input$count.5[i]+1
        else if(counts >=100 & counts <200) bar_input$count.6[i] = bar_input$count.6[i]+1
        else if(counts >=200 & counts <500) bar_input$count.7[i] = bar_input$count.7[i]+1
        else if(counts >=500 & counts <1000) bar_input$count.8[i] = bar_input$count.8[i]+1
        else if(counts >=1000 & counts <5000) bar_input$count.9[i] = bar_input$count.9[i]+1
        else bar_input$count.10[i] = bar_input$count.10[i]+1
    }
}
colnames(bar_input)[1] <- 'cell'
#결과저장
write.table(bar_input ,"counts_bar.txt", sep='\t')
####
#Bar plot 그려보자
#bar_input <- read.table("counts_bar.txt", sep='\t',header=TRUE, row.names=1)
library(RColorBrewer)
sequential <- c("#ECEFF1","#FDD835","#CDDC39","#8BC34A","#4CAF50","#1976D2","#0D47A1","#1A237E","#1A237E","#1A237E","#1A237E")

pdf("bar.pdf")
barplot(t(bar_input[,2:12]), ## 거꾸로 그리고 싶으면 bar_input[,12:2]
 names.arg = bar_input$cell, # x-axis labels
 cex.names = 0.1, # makes x-axis labels small enough to show all
 col = sequential, # colors
 xlab = "cells",
 ylab = "genes",
 ylim = c(0,26500), #전체 유전자 갯수로 y축 기준 세움
 width = 0.5, space=0)
dev.off()

Result

내 데이터 상으로는 cell당 많으면 6000개 정도의 유전자를 발현하고 나머지 유전자는 발현값이 0인듯 하다. 모든 cell에서 발현이 0인 유전자는 제외했으니 남아있는 유전자들은 어떤 cell에서는 1이라도 발현이 되긴한다는 뜻. single cell data가 얼마나 sparse한지 확 와닿는군. 약 1000개가 안되는 유전자가 5 이상의 값을 갖는듯.

bulk RNA와 비교해보기 위해 기준을 동일하게 적용했는데, 좀 단위를 낮춰서 세분화시켜 컬러링하는것도 좋을듯.

뭐 샘플별로 plot을 나눠서 그린다던지, upgrade가능할듯. 지금은 좀 귀찮아서 여기까지.