---
title: "Profile_Matrix_Subset_Regions"
author: "Thomas Goralski"
date: '2022-06-21'
output: html_document
---

```{r}
library(rhdf5)
library(Matrix)
library(data.table)
```




```{r}
h5ls("expression_matrix.hdf5")
```


```{r}
h5closeAll()
```


```{r}



h5read("expression_matrix.hdf5", "data/samples")

```


```{r}
con<- file("Allen_Human_Cortex_Expression_Exon.csv", "r")
A<-read.csv(con, nrows=2)



con<- file("Allen_Human_Cortex_Expression_Exon.csv", "r")
A<-read.csv(con, nrows=2)
A[,1:10]

B<-fread("Allen_Human_Cortex_Expression_Exon.csv", nrows=100)
```



```{r}
#use this to pull only the columns you want from metadata. specify class of column name for column you want, enter null for those you dont in order left to right of matrix
#A<-read.csv(con,  colClasses = c(rep("character", 1), rep("NULL", 4), rep("character", 1), rep("NULL", 6), rep("character", 1), rep("NULL", 44)), header = TRUE)

library(rhdf5)

A<-read.csv("Allen_Human_Cortex_Metadata.csv")


A<-A[,c("exp_component_name", "subclass_label", "region_label")]

unique(A$region_label)
subset<-which(A$region_label== "CgG")

A<-A[subset,]

B<-fread("Allen_Human_Cortex_Expression_Exon.csv")

genes<-B$V1

matches<-A$exp_component_name

subset_index<-match(matches, colnames(B))

B<-B[,..subset_index]

rownames(B)<-genes


type_sub<-which(A$subclass_label=="")

A<-A[-type_sub,]
B<-B[,-..type_sub]

#get average expression
Human_Profile <- create_profile_matrix(mtx = B,            
                                    cellAnnots = A,  
                                    cellTypeCol = "subclass_label" , 
                                    cellNameCol = "exp_component_name",           
                                    matrixName = "Allen_Human_Cortex_Metadata", 
                                    outDir = "/varidata/research/home/thomas.goralski",     
                                    geneList= genes,
                                    normalize = FALSE,                
                                    minCellNum = 3,                   
                                    minGenes = 50,                    
                                    scalingFactor = 1,                
                                    discardCellTypes = TRUE)

genesubs<-as.numeric(rownames(Human_Profile))

genes<-genes[genesubs]

rownames(Human_Profile)<-genes
write.csv(Human_Profile, "Human_Profile.csv")


#if file too large the perform below
#split the data in half

dim(B)
dim(A)
B1<-as.matrix(B[,1:61864])
A1<-A[1:61864,]
B2<-as.matrix(B[,61865:123727])   
A2<-A[61865:123727,]

Mop_Profile1 <- create_profile_matrix(mtx = B1,            
                                    cellAnnots = A1,  
                                    cellTypeCol =  "subclass_label", 
                                    cellNameCol = "sample_name" ,           
                                    matrixName = "MOp_Profile1", 
                                    outDir = "/varidata/research/home/thomas.goralski",     
                                    geneList= genes,
                                    normalize = FALSE,                
                                    minCellNum = 15,                   
                                    minGenes = 100,                    
                                    scalingFactor = 1,                
                                    discardCellTypes = TRUE)




Mop_Profile2 <- create_profile_matrix(mtx = B2,            
                                    cellAnnots = A2,  
                                    cellTypeCol =  "subclass_label", 
                                    cellNameCol = "sample_name",           
                                    matrixName = "MOp_Profile2", 
                                    outDir = "/varidata/research/home/thomas.goralski",     
                                    geneList= genes,
                                    normalize = FALSE,                
                                    minCellNum = 15,                   
                                    minGenes = 100,                    
                                    scalingFactor = 1,                
                                    discardCellTypes = TRUE)





dim(Mop_Profile1)
dim(Mop_Profile2)



subset_index_genes_1<-match(rownames(Mop_Profile2), rownames(Mop_Profile1))
subset_index_genes_2<-match(rownames(Mop_Profile1), rownames(Mop_Profile2))

all_row_to_add_1<-c()
for(i in 1:length(subset_index_genes_1)){row<-subset_index_genes_1[i]
if(is.na(row)==TRUE){
row_to_add<-i
all_row_to_add_1<-c(all_row_to_add_1, row_to_add)}}



Mop_Profile2_row_add<-Mop_Profile2[all_row_to_add_1, ]  ####
Mop_Profile2<-Mop_Profile2[-all_row_to_add_1,]







all_row_to_add_2<-c()
for(i in 1:length(subset_index_genes_2)){row<-subset_index_genes_2[i]
  if(is.na(row)==TRUE){
    row_to_add<-i
    all_row_to_add_2<-c(all_row_to_add_2,row_to_add)
  }}

Mop_Profile1_row_add<-Mop_Profile1[all_row_to_add_2, ] ####
Mop_Profile1<-Mop_Profile1[-all_row_to_add_2,]


reorder_matrix_index<-match(rownames(Mop_Profile1),rownames(Mop_Profile2))


Mop_Profile2<-Mop_Profile2[reorder_matrix_index,]





unique_col<-which(colnames(Mop_Profile1)%in%colnames(Mop_Profile2)==FALSE)
Mop_Profile_col_add<-Mop_Profile1[,unique_col] ####
colname_add<-colnames(Mop_Profile1)[unique_col]




col_reorder_index<-match(colnames(Mop_Profile2), colnames(Mop_Profile1))
Mop_Profile1<-Mop_Profile1[,col_reorder_index]
Mop_Profile1_row_add<-Mop_Profile1_row_add[,c(col_reorder_index,unique_col)]



#Mop_profile1 , Mop_profile2
finished <- matrix(nrow = length(rownames(Mop_Profile1)), ncol = length(colnames(Mop_Profile1)))
for (i in 1:length(rownames(Mop_Profile1))) {
  row <- c()
  for (j in 1:length(colnames(Mop_Profile1))) {
    point1 <- Mop_Profile1[i,j]
    point2 <- Mop_Profile2[i,j]
    
    added <- point1 + point2
    ave <- added/2
    row <- c(row,ave)
  }
  finished[i,] <- row
}

rownames_finished<-rownames(Mop_Profile2)
colnames_finished<-colnames(Mop_Profile2)




col_name_final<-colnames(Mop_Profile1_row_add)



finished2<- rbind(finished, Mop_Profile2_row_add)
add_zeros<-rep(0, length(rownames(Mop_Profile2_row_add)))
Mop_Profile_col_add<-c(Mop_Profile_col_add, add_zeros)
finished2<-cbind(finished2, Mop_Profile_col_add)
finished2<- rbind(finished2, Mop_Profile1_row_add)

write.csv(finished2, "Mop_Profile_Final.csv")



```




```{r}


MOs_subset<-which(A$region_label== "MOp")

A<-A[MOs_subset,]

```



```{r}
B<-h5read("expression_matrix.hdf5", "data/counts", index= list(subset_index, NULL ))
genes<-h5read("expression_matrix.hdf5", "data/gene")
samples<- h5read("expression_matrix.hdf5", "data/samples", index= list(MOs_subset))
B<-t(B)
B<-as.matrix(B)

rownames(B)<- genes
colnames(B)<-samples
fwrite(B,"B_total.csv")
```


```{r}
rhdf5
samples<-A$sample_name

```



```{r}
library(SpatialDecon)
```


```{r}
#had to run on HPC...Big ass files
Mos_Profile <- create_profile_matrix(mtx = B,            # cell x gene count matrix
                                    cellAnnots = A,  # cell annotations with cell type and cell name as columns 
                                    cellTypeCol =  "subclass_label",  # column containing cell type
                                    cellNameCol = "sample_name",           # column containing cell ID/name
                                    matrixName = "MOp_Profile", # name of final profile matrix
                                    outDir = "/varidata/research/home/thomas.goralski",       # path to desired output directory, set to NULL if matrix should not be written
                                    geneList= genes,
                                    normalize = FALSE,                # Should data be normalized? 
                                    minCellNum = 100,                   # minimum number of cells of one type needed to create profile, exclusive
                                    minGenes = 100,                    # minimum number of genes expressed in a cell, exclusive
                                    scalingFactor = 1,                # what should all values be multiplied by for final matrix
                                    discardCellTypes = TRUE)          # should cell types be filtered for types like mitotic, doublet, low quality, unknown, etc.
Mop_Profile2 <- create_profile_matrix(mtx = B2,            
                                    cellAnnots = A2,  
                                    cellTypeCol =  "subclass_label", 
                                    cellNameCol = "sample_name",           
                                    matrixName = "MOp_Profile2", 
                                    outDir = "/varidata/research/home/thomas.goralski",     
                                    geneList= genes,
                                    normalize = FALSE,                
                                    minCellNum = 15,                   
                                    minGenes = 100,                    
                                    scalingFactor = 1,                
                                    discardCellTypes = TRUE)



#no comments




```



```{r}
rownames(MOs_Profile_profileMatrix)<-genes
```































```{r}

label<-unique(A$subclass_label)

```












```{r}
#custom way of doing it
uni<- c()


for (i in label){
  if (i %in% unique(A$subclass_label)){
sub<-which(A$subclass_label==i)
subs<- c(i,sub)
uni<- cbind(uni,subs)
}
}
names<-label
colnames(uni)<-names

uni<-uni[-1,]




```



```{r}
profile_matrix_Mos<- c()
for (i in 1:length(colnames(uni))){
  e<-uni[,i]
vec1<-B[,e]
  vec2<-c(i,vec1)
  profile_matrix_Mos<- cbind(profile_matrix_MOs, vec2)
}
```



```{r}
library(SpatialDecon)
profile_matrix <- download_profile_matrix(species = "Mouse",
                                       age_group = "Adult", 
                                       matrixname = "Brain_AllenBrainAtlas")
```

```{r}
sample_subs<-c()
for(i in 1:length(samnam)) {
  thing1<-samnam[i]
  subs<-which(samples==thing1)
  sample_subs<-c(sample_subs, subs)
   for (ii in seq) {
  if(ii==i) {
    cat(paste0(round(ii / length(samnam) * 100), '% completed'))
 
  if (ii == length(samnam)){ cat(': Done')
  }}}}


for(i in 1:length(samples)){
  thing1<-samples[i]
  subs<-which(types==thing1)
  type_subs<-c(type_subs, subs)

   }

sample_subs2<-c()
for(i in length(samnam)){
  thing1<-samnam[i]
  
  for(j in length(samples)){
    thing2<-samples[j]
   
    if (thing1==thing2){
      sample_subs2<-c(sample_subs2,j)
    }
  }
}

```


```{r}
#Timer for for loop
seq<-seq(1,400000, by=1000)
for (ii in seq){
  if(ii==i){
    cat(paste0(round(ii / length(samnam) * 100), '% completed'))
 
  if (ii == length(samnam)){ cat(': Done')
  else cat('\014')
    }
  }
}


for (ii in 1:length(samnam)) {
  for(iii in seq){
    if(ii==iii){
  cat(paste0(round(ii / length(samnam) * 100), '% completed'))
 
  if (ii == length(samnam)) cat(': Done')
  else cat('\014')
    }
  }
} 
```




```{r}
sample_subs<-rep(NA, length(samnam))

for(i in 1:length(samnam)) {
  thing1<-samnam[i]
  subs<-which(samples==thing1)
  sample_subs[i] = subs
   for (ii in seq) {
  if(ii==i) {
    cat(paste0(round(ii / length(samnam) * 100), '% completed'))
 
  if (ii == length(samnam)){ cat(': Done')
  }}}}




row_subs<-rep(NA, length(row2))
adds<-c()
for(i in 1:length(row1)) {
  thing1<-row1[i]
  subs<-which(row2==thing1)
  if (length(subs) != 0){ row_subs[i] <- subs }else {adds<-c(adds, i)}}
   
```





```{r}
which(rownames(Mos_Profile)=="Etv6")


#Mop_profile1 , Mop_profile2
finished <- matrix(nrow = length(rownames(Mop_Profile1)), ncol = length(colnames(Mop_Profile1)))
for (i in 1:length(rownames(Mop_Profile1))) {
  row <- c()
  for (j in 1:length(colnames(Mop_Profile1))) {
    point1 <- Mop_Profile1[i,j]
    point2 <- Mop_Profile2[i,j]
    
    added <- point1 + point2
    ave <- added/2
    row <- c(row,ave)
  }
  finished[i,] <- row
}


```






Begin Here for more efficient matching algorithm
```{r}

#allocate memory for vector
sample_subs<-rep(NA, length(samnam))

for(i in 1:length(samnam)) {
  index<-samnam[i]
  subs<-which(samples==index)
  sample_subs[i] = subs
   for (ii in seq) {
  if(ii==i) {
    cat(paste0(round(ii / length(samnam) * 100), '% completed'))
 
  if (ii == length(samnam)){ cat(': Done')
  }}}}






#allocate memory for vector
sample_subs<-rep(NA, length(samnam))

for(i in 1:length(samnam)) {
  index<-samnam[i]
for(j in 1:length(sample)){
  while 
}}


samnam<- c("A","B","C","D","F")
sample<- c("B", "C", "A", "D","F")
a<-match(samnam,sample)


samnam<- c(1:150000000)
sample<- sample(samnam, 150000000,replace=FALSE)

#using a sorting function
sample_index<-as.character(1:length(sample))
sample_index
sample_factor<-factor(x=sample_index, levels =sample )
levels(sample_factor)
sample_factor

#use a sorting alorithm with indices stored...can do in c++
```

just use match


