In a study published in Nature Communications in 2019 (Alcala et al. Nat Commun 2019), we performed a molecular characterization of the rare understudied lung neuroendocrine neoplasms (lung NENs) that allowed us to identify clinically relevant molecular groups. To complement this study, we have submitted a data note to journal GigaScience where we firstly describe the preprocessing and the quality control steps performed on data generated by Alcala et al. and propose the first molecular map of lung NENs (pan-LNEN map) using the Uniform Manifold Approximation and Projection (UMAP) on RNA-sequencing data. This map integrates homogenized data from multiple molecular studies on lung NENs (Peifer et al. Nat Genet 2012; Fernandez-Cuesta et al. Nat Commun 2014; George et al. Nature 2015; George et al. Nat Commun 2018; Alcala et al. Nat Commun 2019, Laddha et al. Cancer research 2019) with the aim of increasing the sample size to better understand this rare disease. The map was published on the TumorMap interactive portal.
This notebook complements the data note by providing the code necessary to reproduce the pan-LNEN molecular map. We also propose to assess the quality of the molecular map by: i) validating the biological hypotheses proposed in the previous studies, ii) exploring the preservation of the samples' neighborhood in the map, and iii) validating the preservation of spatial auto-correlations in the map.
This section allows to set up the working environment, install the R packages, and upload the scripts and data sets needed for the analyses.
Install CRAN packages
install.packages(c("ade4", # Principal Component Analysis
"umap", # Uniform manifold analysis projection
"FNN", # K-Neirest Neighbours
"ape", # for Moran index
"reshape2", # Basic tools
"ggpubr", # Layout
"gridExtra", # Layout
"eulerr", # Layout for UMAP
"viridis", # Layout
"foreach", # Parallel computing
"doParallel" # Parallel computing
))
Install Bioconductor packages
install.packages("BiocManager")
BiocManager::install("DESeq2")
library(umap)
library(ade4)
library(FNN)
library(reshape2)
library(ape)
library(ggpubr)
library(foreach)
library(doParallel)
library(gridExtra)
library(eulerr)
library(viridis)
library(DESeq2)
library(RColorBrewer)
sessionInfo()
ls
Note that to avoid running the notebook entirely, the R environment allowing to generate all the results included in the notebook has been saved in the file R_environment.Rdata
and can be loaded by uncommenting the following line.
#load('R_environment.Rdata')
The scripts used to perform and evaluate the dimensionality reduction of the pan-LNEN RNA-seq data are available in the following GitHub repository:
The DRMetrics repository contains (see readme for additional details):
- A folder named "scripts" containing the R scripts enabling the assessment of the quality of low-dimensional representations
- A folder named "data" containing the data required to reproduce the results presented in the GigaScience paper:
- read_counts_all.txt.zip: homogenized gene expression quantification from RNA-seq data from several published manuscripts (Peifer et al. Nat Genet 2012; Fernandez-Cuesta et al. Nat Commun 2014; George et al. Nature 2015; George et al. Nat Commun 2018; Alcala et al. Nat Commun 2019, Laddha et al. 2019). The data are stored as a matrix of raw read counts of 59,607 genes (rows) from the 239 samples (columns, including 1 technical replicate).
- Attributes.txt.zip: main clinical features and the relevant genes (columns) for the 238 unique individuals (rows) reported in the previous manuscripts.
- Coords_umap_nn238.txt.zip: UMAP coordinates (columns) associated with the map shared on the TumorMap interactive portal, for each of 238 unique individuals (rows).
- ref_annot.gtf.gene_spans.zip: genome annotation associated to the GRCh38 human reference genome, version 33, from bundle CTAT from the 6th of April 2020
- hg38_Gencode_V33.bed.zip and hg38_Gencode_V33_merged.bed.zip: bed files used in the RNA-Seq data preprocessing (see below)
Loading the main functions needed to perform the analyses:
source("/DRMetrics/scripts/DR_quality_metrics.r")
Loading the RNA-Seq and clinical data:
The raw RNA-seq fastq files were preprocessed following four nextflow pipelines available in the IARCbioinfo GitHub repository. Note that the pipelines were run using singularity but docker and conda profiles could also be used.
Step 1: we performed RNA sequencing mapping and quality controls using the pipeline IARCbioinfo/RNAseq-nf by running the following command:
nextflow run iarcbioinfo/RNAseq-nf -r v2.3 -profile singularity --input_folder folder_with_fastq_files --output_folder out_RNAseq-nf-2.3 --ref_folder ref_genome.fa.star.idx/ --ref ref_genome.fa --gtf ref_annot.gtf --bed hg38_Gencode_v33.bed --fastq_ext fastq --STAR_mapqUnique 60 --cutadapt
Where folder_with_fastq_files
is a folder containing all the fastq files, ref_genome.fa.star.idx
is a folder with the genome reference files for software STAR, ref_annot.gtf
the annotation file and hg38_Gencode_v33.bed
a bed file with list of intervals for further annotations by RSeQC. The bed file was generated using the UCSC table browser and the genome reference files, the annotation file should be associated to the GRCh38 human reference genome, version 33, from bundle CTAT from the 6th of April 2020 (located in the folder GRCh38_gencode_v33_CTAT_lib_Apr062020.plug-n-play/ctat_genome_lib_build_dir/
).
Step 2: we refined alignments of the BAM files using the pipeline IARCbioinfo/abra-nf by running the following command:
nextflow run iarcbioinfo/abra-nf -r v3.0 -profile singularity --bam_folder out_RNAseq-nf-2.3/BAM --output_folder out_abra-nf-3.0 --ref ref_genome.fa --gtf ref_annot.gtf --bed hg38_Gencode_v33_merged.bed --junctions --rna
The BAM files used as input for ABRA are the outputs of the first alignment step. The reference and annotation files used are the same as the ones used in step 1 with an exception for the bed file which has been modified so that no overlapping intervals remain using bedtools:
sort -k1,1 -k2,2n hg38_Gencode_v33.bed | bedtools merge -i stdin > hg38_Gencode_v33_merged.bed
Step 3: we performed base quality score recalibration using the pipeline IARCbioinfo/BQSR-nf by running the following command:
nextflow run iarcbioinfo/BQSR-nf -r v1.1 -profile singularity --input_folder out_abra-nf-3.0/ --ref ref_genome.fa --snp_vcf dbsnp_146.hg38.vcf.gz --indel_vcf Mills_and_1000G_gold_standard.indels.hg38.vcf.gz --output_folder out_bqsr-nf-1.1
Where the vcf files come from the GATK bundle hg38.
Step 4: we performed gene expression quantification using the pipeline IARCbioinfo/RNAseq-transcript-nf by running the following command:
nextflow run iarcbioinfo/RNAseq-transcript-nf -r v2.1 -profile singularity --input_file input-transcript.txt
--gtf ref_annot.gtf --output_folder out_RNAseq-transcript-nf-2.1
The file input-transcript.txt
contains in rows, the samples ID, paths to BAM files, and read lengths per sample. The latter were computed using the post-trimming fastqc outputs of the first step by computing the average of the two fastq files read length values and rounding it.
For our analysis, we merged the samples raw read counts generated by the fourth step in a single matrix. In the following code cell, we load the read counts associated to all samples as well as the clinical data associated.
Attributes = read.table(unzip('DRMetrics/data/Attributes.txt.zip')[1], sep = '\t', header = T)
Attributes[1:5,1:10]
read_counts = as.matrix(read.table(unzip('DRMetrics/data/read_counts_all.txt.zip')[1], header = T, row.names=1))
read_counts[1:5,1:10]
We removed read counts associated with genes located on sex and mitochondrial chromosomes. After this filtering 56 597 remained.
# Read gene annotations
genespans = read.table(unzip('DRMetrics/data/ref_annot.gtf.gene_spans.zip')[1],stringsAsFactors = F)
dim(genespans)
rownames(genespans) = genespans[,1]
# Remove genes on sex and mitochondrial chromosomes
genenames = rownames(read_counts)
summary(genenames %in% rownames(genespans))
genes.sex = genespans[genespans$V2 %in% c("chrM", "chrX", "chrY"), ]
genes.nosex = genespans[!(genespans$V2 %in% c("chrM", "chrX", "chrY")), ]
read_counts = read_counts[which(rownames(read_counts) %in% rownames(genes.nosex)), ]
dim(read_counts)
Data pre-processing:
The raw read counts needed to be normalized; we used the variance stabilizing transformation from the DESeq2 R package.
colData = as.matrix(colnames(read_counts))
DESeq_object = DESeqDataSetFromMatrix(countData=read_counts,colData = colData,design = ~1)
VST = varianceStabilizingTransformation(DESeq_object)
VST = assay(VST)
dim(VST)
VST[1:5,1:5]
For each analysis, the most variable genes, explaining 50% of the total variance in variance-stabilized read counts, were selected (6,390 out of the 56,597 initial genes).
rvdm = apply(VST,1,var)
ntopdm = rev(which( cumsum( sort(rvdm,decreasing = T) )/sum(rvdm) <= 0.5 ))[1]
selectdm = order(rvdm, decreasing = TRUE)[seq_len(ntopdm)]
vst50 = t(VST[selectdm,])
vst50[1:5,1:5]
dim(vst50)
We kept only one replicate per individual (the column Sample_ID in table Attributes, excluding sample S00716).
vst50 = vst50[as.character(Attributes$Sample_ID),]
dim(vst50)
The analysis of the gene expression matrix is complex due to its high dimensionality. Dimensionality reduction methods are thus essential for hypothesis generation and interpretations. In this section, we present firstly the projections resulting from the principal component analysis (PCA) method and secondly a two-dimensional projection resulting from the Uniform Manifold Analysis Projection (UMAP) algorithm.
In the previously published datasets that we integrated, different molecular groups were identified: carcinoids A1, A2, and B (Alcala et al. Nat Commun 2019), and clusters LC1, LC2 and LC3 (Laddha et al. 2019)-that we found here corresponding respectively to carcinoids A1, B and A2 (see low dimensional projections below)-, LCNEC/TypeI and LCNEC/TypeII (George mun 2018), and SCLC (George et al. Nat Commun 2015). In order to distinguish the six molecular groups, we considered five principal components when evaluating the PCA results.
We used the dudi.pca function to perform a PCA on the lung NEN samples:
pca5D = dudi.pca(vst50, center = T , scale = T, scannf = F , nf=5 )$li
head(pca5D)
However, to visualize the samples in the reduced space generated by the PCA , we needed to consider only two dimensions. The next code cell represents the lung NEN samples in several two dimensional spaces based on the different combination of the first five principal components.
distinctive_cols <- c("#9E0142","#E41A1C","#E95D47","#FF8904","#FDBE6E","#FFD422","#88CFA4","#47A265","#397EB8","#00CCFF","#91569B","#D01C8B","#CCCCCC")
Axis_list = list(c(2,1), c(3,1), c(4,1), c(5,1), c(3,2), c(4,2), c(5,2), c(4,3), c(5,3), c(5,4))
plist <- list()
for(i in 1:length(Axis_list)){
x_index = Axis_list[[i]][1]
y_index = Axis_list[[i]][2]
xlab = paste("PC", x_index)
ylab = paste("PC", y_index)
c_data = data.frame("axisX"= pca5D[,x_index], "axisY" =pca5D[,y_index] )
p = ggplot(c_data, aes(x = axisX , y = axisY, color = Attributes$Molecular_clusters)) + geom_point(size=2, alpha=0.8) + scale_colour_manual(values=distinctive_cols)
if (x_index == 2 & y_index == 1){
plist[[i]] = p + theme( legend.position="right", legend.title= element_blank()) + labs(x = xlab, y = ylab)
}else{
plist[[i]] = p + theme( legend.position="none", legend.title= element_blank()) + labs(x = xlab, y = ylab)
}
}
ggarrange(plist[[1]])
grid.arrange(grobs = lapply(2:10, function(i){plist[[i]]}),
layout_matrix = rbind(c(2:4),c(5:7),c(8:10)), nrow =2)
Another way to visualize data in a two dimensional space is to use UMAP. The UMAP algorithm is a non-linear dimensionality reduction technique based on topology. We applied UMAP to the RNA-seq pan-LNEN data using the umap R package. Using this method we aimed at visualizing the molecular clusters identified in previous studies. UMAP doesn't necessarily preserve the global structure of the data and one important parameter for this is n_neighbors
. Small and high values of n_neighbors
favor the preservation of local structures and global structures, respectively. The pan-LNEN TumorMap was hence generated fixing the n_neighbors
parameter to 238, which is the total number of samples. The other parameters were set to their default values.
UMAP_TM = umap(vst50, n_neighbors = 238)$layout
colnames(UMAP_TM) = c("UMAP1","UMAP2")
UMAP_TM=as.data.frame(UMAP_TM)
Because UMAP is not deterministic, in order to reproduce the exact results presented in the GigaScience paper we upload in the following cell the UMAP coordinates used to create the pan-LNENs TumorMap (available in the IARCbioinfo/DRMetrics github repository):
# Note that the next line replace the UMAP coordinates calculated in the previous cell with our pre-computed coordinates. You should comment it if you are running a customized analysis.
UMAP_TM = read.table(unzip('DRMetrics/data/Coords_umap_nn238.txt.zip')[1], sep = '\t', header = T, row.names=1)
UMAP_TM <- cbind(UMAP_TM, 'Molecular_clusters' = Attributes$Molecular_clusters)
ggplot(UMAP_TM, aes(x = UMAP1 , y = UMAP2, color = Molecular_clusters))+
geom_point(size=4, alpha=0.8)+ scale_colour_manual(values=distinctive_cols)+
labs(x = "UMAP dimension 1" , y = "UMAP dimension 2", color = "")+ theme(legend.position = "bottom")+
guides(col = guide_legend( nrow = 4))
head(UMAP_TM)
Although non-linear dimensionality reduction methods such as UMAP lack interpretability in comparison to PCA, we attempted to show that this simple layout captures the main biological findings highlighted in the previous studies (Fernandez-Cuesta et al. Nat Commun 2014; George et al. Nature 2015; George et al. Nat Commun 2019; Alcala et al. Nat Commun 2019; Laddha et al. Cancer research 2019).
In previous studies, some subgroups of lung NENs have been identified as having the molecular features of a histopathological type that is different from their initial histopathological classification. We checked whether our tumor map captures this phenomenon by comparing the distances from these discordant tumors to their molecular and histopathological groups, respectively. The following code chunk corresponds to the code used to perform these comparisons using a one-sided Wilcoxon's test.
clusters_distances_F <- function(Coords, evals, ref1, ref2, alt_hyp){
# Coords: a data frame with the following columns: i) the samples coordinates on the dimensionality reduction axes, and ii) the molecular cluster associated to each sample
# evals: the molecular cluster to analyse
# ref1: the supposed molecular cluster(s)
# ref2: the supposed histopathological group(s)
# alt_hyp: the alternative hypothesis for the wilcoxon test
Coords_evals = Coords[Coords$Molecular_clusters==evals,]
T1 = Coords[which(Coords$Molecular_clusters %in% ref1),]
T1_centroid = colMeans(T1[,1:2])
dist_1 = sqrt( (Coords_evals[,1] - T1_centroid[1])^2 + ( Coords_evals[,2] - T1_centroid [2])^2 )
T2 = Coords[which(Coords$Molecular_clusters %in% ref2),]
T2_centroid = colMeans(T2[,1:2])
dist_2 = sqrt( (Coords_evals[,1] - T2_centroid[1])^2 + ( Coords_evals[,2] - T2_centroid [2])^2 )
return(wilcox.test(dist_1 , dist_2, alt_hyp)$p.value)
}
As an example, we can consider the LCNEC/SCLC-like samples identified by George et al. in 2019 as being LCNEC samples that have the molecular profile of the SCLC samples. Using the function above, we showed that in the pan-LNEN TumorMap, these LCNEC/SCLC-like samples were indeed significantly closer to the SCLC samples than the LCNEC samples:
wilcoxon_res = clusters_distances_F(UMAP_TM, "LCNEC/SCLC-like","SCLC/SCLC-like",c("LCNEC/TypeI", "LCNEC/TypeII"), 'less')
paste('Wilcoxon p-value = ', format(wilcoxon_res, digits = 2))
We evaluated the preservation of the samples' neighborhood using the sequence difference view (SD) dissimilarity metric proposed by Martins et al. in 2015. This measure compares samples' neighborhoods in two different spaces, $$$\small
D_1$ and
Parameter
The SD metric firstly evaluates, for each point, its ranking in the two neighborhoods (
The following code chunk uses the compute_SD
function to compute the SD values associated to each sample for different values of
distpca5D = as.matrix(dist(pca5D, method = "euclidean", diag = TRUE, upper = TRUE))
distpca5D[1:5,1:5]
k_values = seq(2,238,5)
SD_pca5D = compute_SD_allSamples(distRef = distpca5D, List_projection = list("UMAP_TM" = UMAP_TM[,1:2]), k_values = k_values,colnames_res_df = "UMAP-nn238")
head(SD_pca5D)
The SD_map_f
function is then called to represent the lung NEN samples in the UMAP projection and show, with a color gradient, the SD values averaged over the
sd_map_lnen = SD_map_f(SD_pca5D , UMAP_TM, "bottom")
sd_map_res = sd_map_lnen[[1]]
head(sd_map_res)
sd_map_lnen[[2]]
The above representation allows the identification of the samples with the lowest sample neighborhood preservation between UMAP and PCA-5D. We chose to explore the samples' neighborhood of S02340 which have the highest SD metric. This sample is an atypical sample that clusters with the Carcinoid-A1 cluster on the UMAP projection. To explain this high SD value, we computed the euclidean distances between each sample and S02340 in the two spaces.
plot_dist <- function(coords,projection_coords,focal_sample){
dist1 = as.matrix(dist(coords[, 1:(ncol(coords)-1)], method = "euclidian", diag = TRUE, upper = TRUE))
Sample_names = rownames(coords)
dist1 = data.frame("Sample_ID" = rownames(dist1)[order(dist1[,focal_sample])], "rank" = 1:length(rownames(dist1)) )
merged_df = merge(cbind(Sample_ID = rownames(projection_coords), projection_coords), dist1, by= "Sample_ID")
colnames(merged_df)[2:3] = c("V1","V2")
proj_dist = ggplot(merged_df, aes(x=V1, y=V2, color=rank, label=merged_df$Sample_ID)) +
geom_point(size=2, alpha =.8) + scale_color_distiller(palette = "Spectral") +
geom_point(aes(merged_df$V1[which(merged_df$Sample_ID == focal_sample)],
merged_df$V2[which(merged_df$Sample_ID == focal_sample)] ),
size=2, col = "black") +
geom_text(aes(merged_df$V1[ merged_df$Sample_ID== focal_sample],
merged_df$V2[ merged_df$Sample_ID== focal_sample],
label=merged_df$Sample_ID[merged_df$Sample_ID==focal_sample]),
col = 'black' , hjust=0, vjust=0 ) +
labs(title="", y="dim2", x="dim1") + theme( legend.position = "bottom")
return(proj_dist)
}
p1 <- plot_dist(UMAP_TM,UMAP_TM,"S02340")
p2 <- plot_dist(pca5D,UMAP_TM,"S02340")
grid.arrange(grobs = list(p1,p2), layout_matrix = rbind(c(1,2)), nrow =1)
In the above figure, the left and right panels represent the UMAP projection of lung NEN samples, with a color gradient corresponding to the rank of each sample in the neighborhoods of S02340 in UMAP and PCA-5D projections, respectively. In the UMAP projection (left panel), the nearest neighbors of S02340 belong to the carcinoids cluster (top cluster). However, in the PCA projection (right panel) S02340 is also close to LCNEC samples, more specifically the LCNEC Type I samples (bottom cluster), suggesting that the UMAP projection for that sample might be misleading.
Exploring the SD metric for each sample allows us to identify samples with questionable projections.
Using the SD metric, we also evaluated the global preservation of the samples neighborhoods of different dimensionality reductions: i) PCA-5D: the projection based on the five first principal components of the PCA, ii) PCA-2D: the projection based on the two first principal components of the PCA, iii) projections generated by UMAP based on different values of the n_neighbors
parameter (15 (default value), 30, 90, 180 and 238).
n_neighborsL <- c(15,30,90,180,238)
List_projection <- list("pca2D" = pca5D[1:2], "pca5D" = pca5D)
for (j in 1:length(n_neighborsL)){
umap_c = umap(vst50, n_neighbors = n_neighborsL[j])$layout
List_projection[[j+2]] <- umap_c
names(List_projection)[j+2] <- paste("umap_nn", n_neighborsL[j])
}
We computed the SD metric for each sample and each projection while varying the parameter n_iter
parameter below) but this number can be increased by the user. In the GigaScience paper, n_iter
was fixed to 1000 to generate Figure 4.
distRef = as.matrix(dist(vst50, method = "euclidean", diag = TRUE, upper = TRUE))
distRef[1:5,1:5]
k_values = seq(2,238,25)
n_neighborsL = c(15,30,90,180,238)
colnames_res_df = c('PCA-2D', 'PCA-5D', paste0('UMAP-nn',n_neighborsL) )
n_iter = 10
no_cores = 2
cl = makeCluster(no_cores)
registerDoParallel(cl)
SD_means <- foreach(i=1:n_iter, .packages = c("umap","doParallel"))%dopar%{
List_projection = list("PCA-2D" = pca5D[,1:2], "PCA-5D" = pca5D)
for (j in 1:length(n_neighborsL)){
List_projection[[j+2]] = umap(vst50, n_neighbors = n_neighborsL[j])$layout
names(List_projection)[j+2] = paste("umap_nn", n_neighborsL[j])
}
SD_comp = compute_SD_allSamples(distRef, List_projection, k_values, colnames_res_df)
# For each k compute the mean value over all samples
SD_k_mean = sapply(1:length(colnames_res_df), function(j) aggregate(SD_comp[,colnames_res_df[j]],list(SD_comp$K),mean)$x)
colnames(SD_k_mean) = colnames_res_df
rownames(SD_k_mean) = k_values
return(SD_k_mean)
}
stopCluster(cl)
For each of the 10 replicate iterations (elements of list SD_means), for each dimensionality reduction technique (columns), we computed the mean SD value across samples for different *$\small k$ *(rows):
#mean SD as function of k for the first iteration
head(SD_means[[1]])
Mean_SD = melt( apply(simplify2array(SD_means), 1:2, mean), varnames = c("k","DRM"), value.name = "mean" )
sd_SD = melt( apply(simplify2array(SD_means), 1:2, sd), varnames = c("k","DRM"), value.name = "sd" )
Main_SD_DF = cbind(Mean_SD, sd = sd_SD$sd)
head(Main_SD_DF)
p = ggplot(Main_SD_DF, aes(x = k, y = mean, color = DRM)) +
geom_line() + geom_point()+
scale_color_viridis(discrete=TRUE) +
geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd))+
labs(title = "", y = "mean(SD)", x = "k")
p
To facilitate the visualization and to compare the sample neighborhood preservation of the different projections, we normalized the mean SD values so that the SD of the PCA-2D is 1 and that of the PCA-5D is 0, as described in the GigaScience paper (equation 1).
SD_means_norm = lapply( SD_means, function(x) (x[,2] - x)/(x[,2] - x[,1] ) )
Mean_SD_norm = melt( apply(simplify2array(SD_means_norm), 1:2, mean), varnames = c("k","DRM"), value.name = "mean" )
sd_SD_norm = melt( apply(simplify2array(SD_means_norm), 1:2, sd), varnames = c("k","DRM"), value.name = "sd" )
Main_SD_norm_DF = cbind(Mean_SD_norm, sd = sd_SD_norm$sd)
theme_set(theme_bw())
ggplot(Main_SD_norm_DF, aes(x = k, y = mean, color = DRM)) +
geom_line() + geom_point()+
scale_color_viridis(discrete=TRUE) +
geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd))+
labs(title = "", y = "mean(SD')", x = "k") +
theme( axis.title.x=element_text(size=20), # X axis title
axis.title.y=element_text(size=20), # Y axis title
axis.text.x=element_text(size=16), # X axis text
axis.text.y=element_text(size=16),
legend.text = element_text(size = 20) ,
legend.title = element_blank())
The above figure shows that varying the n_neighbors
parameter has an impact on the sample neighborhood preservation. While lower values of n_neighbors
improve the preservation of the local sample neighborhood (low SD for small n_neighbors
improve the global preservation of the neighborhood (low SD large n_neighbors=238
is a good compromise between local and global preservation, having preservations at least as good as PCA-2D for most
We can evaluate the spatial auto-correlation for focal genes using the Moran index. This metric ranges from -1 (inverse spatial autocorrelation) to 1 (perfect spatial autocorrelation). For a given gene, the index considers the similarity between the gene expression levels of samples in the
n_genes = ncol(vst50)
k_step = seq(2,238,25)
MI_GeneExprMat = moran_I_knn(vst50[,1:n_genes], vst50,k_step)
MI_UMAP = moran_I_knn(vst50[,1:n_genes], UMAP_TM[,1:2],k_step)
MI_PCA5D = moran_I_knn(vst50[,1:n_genes], pca5D,k_step)
head(MI_UMAP)
We computed first the average of the Moran's indexes over
MI_PCA5D_order = order(colMeans(MI_PCA5D[,1,]), decreasing=T)
MI_UMAP_order = order(colMeans(MI_UMAP[,1,]), decreasing=T)
MI_GeneExprMat_order = order(colMeans(MI_GeneExprMat[,1,]), decreasing=T)
PCA5D_GeneExprMat_overlap = sapply(1:n_genes, function(i) mean( MI_PCA5D_order[1:i] %in% MI_GeneExprMat_order[1:i] ) )
PCA5D_UMAP_overlap = sapply(1:n_genes, function(i) mean( MI_PCA5D_order[1:i] %in% MI_UMAP_order[1:i] ) )
UMAP_GeneExprMat_overlap = sapply(1:n_genes, function(i) mean( MI_UMAP_order[1:i] %in% MI_GeneExprMat_order[1:i] ) )
PCA5D_UMAP_GeneExprMat_overlap = sapply(1:n_genes, function(i) mean( MI_PCA5D_order[1:i] %in% MI_GeneExprMat_order[1:i] & MI_PCA5D_order[1:i] %in% MI_UMAP_order[1:i] ) )
Expected = sapply(1:n_genes, function(i) (i/n_genes) )
cols <- c("Expected" = "#FDE725FF", "PCA 5D - GeneExprMat" = "#440154FF","PCA 5D - UMAP 238 - GeneExprMat" = "#5DC863FF", "PCA 5D - UMAP 238" = "#21908CFF", "UMAP 238 - GeneExprMat" = "#3B528BFF")
ggplot()+
geom_point(aes(x=seq(1:n_genes), y=PCA5D_GeneExprMat_overlap , colour="PCA 5D - GeneExprMat"), size = 0.1)+
geom_point(aes(seq(1:n_genes), y=UMAP_GeneExprMat_overlap, colour="UMAP 238 - GeneExprMat"), size = 0.1)+
geom_point(aes(seq(1:n_genes), y=PCA5D_UMAP_overlap , colour="PCA 5D - UMAP 238"), size = 0.1)+
geom_point(aes(seq(1:n_genes), y=PCA5D_UMAP_GeneExprMat_overlap, colour="PCA 5D - UMAP 238 - GeneExprMat"), size = 0.1)+
labs(title="", y="Overlapping %", x="N", caption = "") +
scale_colour_manual(values=cols, labels =c("Expected", "PCA 5D - GeneExprMat","PCA 5D - UMAP 238 - GeneExprMat", "PCA 5D - UMAP 238", "UMAP 238 - GeneExprMat")) +
theme( legend.position = "bottom")+
geom_line(aes(seq(1:n_genes), y=Expected, colour="Expected"), size = 1)+
geom_vline(xintercept = 1000, linetype="dashed")
We finally represented the Euler diagram of the gene sets corresponding to the top 1000 Moran indexes:
n_genes_sel = 1000
abc = PCA5D_UMAP_GeneExprMat_overlap[n_genes_sel]*n_genes_sel
ab = PCA5D_UMAP_overlap[n_genes_sel]*n_genes_sel
ac = PCA5D_GeneExprMat_overlap[n_genes_sel]*n_genes_sel
bc = UMAP_GeneExprMat_overlap[n_genes_sel]*n_genes_sel
Venn <- c(c("A" = n_genes_sel - ((ab -abc) + abc + (ac -abc)) ,
"B" = n_genes_sel - ((ab -abc) + abc + (bc -abc)),
"C" = n_genes_sel - ((ac -abc) + abc + (bc -abc)),
"A&B" = ab -abc,
"A&C" = ac -abc,
"B&C" = bc -abc ,
"A&B&C" = abc))
PlotEuler <- plot(euler(Venn), fills = c("#5DC863FF", "#FDE725FF", "#21908CFF"), quantities = list(cex = 1), alpha= .7, cex = 2, legend = list(labels = c("PCA 5D", "UMAP nn = 238", "GenesExprMat"), alpha=.7 ,cex =1.2))
PlotEuler
The Euler diagram above represents the overlap between the top 1000 genes in the three lists and shows that 88.8% of the genes are in common. This result shows that the UMAP projection preserves genes spatial auto-correlation.
save.image('R_environment.Rdata')
This notebook was exported from https://nextjournal.com/a/LjZ9gwQSLfU5EMBgXRFrJ?change-id=CiJotpJrhqgchpzJGYM3Qx
{:article
{:settings {:subtitle? false, :numbered? true, :authors? true},
:nodes
{"005f73af-71f3-4108-82bf-060ab3fb9010"
{:compute-ref #uuid "789f6d90-a9e7-4057-b51c-58819b8d30d8",
:exec-duration 275,
:id "005f73af-71f3-4108-82bf-060ab3fb9010",
:kind "code",
:locked? false,
:output-log-lines {},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"0c6fa5c3-4404-4579-a8d2-46f5edc2218c"
{:id "0c6fa5c3-4404-4579-a8d2-46f5edc2218c",
:kind "github-repository",
:ref "4fb3527a37b3bb48f4b88e059a72d5ced6e3fa3a"},
"13d01cef-053c-4588-a8ed-da1685c608d0"
{:compute-ref #uuid "31a6b5a7-24a1-4214-8d95-bdfa69b0eb10",
:exec-duration 214,
:id "13d01cef-053c-4588-a8ed-da1685c608d0",
:kind "code",
:locked? false,
:output-log-lines {},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"1ac3cd60-b22d-4e6c-be18-5e8a1c1b6529"
{:compute-ref #uuid "7d7d70e1-12b7-48f2-bd56-71991346a275",
:exec-duration 2548,
:id "1ac3cd60-b22d-4e6c-be18-5e8a1c1b6529",
:kind "code",
:locked? false,
:output-log-lines {},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"1c82ce51-9c71-4b29-a0d8-849a199dadc3"
{:compute-ref #uuid "dce8b76b-dc28-47ff-95e8-dea0ece80088",
:exec-duration 10143,
:id "1c82ce51-9c71-4b29-a0d8-849a199dadc3",
:kind "code",
:locked? false,
:output-log-lines {:stdout 156},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"],
:stdout-collapsed? true},
"40b1fcb8-eeb7-4c3c-9a00-1da6a4f69f93"
{:compute-ref #uuid "cd00cc98-eba2-436e-9a19-83c43ccdd03e",
:exec-duration 697,
:id "40b1fcb8-eeb7-4c3c-9a00-1da6a4f69f93",
:kind "code",
:locked? false,
:output-log-lines {},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"410d0763-12d4-41f5-92d5-7150cadecc38"
{:compute-ref #uuid "94208b6f-9b13-4675-956d-540bc5519767",
:exec-duration 14191,
:id "410d0763-12d4-41f5-92d5-7150cadecc38",
:kind "code",
:locked? false,
:output-log-lines {},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"45682072-c84b-4a42-b632-5aee00034ade"
{:compute-ref #uuid "58ed5e5f-6afa-4c0c-b021-6c2a000e7050",
:exec-duration 55917,
:id "45682072-c84b-4a42-b632-5aee00034ade",
:kind "code",
:locked? true,
:output-log-lines {:stdout 71},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"],
:stdout-collapsed? true},
"467f7e8a-9e06-438d-81da-0fa880ef1297"
{:compute-ref #uuid "bcc1fffc-f13a-4c11-968f-dc4151962465",
:exec-duration 6150,
:id "467f7e8a-9e06-438d-81da-0fa880ef1297",
:kind "code",
:locked? false,
:output-log-lines {},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"4dcc91d9-74ce-4890-a8ba-0e2a498fdf6b"
{:compute-ref #uuid "3796ec64-591c-4ae3-a7f2-59c3cfd55cab",
:exec-duration 662,
:id "4dcc91d9-74ce-4890-a8ba-0e2a498fdf6b",
:kind "code",
:output-log-lines {:stdout 3},
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"509e262c-1ace-4cf1-b841-f016f30815cd"
{:compute-ref #uuid "4a9334a3-8d59-41e5-ba07-044be63f1e53",
:exec-duration 7973,
:id "509e262c-1ace-4cf1-b841-f016f30815cd",
:kind "code",
:locked? false,
:output-log-lines {},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"509f01f5-a364-4863-ab0a-5dd7be5d70e9"
{:compute-ref #uuid "ddf10f89-7c8e-4544-82fc-e69a68517e25",
:exec-duration 539368,
:id "509f01f5-a364-4863-ab0a-5dd7be5d70e9",
:kind "code",
:locked? true,
:output-log-lines {:stdout 1935},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"],
:stdout-collapsed? true},
"5a7f58c2-aaf6-4c11-a5fb-2354339e0c67"
{:compute-ref #uuid "b70bf8f4-c0e8-4c89-8118-d6f842d58922",
:exec-duration 253,
:id "5a7f58c2-aaf6-4c11-a5fb-2354339e0c67",
:kind "code",
:locked? false,
:output-log-lines {},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"6414702a-ec41-40cb-a026-10c15e1542a6"
{:compute-ref #uuid "8d1090ef-b9bd-46d9-aaff-6fb860703385",
:exec-duration 374,
:id "6414702a-ec41-40cb-a026-10c15e1542a6",
:kind "code",
:locked? false,
:output-log-lines {:stdout 2},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"66373f61-0971-4db0-9192-626600962918"
{:compute-ref #uuid "7a559477-d3b2-4b1c-ae20-2826d43f7f6e",
:exec-duration 395,
:id "66373f61-0971-4db0-9192-626600962918",
:kind "code",
:locked? false,
:output-log-lines {:stdout 2},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"74ba1b1d-7dc6-4582-9425-acf47845f01c"
{:compute-ref #uuid "c76d4e7d-b3e4-42d6-87b9-e01c3a6de409",
:exec-duration 20075,
:id "74ba1b1d-7dc6-4582-9425-acf47845f01c",
:kind "code",
:output-log-lines {},
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"7700f555-a7f5-45c8-aa0c-00d416412ede"
{:compute-ref #uuid "fac2de99-a2f0-46f0-8215-6497bf7fa2e4",
:exec-duration 113015,
:id "7700f555-a7f5-45c8-aa0c-00d416412ede",
:kind "code",
:locked? false,
:output-log-lines {:stdout 2},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"77c28fb5-2eb9-42e7-b56d-fdde93527d0e"
{:compute-ref #uuid "93287198-dc3c-4a2e-80fd-ff2e8e0e0489",
:exec-duration 1047,
:id "77c28fb5-2eb9-42e7-b56d-fdde93527d0e",
:kind "code",
:output-log-lines {:stdout 32},
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"797d97b0-dcbc-4135-896c-83f3185e0f18"
{:compute-ref #uuid "88515cae-51da-4774-82fa-bf790290ab85",
:exec-duration 5044,
:id "797d97b0-dcbc-4135-896c-83f3185e0f18",
:kind "code",
:locked? false,
:output-log-lines {},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"7a6270fd-36e0-4962-9e7f-70e1071568e3"
{:compute-ref #uuid "0b43780b-767b-45e8-ad77-29a258de598f",
:exec-duration 174451,
:id "7a6270fd-36e0-4962-9e7f-70e1071568e3",
:kind "code",
:locked? false,
:output-log-lines {},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"85b8d09e-d204-4b1b-8371-94279acb9f0a"
{:compute-ref #uuid "d1e04456-0109-4f94-97a9-d63e2ff9c75c",
:exec-duration 1915,
:id "85b8d09e-d204-4b1b-8371-94279acb9f0a",
:kind "code",
:locked? false,
:output-log-lines {},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"9f8ae330-19eb-47f6-a2f7-bb6ef9e96f2e"
{:compute-ref #uuid "8aba8102-3c21-4357-940a-a9abc38696bf",
:exec-duration 1822,
:id "9f8ae330-19eb-47f6-a2f7-bb6ef9e96f2e",
:kind "code",
:locked? false,
:output-log-lines {:stdout 2},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"acdf9fd1-15a4-403a-9946-2a2b4df8ee03"
{:compute-ref #uuid "3ef6cfed-cee1-4d3c-b54d-a2cea1a9294b",
:exec-duration 980,
:id "acdf9fd1-15a4-403a-9946-2a2b4df8ee03",
:kind "code",
:locked? false,
:output-log-lines {},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"ae19c779-93b2-4221-86c7-0198423848ad"
{:compute-ref #uuid "b071340b-d8b0-4a04-bc32-3c478ea0beda",
:exec-duration 449,
:id "ae19c779-93b2-4221-86c7-0198423848ad",
:kind "code",
:locked? false,
:output-log-lines {},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"c9a14654-44c5-482d-a55f-a88c5037397c"
{:compute-ref #uuid "434995f3-c214-4574-af89-6e8ae48287ff",
:exec-duration 581,
:id "c9a14654-44c5-482d-a55f-a88c5037397c",
:kind "code",
:locked? false,
:output-log-lines {},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"d73c09e3-b33f-48a7-9bde-ba9a90c8bb69"
{:compute-ref #uuid "9afe16b2-20fe-4a96-8502-4c8fea4cb203",
:exec-duration 151,
:id "d73c09e3-b33f-48a7-9bde-ba9a90c8bb69",
:kind "code",
:locked? false,
:output-log-lines {},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"e51d3446-16ca-4c64-a58f-66b5219318ae"
{:compute-ref #uuid "203fbfe6-3074-43ca-8a85-854a2e2ce866",
:exec-duration 475,
:id "e51d3446-16ca-4c64-a58f-66b5219318ae",
:kind "code",
:locked? false,
:output-log-lines {},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"e51ebb08-cde7-4c31-aafb-939b40c24b82"
{:compute-ref #uuid "6b3d4bd5-910a-4031-b05e-11a07b4acf64",
:exec-duration 1314,
:id "e51ebb08-cde7-4c31-aafb-939b40c24b82",
:kind "code",
:locked? false,
:output-log-lines {:stdout 5},
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"eaf11b0d-66d5-4547-aba2-0cb76239323a"
{:compute-ref #uuid "f1fb6490-5b8e-4fe8-a420-9e2a2e49c45b",
:exec-duration 388,
:id "eaf11b0d-66d5-4547-aba2-0cb76239323a",
:kind "code",
:locked? false,
:output-log-lines {},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"f152ae54-239a-4d2a-a0ba-1a1a6665e7b7"
{:compute-ref #uuid "0ae0fb32-5ace-44a3-94c1-707f95ee5e8f",
:exec-duration 6133,
:id "f152ae54-239a-4d2a-a0ba-1a1a6665e7b7",
:kind "code",
:locked? false,
:output-log-lines {},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"],
:stdout-collapsed? false},
"f44a00e9-f28b-48e1-96f2-82271d746bad"
{:compute-ref #uuid "b984041f-71a5-404e-89c6-d3e3fcfbe047",
:exec-duration 643825,
:id "f44a00e9-f28b-48e1-96f2-82271d746bad",
:kind "code",
:locked? false,
:output-log-lines {:stdout 2},
:refs (),
:runtime [:runtime "fbb04fb9-e412-4a90-a464-b3950395ffd6"]},
"fbb04fb9-e412-4a90-a464-b3950395ffd6"
{:environment
[:environment
{:article/nextjournal.id
#uuid "02c15b4d-a3ef-4b83-b807-84e3d2f75ba3",
:change/nextjournal.id
#uuid "5ece4061-5872-4aac-902c-8354fc3e9f24",
:node/id "fbb04fb9-e412-4a90-a464-b3950395ffd6"}],
:environment? true,
:id "fbb04fb9-e412-4a90-a464-b3950395ffd6",
:kind "runtime",
:language "r",
:name "R+Libs+DESeq2",
:resources {:machine-type "n1-highmem-4"},
:type :nextjournal,
:docker/environment-image
"docker.nextjournal.com/environment@sha256:ac6986cbc3c4610eb0275bbeebc329b42b69e7f92f135a6191d63a4ba9b45b95",
:runtime/checkpoint? false,
:runtime/mounts
[{:src [:node "0c6fa5c3-4404-4579-a8d2-46f5edc2218c"],
:dest "/DRMetrics"}]}},
:nextjournal/id #uuid "02c15b4d-a3ef-4b83-b807-84e3d2f75ba3",
:article/change
{:nextjournal/id #uuid "5ed8fea5-7093-467b-86ed-435c324227e5"}}}