knitr::opts_chunk$set(echo = TRUE)
library("Seurat")
## Loading required package: ggplot2
## Loading required package: reshape2
## Loading required package: useful
## Loading required package: gridExtra
## Loading required package: gplots
## 
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
## 
##     lowess
## Loading required package: gdata
## gdata: read.xls support for 'XLS' (Excel 97-2004) files ENABLED.
## 
## gdata: read.xls support for 'XLSX' (Excel 2007+) files ENABLED.
## 
## Attaching package: 'gdata'
## The following object is masked from 'package:gridExtra':
## 
##     combine
## The following objects are masked from 'package:useful':
## 
##     left, right
## The following object is masked from 'package:stats':
## 
##     nobs
## The following object is masked from 'package:utils':
## 
##     object.size
## The following object is masked from 'package:base':
## 
##     startsWith
library("ggplot2")
barplotReAxis<-function(mp){
mp<-c(0,mp)
np<-c()
end<-length(mp)-1
for(i in 1:end){
  np.tmp<-(mp[i]+mp[i+1])/2
  np<-c(np,np.tmp)
}
np<-round(np)
names(np)=names(mp)[2:length(mp)]
return(np)
}
load("C:\\Users\\shicheng\\Downloads\\Barplot_Workspace.RData")
cell.labels <- substr(colnames(nbt.data),0,7)
cell.labels<-unlist(lapply(rownames(t.nbt.data),function(x) unlist(strsplit(x,"_"))[1]))
group.labels<-unlist(lapply(strsplit(cell.labels,split="BA"),function(x) x[[1]]))
input<-data.frame(GEX=t.nbt.data[,1],cell.labels=cell.labels,group.labels=group.labels)
input<-input[order(input$cell.labels),]
quantile(input$GEX)
##       0%      25%      50%      75%     100% 
## 0.000000 0.000000 0.000000 1.056400 6.808702
par(mar=c(3,3,2,4))
barplot(input$GEX,border=input$cell.labels)
legend("topright",legend=unique(input$group.labels),lty=rep(1,length(unique(input$group.labels))),col=as.numeric(as.factor(unique(input$group.labels))),bty="n",cex=0.5)

t.nbt.data<-t(nbt.data)
cell.labels <- substr(colnames(nbt.data),0,7)
cell.labels<-unlist(lapply(rownames(t.nbt.data),function(x) unlist(strsplit(x,"_"))[1]))
sort(table(cell.labels))
## cell.labels
## Ex7BA17 In3BA17  In2BA8 In5BA17 Ex8BA17 In7BA17 Ex2BA41  Ex3BA8  In7BA8 
##       1       1       2       2       3       3       4       4       4 
## Ex3BA21 In7BA21 Ex2BA22 Ex8BA22 Ex8BA41 In1BA17 In3BA22 In7BA10 Ex8BA21 
##       5       5       6       6       6       6       6       6       7 
## In2BA17 In2BA21 Ex6BA41  Ex8BA8 In5BA41 In2BA22 In4BA17 In4BA21 In5BA22 
##       7       7       8       8       8       9       9       9       9 
## Ex3BA22 Ex6BA17 Ex7BA41  In3BA8 In5BA10 In7BA22 Ex5BA17 Ex6BA22 In5BA21 
##      10      10      10      10      10      10      11      11      11 
## In2BA10 In8BA22 In3BA10 In3BA21 In4BA22 Ex3BA10 In6BA17 In8BA17  Ex4BA8 
##      12      12      13      14      14      15      16      16      17 
## Ex7BA22 Ex8BA10 In1BA21 In1BA22 In2BA41 In6BA10 In6BA22 Ex2BA21 In8BA21 
##      17      17      17      18      18      18      18      19      20 
##  Ex7BA8  In8BA8  In5BA8 Ex4BA22 In4BA10 Ex5BA22 In3BA41 In4BA41 In8BA10 
##      21      21      22      25      25      27      27      27      28 
## Ex7BA21 Ex6BA21 Ex2BA10 Ex5BA41  In1BA8 In7BA41  Ex2BA8 Ex4BA21 Ex6BA10 
##      29      31      32      34      34      34      35      35      37 
## Ex7BA10  In4BA8 In6BA21 In8BA41 Ex4BA41 In1BA10  Ex6BA8 In1BA41 Ex5BA21 
##      37      37      37      38      41      41      42      44      49 
## Ex4BA10 Ex1BA17 Ex5BA10  Ex5BA8  Ex1BA8 In6BA41  In6BA8 Ex3BA41 Ex1BA22 
##      56      62      64      65      70      73      77      84     131 
## Ex1BA10 Ex3BA17 Ex1BA21 Ex1BA41 
##     153     181     218     424
sort(table(group.labels))
## group.labels
##  Ex8  In2  In5  In7  In3  Ex2  Ex7  In4  In8  Ex6  In1  Ex4  In6  Ex5  Ex3 
##   47   55   62   62   71   96  115  121  135  139  160  174  239  250  299 
##  Ex1 
## 1058
sum(table(cell.labels)>30)
## [1] 30
table(cell.labels)[table(cell.labels)>30]
## cell.labels
## Ex1BA10 Ex1BA17 Ex1BA21 Ex1BA22 Ex1BA41  Ex1BA8 Ex2BA10  Ex2BA8 Ex3BA17 
##     153      62     218     131     424      70      32      35     181 
## Ex3BA41 Ex4BA10 Ex4BA21 Ex4BA41 Ex5BA10 Ex5BA21 Ex5BA41  Ex5BA8 Ex6BA10 
##      84      56      35      41      64      49      34      65      37 
## Ex6BA21  Ex6BA8 Ex7BA10 In1BA10 In1BA41  In1BA8  In4BA8 In6BA21 In6BA41 
##      31      42      37      41      44      34      37      37      73 
##  In6BA8 In7BA41 In8BA41 
##      77      34      38
In1<-c("In1BA10","In1BA8","In1BA21","In1BA22","In1BA41","In1BA17")
In2<-c("In2BA10","In2BA8","In2BA21","In2BA22","In2BA41","In2BA17")
In3<-c("In3BA10","In3BA8","In3BA21","In3BA22","In3BA41","In3BA17")
In4<-c("In4BA10","In4BA8","In4BA21","In4BA22","In4BA41","In4BA17")
In5<-c("In5BA10","In5BA8","In5BA21","In5BA22","In5BA41","In5BA17")
In6<-c("In6BA10","In6BA8","In6BA21","In6BA22","In6BA41","In6BA17")
In7<-c("In7BA10","In7BA8","In7BA21","In7BA22","In7BA41","In7BA17")
In8<-c("In8BA10","In8BA8","In8BA21","In8BA22","In8BA41","In8BA17")
Ex1<-c("Ex1BA10","Ex1BA8","Ex1BA21","Ex1BA22","Ex1BA41","Ex1BA17")
Ex2<-c("Ex2BA10","Ex2BA8","Ex2BA21","Ex2BA22","Ex2BA41","Ex2BA17")
Ex3<-c("Ex3BA10","Ex3BA8","Ex3BA21","Ex3BA22","Ex3BA41","Ex3BA17")
Ex4<-c("Ex4BA10","Ex4BA8","Ex4BA21","Ex4BA22","Ex4BA41","Ex4BA17")
Ex5<-c("Ex5BA10","Ex5BA8","Ex5BA21","Ex5BA22","Ex5BA41","Ex5BA17")
Ex6<-c("Ex6BA10","Ex6BA8","Ex6BA21","Ex6BA22","Ex6BA41","Ex6BA17")
Ex7<-c("Ex7BA10","Ex7BA8","Ex7BA21","Ex7BA22","Ex7BA41","Ex7BA17")
Ex8<-c("Ex8BA10","Ex8BA8","Ex8BA21","Ex8BA22","Ex8BA41","Ex8BA17")
groups <- factor(cell.labels,levels=c(In1,In2,In3,In4,In5,In6,In7,In8,Ex1,Ex2,Ex3,Ex4,Ex5,Ex6,Ex7,Ex8))
NewGroup<-names(sort(table(groups),decreasing=T))[1:30]
bar.len<-c(cumsum(table(groups)))
newaxis<-barplotReAxis(bar.len)
nbt.data<-data.frame(nbt.data)
nGenes <- nrow(nbt.data)
coverage <- colSums(nbt.data)/nGenes
ord <- order(groups)

fig.width=8

bar.positions <- barplot(coverage[ord],border=groups[ord],xaxt='n',ylab="Counts per gene",beside=T,space=10)
axis(side=1,at=c(bar.positions[newaxis[match(NewGroup,names(newaxis))]]),labels=NewGroup,tick=FALSE,las=2,cex.axis=0.7)

fig.width=16

bar.positions <- barplot(coverage[ord],border=groups[ord],xaxt='n',ylab="Counts per gene",beside=T,space=10)
axis(side=1,at=c(bar.positions[newaxis[match(NewGroup,names(newaxis))]]),labels=NewGroup,tick=FALSE,las=2,cex.axis=0.7)

fig.width=64

bar.positions <- barplot(coverage[ord],border=groups[ord],xaxt='n',ylab="Counts per gene",beside=T,space=10)
axis(side=1,at=c(bar.positions[newaxis[match(NewGroup,names(newaxis))]]),labels=NewGroup,tick=FALSE,las=2,cex.axis=0.7)