DataDive.Rmd

---
title: "Imperial Data Dive: single cell alignment project"
author: "Nathan Skene"
date: "01/07/2019"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
R.utils::sourceDirectory("R")
library(tidyverse)
```

## Notes:

Tasic vs Zeisel2015:
  - Some cell types are not shared by both datasets, identifying these will be interesting challenge
  - Many cells match to multiple cell types

## Which datasets?

The following datasets will be provided

* Mouse single cell pair of datasets: Tasic vs Zeisel 2015
* Human vs mouse nuclei: Allan Human MTG vs Allan Mouse Visual (Tasic)
* Intron vs Exon: Allan MTG
* Two additional datasets: Zeisel 2018 & controls from ours

The datasets will be provided as a list. The expression matrix will be provided in sparse matrix format. An annotation table will accompany detailing the associated cell types. Only genes shared across all datasets will be used. Mouse genes will be mapped to human using the 1:1 package.

## Dataset: Zeisel 2015 (mouse somatosensory + hippocampus)

A function already exists in the EWCE package to load this dataset: 

```{r pressure, echo=FALSE}
library(devtools)
install_github("nathanskene/ewce")
library(EWCE)
linURL = "goo.gl/r5Y24y"
path = "expression_mRNA_17-Aug-2014.txt"
download.file(linURL, destfile=path)
zeisel2015 = EWCE::load.linnarsson.sct.data(path)
zeisel2015$exp = Matrix::Matrix(zeisel2015$exp)
file.remove(path)
zeisel2015$annot = zeisel2015$annot %>% dplyr::rename(cellID=cell_id,celltype=level2class) %>% dplyr::select(cellID,celltype)
save(tasic,file="Output/zeisel2015.rda")
```

## Dataset: Tasic (mouse visual cortex)

The Tasic data originally comes from here: http://casestudies.brain-map.org/celltax#section_explorea

It is associated with this paper: https://www.nature.com/articles/nn.4216

The Hemberg group downloaded the data and put it on AWS, which is what we download here (https://github.com/hemberg-lab/scRNA.seq.datasets/blob/master/bash/tasic.sh). 

```{r }
tasic = get_tasic_data()
tasic$annot = tasic$annot %>% dplyr::rename(cellID=long_name,celltype=Tasic_et_al_2016_label) %>% dplyr::select(cellID,celltype)
save(tasic,file="Output/tasic.rda")
```

## Mapping between Tasic and Zeisel2015

This is given in Supplementary Table 17 of the Tasic paper: https://www.nature.com/articles/nn.4216

https://media.nature.com/original/nature-assets/neuro/journal/v19/n2/extref/nn.4216-S1.pdf

```{r }
library("openxlsx")
mapping_tasic_zeisel2015 = read.xlsx("Data/zeisel_tasic_mapping.xlsx")
mapping_tasic_zeisel2015$comparison = "tasic - zeisel2015"
save(mapping_tasic_zeisel2015,file="Output/mapping_tasic_zeisel2015.rda")

check_if_all_celltypes_mapped(annotCTs=tasic$annot$celltype,mappingCTs=mapping_tasic_zeisel2015$Tasic)
check_if_all_celltypes_mapped(annotCTs=zeisel2015$annot$celltype,mappingCTs=mapping_tasic_zeisel2015$Zeisel)
```

## Intron vs Exon data

Data comes from here: http://celltypes.brain-map.org/download#transcriptomics

```{r }
library(data.table)
linURL = "http://celltypes.brain-map.org/api/v2/well_known_file_download/694416044"
path = "Allan_MTG.zip"
download.file(linURL, destfile=path)
unzip(path)
file.remove(path)

allanMTG_intron = load_allan_exp_matrix_with_hgnc_symbols(path="human_MTG_2018-06-14_intron-matrix.csv")
allanMTG_exon   = load_allan_exp_matrix_with_hgnc_symbols(path="human_MTG_2018-06-14_exon-matrix.csv")

# Check all intron cell types have corresponding exon cell types
sum(!unique(allanMTG_exon$annot$celltype) %in% unique(allanMTG_intron$annot$celltype))

# Delete the files
lapply(list(list.files(pattern="human_MTG_2018-06-14")),FUN=file.remove)

# Create the mapping file
mapping_allanMTG_intron_exon = data.frame(allanMTG_intron=unique(allanMTG_exon$annot$celltype),allanMTG_exon=unique(allanMTG_exon$annot$celltype))
mapping_allanMTG_intron_exon$comparison = "allanMTG_intron - allanMTG_exon"
save(mapping_allanMTG_intron_exon,file="Output/mapping_allanMTG_intron_exon.rda")

# Drop genes which are not expressed in one of the cell types
unexpressed_intron = Matrix::rowSums(allanMTG_intron$exp)==0
unexpressed_exon   = Matrix::rowSums(allanMTG_exon$exp)==0 
unexpressed        = unexpressed_intron | unexpressed_exon
allanMTG_intron$exp = allanMTG_intron$exp[!unexpressed,]
allanMTG_exon$exp   = allanMTG_exon$exp[!unexpressed,]

# Save
save(allanMTG_intron,file="Output/allanMTG_intron.rda")
save(allanMTG_exon,file="Output/allanMTG_exon.rda")
```

## Allan Mouse to Human mapping

A large single cell dataset from mouse visual cortex is available. While there is no authoratitive mapping available, it should form a good comparison set (and the cell type names are relatively informative). The dataset has 15,413 cells downloadable from here: http://celltypes.brain-map.org/download#transcriptomics

```{r }
library(data.table)
linURL = "http://celltypes.brain-map.org/api/v2/well_known_file_download/694413985"
path = "Allan_Mouse_V1.zip"
download.file(linURL, destfile=path)
unzip(path)
file.remove(path)

Allan_Mouse_V1_intron = load_allan_exp_matrix_with_hgnc_symbols(path="mouse_VISp_2018-06-14_intron-matrix.csv",prelim="mouse_VISp_2018-06-14")
Allan_Mouse_V1_exon   = load_allan_exp_matrix_with_hgnc_symbols(path="mouse_VISp_2018-06-14_exon-matrix.csv",prelim="mouse_VISp_2018-06-14")
Allan_Mouse_V1        = Allan_Mouse_V1_intron
Allan_Mouse_V1$exp    = Allan_Mouse_V1_intron$exp + Allan_Mouse_V1_exon$exp

# Delete the files
lapply(list(list.files(pattern="mouse_VISp_2018-06-14")),FUN=file.remove)

# Drop genes which are not expressed in one of the cell types
unexpressed   = Matrix::rowSums(Allan_Mouse_V1$exp)==0 
Allan_Mouse_V1$exp = Allan_Mouse_V1$exp[!unexpressed,]

# Save
save(Allan_Mouse_V1,file="Output/Allan_Mouse_V1.rda")
```

## Zeisel 2018

The mapping against Zeisel 2015 (and various other datasets is given in Supplementary Table 2): https://www.sciencedirect.com/science/article/pii/S009286741830789X?via%3Dihub#mmc2

https://ars.els-cdn.com/content/image/1-s2.0-S009286741830789X-mmc2.xlsx

The full Zeisel 2018 dataset is huge (19GB). So here we'll map only against the average level of expression in each celltype.

Should you wish to try against the full dataset it is available here: http://mousebrain.org/downloads.html

Another challenge with this dataset is that it contains huge numbers of cell types which are not found in the other datasets.

```{r }
zeisel2018 = load_zeisel_dataset()
mean_exp = convert_zeisel_tibble_to_exp_matrix(tibbleIN=zeisel2018,level=5)

zeisel2018 = list()
zeisel2018$exp = Matrix::Matrix(mean_exp)
zeisel2018$annot = data.frame(cellID=colnames(zeisel2018$exp),celltype=colnames(zeisel2018$exp))
save(zeisel2018,file="Output/zeisel2018.rda")

mapping_zeisel2018_zeisel2015 = read.xlsx("https://ars.els-cdn.com/content/image/1-s2.0-S009286741830789X-mmc2.xlsx")
mapping_zeisel2018_zeisel2015 = mapping_zeisel2018_zeisel2015[grepl("Cortex-",mapping_zeisel2018_zeisel2015$Published.cell.type),]
mapping_zeisel2018_zeisel2015$Published.cell.type = gsub("Cortex-","",mapping_zeisel2018_zeisel2015$Published.cell.type)
mapping_zeisel2018_zeisel2015$`This.study.cell.type.(#.cells)` = gsub("\\(.*","",mapping_zeisel2018_zeisel2015$`This.study.cell.type.(#.cells)`)
mapping_zeisel2018_zeisel2015 = mapping_zeisel2018_zeisel2015 %>% dplyr::rename(zeisel2015=Published.cell.type,zeisel2018=`This.study.cell.type.(#.cells)`) %>%
  dplyr::select(zeisel2015,zeisel2018)
mapping_zeisel2018_zeisel2015$comparison = "zeisel2015 - zeisel2018"
save(mapping_zeisel2018_zeisel2015,file="Output/mapping_zeisel2018_zeisel2015.rda")
```

## Unannotated dataset from Imperial

The next challenge is to map cells from an existing unannotated / unclustered dataset (generated at Imperial) onto the datasets given above. This dataset is human nuclei from similar brain regions to most the datasets given above.

```{r }
load("Data/seurat_PDC91.rda")
imperial_unannotated = list()
imperial_unannotated$exp = seurat_individual@assays$RNA
imperial_unannotated$exp = imperial_unannotated$exp[,seurat_individual$doublets=="Singlet"]
imperial_unannotated$annot = data.frame(cellID=colnames(imperial_unannotated$exp),celltype="Unknown")
save(imperial_unannotated,file="Output/imperial_unannotated.rda")
```

## Reduce datasets to all have the same gene symbols

```{r }
allDatasets = list()
allDatasets[["zeisel2015"]] = zeisel2015
allDatasets[["tasic"]] = tasic
allDatasets[["allanMTG_intron"]] = allanMTG_intron # HUMAN
allDatasets[["allanMTG_exon"]] = allanMTG_exon  # HUMAN
allDatasets[["Allan_Mouse_V1"]] = Allan_Mouse_V1
allDatasets[["zeisel2018"]] = zeisel2018
allDatasets[["imperial_unannotated"]] = imperial_unannotated

# Two of the datasets are from human, so convert HNGNC symbols to MGI
devtools::install_github("nathanskene/One2one")
keepGenesA = One2One::ortholog_data_Mouse_Human$orthologs_one2one$human.symbol %in% rownames(allDatasets[["allanMTG_intron"]]$exp)
keepGenesB = One2One::ortholog_data_Mouse_Human$orthologs_one2one$human.symbol %in% rownames(allDatasets[["allanMTG_exon"]]$exp)
keepGenesC = One2One::ortholog_data_Mouse_Human$orthologs_one2one$human.symbol %in% rownames(allDatasets[["imperial_unannotated"]]$exp)
keepGenes  = keepGenesA & keepGenesB & keepGenesC
keptGenes  = One2One::ortholog_data_Mouse_Human$orthologs_one2one$human.symbol[keepGenes]
o2o = One2One::ortholog_data_Mouse_Human$orthologs_one2one
rownames(o2o) = One2One::ortholog_data_Mouse_Human$orthologs_one2one$human.symbol
o2o = o2o[keptGenes,]
allDatasets[["allanMTG_intron"]]$exp = allDatasets[["allanMTG_intron"]]$exp[keptGenes,]
allDatasets[["allanMTG_exon"]]$exp = allDatasets[["allanMTG_exon"]]$exp[keptGenes,]
allDatasets[["imperial_unannotated"]]$exp = allDatasets[["imperial_unannotated"]]$exp[keptGenes,]
rownames(allDatasets[["allanMTG_intron"]]$exp) = o2o$mouse.symbol
rownames(allDatasets[["allanMTG_exon"]]$exp) = o2o$mouse.symbol
rownames(allDatasets[["imperial_unannotated"]]$exp) = o2o$mouse.symbol
keptGenes_ms = rownames(allDatasets[["allanMTG_exon"]]$exp)

# Find genes used across all the datasets
for(i in 1:length(allDatasets)){ keptGenes_ms = keptGenes_ms[keptGenes_ms %in%rownames(allDatasets[[i]]$exp)] }
for(i in 1:length(allDatasets)){ allDatasets[[i]]$exp = allDatasets[[i]]$exp[keptGenes_ms,] }

save(allDatasets,file="Output/allDatasets.rda")
```

Now merge all the mappings into one file

```{r }
mappings = list()
mappings[["mapping_tasic_zeisel2015"]]     = mapping_tasic_zeisel2015
mappings[["mapping_allanMTG_intron_exon"]] = mapping_allanMTG_intron_exon
mappings[["zeisel2018_zeisel2015"]]        = mapping_zeisel2018_zeisel2015
save(mappings,file="Output/mappings.rda")
```