R/visualization.r

In JuanXie19/TCRlineageClust:

#' Clonotype distribution plot 
#'
#' This function shows how each clonotype is distributed on the UMAP.
#'
#' @param combinedData The Seurat object with clonotype information.
#' @param I The index of the clonotype which you'd like to plot its distribution
#'
#' @details The combinedData must contain UMAP coordinates for cells and cell barcodes saved as 'cell.id'. 
#' It must also contain the clonotype information for each cell. 
#'
#' @import Seurat
#' @import dplyr
#' @import ggplot2
#' @import tools
#' 
#' @examples
#' TCR <-read.csv('C:/TCR-seq/Gang/GSE158896_RAW/combinedTCR_2cells.csv',header=T')
#' load('seurat.rda')  # load seurat object with UMAP coordinates for cells
#' CombinedData <-getCombinedDataSet(TCR,seurat.data)
#' plot(CombinedData,1)  # plot the first clonotype

if(FALSE){
setMethod(
	f = 'plot',
	signature = c(object='CombinedDataSet',I='numeric'),
	definition = function(object,I,...){
		UMAP1 <- object@seurat@reductions$umap@cell.embeddings[,1]
		UMAP2 <- object@seurat@reductions$umap@cell.embeddings[,2]
		G <- object@seurat$clusters
		G.df <- data.frame(UMAP1,UMAP2,G,cdr3=object@seurat$cdr3)
		
		CLONE <- unique(G.df$cdr3)

		
		INDEX <- which(G.df$cdr3==CLONE[I])

		G.df.highlight <- G.df[INDEX,]
			
		ggplot(G.df,aes(UMAP1,UMAP2,color=G))+geom_point(alpha=0.6,shape=1,size=0.6)+
				geom_point(data=G.df.highlight,aes(UMAP1,UMAP2),size=3,shape=17)+
				theme_bw()+
				labs(title=paste0('Cell_',length(INDEX),'_cluster_',length(unique(G.df.highlight$G))))
		if(length(INDEX)<3){
			print('too few cells')
			}			
	}
)

}

	.plotCombinedDataSet <- function(x,I,...){
		UMAP1 <- x@seurat@reductions$umap@cell.embeddings[,1]
		UMAP2 <- x@seurat@reductions$umap@cell.embeddings[,2]
		G <- x@seurat$clusters
		G.df <- data.frame(UMAP1,UMAP2,G,cdr3=x$cdr3)
	
		CLONE <- unique(G.df$cdr3)

		
		INDEX <- which(G.df$cdr3==CLONE[I])

		G.df.highlight <- G.df[INDEX,]
			
		ggplot(G.df,aes(UMAP1,UMAP2,color=G))+geom_point(alpha=0.6,shape=1,size=0.6)+
			geom_point(data=G.df.highlight,aes(UMAP1,UMAP2),size=3,shape=17)+
			theme_bw()+
			labs(title=paste0('Cell_',length(INDEX),'_cluster_',length(unique(G.df.highlight$G))))
		if(length(INDEX)>=3){
			ggsave(paste0(dir.out,CLONE[i],'.pdf'),width=10,height=6)
		}else{
			print('too few cells')
		}			
	}

#'
#' Overlay lineages
#'
#' This function displays all clonotype lineages.
#' @param data The \code{TrajectoryDataSet}object
#' @import slingshot
#' @import Seurat
#' @import dplyr
#' @import tidyverse
#'
#' @export
#'
#' @return A plot of overlaied MST
#' 
#' @examples
#' TrajectoryData <- getClonotypeLineages(CombinedDataSet)
#' plot(TrajectoryData)

if(FALSE){
setMethod(
	f = 'plot',
	signature = c('TrajectoryDataSet'),
	definition = function(object,...){
		df <- object@lineages
		
		t <- lapply(df,function(x) x%>%group_split(Lineage)) # some MST may contain multiple lineages, here split
		t <- unlist(t,recursive=FALSE)

		t <- lapply(t,function(x) column_to_rownames(x, var = "Cluster"))
		
		temp <-lapply(t, function(x) x[,1:2])
		
		UMAP1_min <- floor(min(sapply(temp,function(x) min(x$UMAP_1))))-1
		UMAP1_max <- ceiling(max(sapply(temp,function(x) max(x$UMAP_2))))+1
		UMAP2_min <- floor(min(sapply(temp,function(x) min(x$UMAP_2))))-1
		UMAP2_max <- ceiling(max(sapply(temp,function(x) max(x$UMAP_2))))+1
		
		plot(temp[[1]]$UMAP_1,temp[[1]]$UMAP_2,type='l',xlab='UMAP1',col='blue',ylab='UMAP2',xlim=c(UMAP1_min,UMAP1_max),ylim=c(UMAP2_min,UMAP2_max),pch=19,lwd=1.5)
		text(temp[[1]]$UMAP_1,temp[[1]]$UMAP_2,labels=rownames(temp[[1]]),cex=0.5)
	
		for (i in 2:length(temp)){
			lines(temp[[i]]$UMAP_1,temp[[i]]$UMAP_2,type='l',col="blue",pch=19,cex=0.5,lwd=1.5)
			text(temp[[i]]$UMAP_1,temp[[i]]$UMAP_2,labels=rownames(temp[[i]]),cex=0.5)
		}	
	}
)
}


	.plotTrajectoryDataSet <- function(x,...){
		df <- x@lineages
		
		t <- lapply(df,function(x) x%>%group_split(Lineage)) # some MST may contain multiple lineages, here split
		t <- unlist(t,recursive=FALSE)

		t <- lapply(t,function(x) column_to_rownames(x, var = "Cluster"))
		
		temp <-lapply(t, function(x) x[,1:2])
		
		UMAP1_min <- floor(min(sapply(temp,function(x) min(x$UMAP_1))))-1
		UMAP1_max <- ceiling(max(sapply(temp,function(x) max(x$UMAP_2))))+1
		UMAP2_min <- floor(min(sapply(temp,function(x) min(x$UMAP_2))))-1
		UMAP2_max <- ceiling(max(sapply(temp,function(x) max(x$UMAP_2))))+1
		
		plot(temp[[1]]$UMAP_1,temp[[1]]$UMAP_2,type='l',xlab='UMAP1',col='blue',ylab='UMAP2',xlim=c(UMAP1_min,UMAP1_max),ylim=c(UMAP2_min,UMAP2_max),pch=19,lwd=1.5)
		text(temp[[1]]$UMAP_1,temp[[1]]$UMAP_2,labels=rownames(temp[[1]]),cex=0.5)
	
		for (i in 2:length(temp)){
			lines(temp[[i]]$UMAP_1,temp[[i]]$UMAP_2,type='l',col="blue",pch=19,cex=0.5,lwd=1.5)
			text(temp[[i]]$UMAP_1,temp[[i]]$UMAP_2,labels=rownames(temp[[i]]),cex=0.5)
		}	

	}


#' This function displays clusters of clonotype lineages.
#' @param object The \code{ClusterInput}object
#' @param TSCLUST The clustering results from \code{clustLineages}
#' @param I which cluster you'd like to plot
#'
#' @import slingshot
#' @import Seurat
#' @import dplyr
#'
#' @export
#'
#' @return A plot of overlaied MST
#' 
#' @examples
#' TrajectoryData <- getClonotypeLineages(CombinedDataSet,start.clus = NULL, end.clus = NULL, dist.method = 'simple', use.median = TRUE)
#' ClusterInput <-getClusteringInput(TrajectoryData)
#' p.dtw <-clustLineages(ClusterInput,type='p',seed=123456,iter.max=20,k=7L)
#' plot(ClusterInput,p.dtw)
#'

if(FALSE){
setMethod(
	f = 'plot',
	signature = c(object='ClusterInput', TSCLUST='ANY',I='numeric'),
	definition = function(object,TSCLUST,I,...){
		df <- object@lineages
		clusters <-TSCLUST@cluster
		## plot the trajectories within each cluster
		
		UMAP1_min <- floor(min(sapply(df,function(x) min(x$UMAP_1))))-1
		UMAP1_max <- ceiling(max(sapply(df,function(x) max(x$UMAP_2))))+1
		UMAP2_min <- floor(min(sapply(df,function(x) min(x$UMAP_2))))-1
		UMAP2_max <- ceiling(max(sapply(df,function(x) max(x$UMAP_2))))+1
			
		
		INDEX <-which(TSCLUST@cluster==I)
			
			
		plot(df[[INDEX[1]]][,1:2],type='l',xlab = "UMAP1",col="blue", ylab = "UMAP2",xlim=c(UMAP1_min,UMAP1_max),ylim=c(UMAP2_min,UMAP2_max),pch=19,lwd=1.5)
		text(df[[INDEX[1]]][,1],df[[INDEX[1]]][,2],labels=rownames(df[[INDEX[1]]]),cex=0.5)
		for (j in 2:length(INDEX)){
			lines(df[[INDEX[j]]][,1],df[[INDEX[j]]][,2],type='l',col="blue",pch=19,cex=0.5,lwd=1.5)
			text(df[[INDEX[j]]][,1],df[[INDEX[j]]][,2],labels=rownames(df[[INDEX[j]]]),cex=0.5)
			}
		dev.off()

	}

)
}


	.plotClusterInput <-function(x,I,...){
		df <- ClusterInput@lineages
		
		
		
		clusters <-TSCLUST@cluster
		## plot the trajectories within each cluster
		
			
		UMAP1_min <- floor(min(sapply(temp,function(x) min(x$UMAP_1))))-1
		UMAP1_max <- ceiling(max(sapply(temp,function(x) max(x$UMAP_2))))+1
		UMAP2_min <- floor(min(sapply(temp,function(x) min(x$UMAP_2))))-1
		UMAP2_max <- ceiling(max(sapply(temp,function(x) max(x$UMAP_2))))+1
			
		
		INDEX <-which(TSCLUST@cluster==I)
			
			
		plot(df[[INDEX[1]]][,1:2],type='l',xlab = "UMAP1",col="blue", ylab = "UMAP2",xlim=c(UMAP1_min,UMAP1_max),ylim=c(UMAP2_min,UMAP2_max),pch=19,lwd=1.5)
		text(df[[INDEX[1]]][,1],df[[INDEX[1]]][,2],labels=rownames(df[[INDEX[1]]]),cex=0.5)
		for (j in 2:length(INDEX)){
			lines(df[[INDEX[j]]][,1],df[[INDEX[j]]][,2],type='l',col="blue",pch=19,cex=0.5,lwd=1.5)
			text(df[[INDEX[j]]][,1],df[[INDEX[j]]][,2],labels=rownames(df[[INDEX[j]]]),cex=0.5)
		}
		dev.off()

	}