---
title: "Build counts table"
output: html_document
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, tidy = T, cache = T)
print(getwd())

wd = getwd()
library(DESeq2)
library(ggplot2)
library(pheatmap)
library(RColorBrewer)
library(gplots)
library(ggrepel)
library(reshape2)
library(plyr)
library(dplyr)
library(stringr)
library(cowplot)
library(expss)
library(readr)
library(rhdf5)
library(tximport)

#library(reticulate)
#use_python("/usr/bin/python")
```


# get the raw data

To get the data:

```{bash, eval = F}
#BASH
mkdir -p online_cds2019
rsync -avz E5_comp1_psmn:/Xnfs/lbmcdb/Semon_team/data/Convergenomix/Analyses/ropipe_pipeline2_output/online_cds/kallisto online_cds --exclude=*idx
rsync -avz E5_comp1_psmn:/Xnfs/lbmcdb/Semon_team/data/Convergenomix/Analyses/ropipe_pipeline2_output/online_cds/blastx_eggnog online_cds
```

```{bash, eval = F}
#BASH
mkdir -p data_datalocal_cds2019
rsync -avz E5_comp1_psmn:/Xnfs/lbmcdb/Semon_team/data/Convergenomix/Analyses/ropipe_pipeline2_output/datalocal_cds/kallisto data_datalocal_cds2019 --exclude=*idx
rsync -avz E5_comp1_psmn:/Xnfs/lbmcdb/Semon_team/data/Convergenomix/Analyses/ropipe_pipeline2_output/datalocal_cds/blastx_eggnog data_datalocal_cds2019
```

```{bash, eval = F}
#BASH
mkdir -p extra_files_for_annotation
rsync -avz E5_comp1_psmn:/Xnfs/lbmcdb/Semon_team/data/Convergenomix/Analyses/extra_files_for_annotation/* extra_files_for_annotation
```

# Function to build a counts table

```{r chargement tables}
#R

###Preparation of the clNOG table (reorganize EggNOG assignation to split column 6 and obtain 1:1 cluster:gene)
#Download roNog_Member table
#wget http://eggnog5.embl.de/download/eggnog_5.0/per_tax_level/9989/9989_members.tsv.gz

#Cluster EggNOG assignation to Ensembl gene ID

clNOG = read.csv('extra_files_for_annotation/9989_members.tsv', h=F, sep = "\t") #Change the member version 9989_members.tsv in June 2020. Was before roNOG_member.tsv

#colnames(clNOG) = c("Class","ClusterID","NbOfProteins","NbOfspecies","NoInfo","GeneID") #for older version
colnames(clNOG) = c("dataset","ClusterID","NbOfProteins","NbOfspecies","GeneID","TaxonID")
dim(clNOG)

#To detrmine number of mouse sequences in one cluster (to remove paralogs)
clNOG$NbOfMus = 0
clNOG$NbOfMus = str_count(clNOG$GeneID, "MUS") 
a = table(clNOG$NbOfMus)
clNOGortho <- clNOG[clNOG$NbOfMus == 1,]


#Creation of a new clNOG table with splitting the col 6
v=gsub( ","," ",as.character(clNOGortho$"GeneID")) 
vl=lapply(v,function(x){strsplit(x," ")})
n=lapply(v,function(x){length(strsplit(x," ")[[1]])})
ul=unlist(vl)
g=rep(clNOGortho$ClusterID,n)
clNOGok=data.frame(ClusterID=g,GeneID=ul)
head(clNOGok)

#enleve le numero devant  GeneEnsProtID

cm=gsub( "^[0-9]*.","",as.character(clNOGok$"GeneID"))
cml=lapply(cm,function(x){strsplit(x," ")})
cmn=lapply(cm,function(x){length(strsplit(x," ")[[1]])})
cmul=unlist(cml)
cmg=rep(clNOGok$ClusterID,cmn)
clNOGokMus2=data.frame(ClusterID=cmg,GeneEnsProtID=cmul)
head(clNOGokMus2)

#Use of Ensembl ID converter to get the IDs from v70 to v88 (current one)
#Import Ensembl IDs, lengths and MGI IDs (downloaded from Core Ensembl release 70)

#BiomartMus = read.csv('/extra_files_for_annotation/table_biomart_mus_v99.tsv', h=T, sep = "\t")

BiomartMus = read.csv('extra_files_for_annotation/Biomart_GRCm38.p6.txt', h=F, sep = "\t")

#WARNING Biomart version
colnames(BiomartMus) = c("EnsGeneID","EnsTranscriptID","MGI","EnsTranscriptLen","GeneEnsProtID")

length(grep("ENSMUSP", BiomartMus$GeneEnsProtID)) #Number of prot ID in Biomart 68348


## HERE : add patch to remove cases were 1 MGI != 1 cluster 
correspMGICluster = merge(x=clNOGokMus2, y=BiomartMus, by=c("GeneEnsProtID"), all=F, incomparables = NULL)

MusAssignIDs_short <- unique(subset(correspMGICluster, select = c(ClusterID, MGI)))
MusAssignIDs_short$MGI = as.character(MusAssignIDs_short$MGI)
MusAssignIDs_short$ClusterID = as.character(MusAssignIDs_short$ClusterID ) 

a = table(MusAssignIDs_short$ClusterID)
b = table(MusAssignIDs_short$MGI)

MusAssignIDs_short$unique_mgi = F
MusAssignIDs_short$unique_mgi[ MusAssignIDs_short$ClusterID %in% names(a[a==1])] = T

MusAssignIDs_short$several_cluster = T
MusAssignIDs_short$several_cluster[ MusAssignIDs_short$MGI %in% names(b[b==1])] = F

table(severalClus=MusAssignIDs_short$several_cluster,uniqMGI=MusAssignIDs_short$unique_mgi  )

MusAssignIDs_short_clean=MusAssignIDs_short[MusAssignIDs_short$unique_mgi==T&MusAssignIDs_short$several_cluster==F,]
table(MusAssignIDs_short_clean$several_cluster,MusAssignIDs_short_clean$unique_mgi  )

clNOGok2= merge(x=clNOGok, y=MusAssignIDs_short_clean, by=c("ClusterID"), all=F, incomparables = NULL)

#Define function to get path and run tximport

get_blastx_path <- function(sp_name, chemin_dir_blast){
  #print(paste0(chemin_dir_blast,"/",sp_name,"_*"))
  unlist(Sys.glob(paste0(chemin_dir_blast,sp_name,"/",sp_name,"_blastx.out")))
}

get_kallisto_path <- function(sp_name, chemin_dir_kallisto){
  #print(paste0(chemin_dir_kallisto,"/",sp_name,"/*/abundance.h5"))
  unlist(Sys.glob(paste0(chemin_dir_kallisto,sp_name,"/*/abundance.h5")))
}

# fonction qui avec un blast et le chemin vers les fichier de quantif kallisto lance tximport 
annot <- function(sp_name, chemin_dir_blast , chemin_dir_kallisto,prefix) {
  print(sp_name)
  print("blastx")
  blastx = get_blastx_path(sp_name, chemin_dir_blast)
  print(blastx)

  if (length(blastx) > 0 ){
    #blastx de rongeurs 
    rodblastx = read.csv(blastx , h=F, sep = "\t")
    colnames(rodblastx) = c("ContigID","GeneID","Evalue","BitScore","Length","pident","qstart","qend","sstart","send" , "qlen")
    head(rodblastx) 
    MusAssignIDsInt = merge(x=clNOGok2, y=rodblastx[ , c("ContigID", "GeneID","Evalue","BitScore","Length" ,"qlen", "pident","qstart","qend","sstart","send")], by=c("GeneID"), all=F)
    length(unique(MusAssignIDsInt$GeneID)) 
    ####Pb is here dans blastx GeneID is the one corresponding to clNOGok2
    
    #table qui contient le nom du contig le cluster eggnog associer la longeur du contig assemble et le mgi associer au cluster eggnog
    MusAssignIDsInt$MGI_ClusterID=paste(MusAssignIDsInt$MGI,MusAssignIDsInt$ClusterID, sep="_")
    table_info_long_MGI=MusAssignIDsInt[,c("ClusterID","ContigID","qlen","MGI_ClusterID")]
    #table_info_long_MGI=MusAssignIDsInt[,c("ClusterID","ContigID","qlen","MGI")]

    # lecture dans tximport à partir de mapping kallisto
    
    files_quantif <- get_kallisto_path(sp_name, chemin_dir_kallisto)
    print("kallisto")
    print(files_quantif)
    
    #recupere les nom du srr et de l'espèce 
    files_indiv=sapply(files_quantif,function(x){strsplit(x,"/")[[1]][5]}) ##Warning path level to abundance file
    indiv_name=paste0(prefix,"_",files_indiv,"_",sp_name)
    names(files_quantif)=indiv_name
    
    tx2gene=MusAssignIDsInt[MusAssignIDsInt$pident > 60 & MusAssignIDsInt$Length > 50 , c("ContigID","MGI")]
    
    kallisto.cells=tximport(files_quantif, type = "kallisto",txOut=TRUE)
    kallisto.cells.counts=kallisto.cells$"counts"
    ## modif marie pour corriger pour la taille
    kallisto.cells.tx2gene=tximport(files_quantif, type = "kallisto",tx2gene = tx2gene,countsFromAbundance="lengthScaledTPM")
    kallisto.cells.tx2gene.counts=kallisto.cells.tx2gene$"counts"
    kallisto.cells.tx2gene.abundance=kallisto.cells.tx2gene$"abundance"
    
    
    AssignIDsIntCounts = merge(x=MusAssignIDsInt, y=kallisto.cells.counts, by.x=c("ContigID"),by.y=0, all=F)
    
    
    #retourne les comptes bruts et l'abondance 
    # return(kallisto.cells)
    return(list(counts_per_genes = kallisto.cells.tx2gene.counts,
                abundance_per_genes = kallisto.cells.tx2gene.abundance,
                blast_counts=AssignIDsIntCounts,
                contig_info_len=table_info_long_MGI))
  }
}

```

# Build the counts table for online_cds2019

```{r}
#R

#LAUNCHING

data_dir="data_online_cds2019/"
prefix= "onl_cds"
outdir=gsub("data","counts_table", data_dir)

#Get species list
list_sp_online <- Sys.glob(paste0(data_dir,"kallisto/*"))
list_sp_online <- unlist(lapply(list_sp_online, basename))

#Start annot function
list_online_tximp <- lapply(list_sp_online, annot, 
                            chemin_dir_blast= paste0(data_dir,"blastx_eggnog/"), 
                            chemin_dir_kallisto=  paste0(data_dir,"kallisto/"),
                            prefix = prefix)
names(list_online_tximp) = list_sp_online
list_online_tximp = list_online_tximp[! sapply(list_online_tximp, is.null)]


system(paste0("mkdir -p ", outdir, "/Expression_blast_per_species_per_contig"))
system(paste0("mkdir -p ", outdir, "/Counts_per_species_per_MGI"))
system(paste0("mkdir -p ", outdir, "/Abundances_per_species_per_MGI"))
system(paste0("mkdir -p ", outdir, "/tableau_info_contig_length"))

messages = lapply(list_sp_online[list_sp_online %in%names(list_online_tximp) ],function(x){
  write.table(list_online_tximp[[x]]$counts_per_genes,file=paste0(outdir, "/Counts_per_species_per_MGI/",prefix,"_",x,".tsv"),quote=F,row.names=T, sep="\t");
  write.table(list_online_tximp[[x]]$abundance_per_genes,file=paste0(outdir, "/Abundances_per_species_per_MGI/",prefix,"_",x,".tsv"),quote=F,row.names=T,sep="\t");
  write.table(list_online_tximp[[x]]$blast_counts,file=paste0(outdir,"/Expression_blast_per_species_per_contig/",prefix,"_",x,".tsv"),quote=F,row.names=T, sep="\t");
  write.table(list_online_tximp[[x]]$contig_info_len,file=paste0(outdir,"/tableau_info_contig_length/",prefix,"_",x,".tsv"),quote=F,row.names=T, sep="\t")
  }
)

```


```{python}
# PYTHON (run only when you compile the Rmd with knitr (not working if you manually run the chunk separately))

import pandas as pd
import glob
import os
os.chdir( "counts_table_online_cds2019")
print(os.getcwd())

a=glob.glob("Counts_per_species_per_MGI/*.tsv")

#os.chdir( "/home/crey02/Documents/Projets/Convergences/Online_rodents/Convergent-xeric-rodent_online-dataset/comparaison_DE_2018_2019")

list_table=[]
for i in a:
  df = pd.read_csv(i , sep='\t')
  list_table.append(df)

result=pd.concat(list_table,  axis=1, join='outer' )

result.to_csv("online_cds_counts_fev2022.tsv" ,index=True, sep="\t")

a=glob.glob("Abundances_per_species_per_MGI/*.tsv")

list_table=[]
for i in a:
  list_table.append(pd.read_csv(i , sep='\t'))

result=pd.concat(list_table,  axis=1, join='outer' )

result.to_csv("online_cds_abundances_fev2022.tsv" , sep="\t")

```

# Build the counts table for datalocal_cds2019

```{r}
# R

#LAUNCHING

data_dir="data_local_cds2019/"
prefix= "loc_cds"
outdir=gsub("data_","counts_table_", data_dir)

#Get species list
list_sp_local <- Sys.glob(paste0(data_dir,"kallisto/*"))
list_sp_local <- unlist(lapply(list_sp_local, basename))

#Start annot function
# MARIE : était écrit "assemblage 2020
list_local_tximp <- lapply(list_sp_local, annot, 
                            chemin_dir_blast= paste0(data_dir,"blastx_eggnog/"), 
                            chemin_dir_kallisto=  paste0(data_dir,"kallisto/"),
                            prefix = prefix)
names(list_local_tximp) = list_sp_local
list_local_tximp = list_local_tximp[! sapply(list_local_tximp, is.null)]


system(paste0("mkdir -p ", outdir, "/Expression_blast_per_species_per_contig"))
system(paste0("mkdir -p ", outdir, "/Counts_per_species_per_MGI"))
system(paste0("mkdir -p ", outdir, "/Abundances_per_species_per_MGI"))
system(paste0("mkdir -p ", outdir, "/tableau_info_contig_length"))

messages = lapply(list_sp_local[list_sp_local %in%names(list_local_tximp) ],function(x){
  write.table(list_local_tximp[[x]]$counts_per_genes,file=paste0(outdir, "/Counts_per_species_per_MGI/",prefix,"_",x,".tsv"),quote=F,row.names=T, sep="\t");
  write.table(list_local_tximp[[x]]$abundance_per_genes,file=paste0(outdir, "/Abundances_per_species_per_MGI/",prefix,"_",x,".tsv"),quote=F,row.names=T,sep="\t");
  write.table(list_local_tximp[[x]]$blast_counts,file=paste0(outdir,"/Expression_blast_per_species_per_contig/",prefix,"_",x,".tsv"),quote=F,row.names=T, sep="\t");
  write.table(list_local_tximp[[x]]$contig_info_len,file=paste0(outdir,"/tableau_info_contig_length/",prefix,"_",x,".tsv"),quote=F,row.names=T, sep="\t")
  }
)
```



```{python}

# PYTHON (run only when you compile the Rmd with knitr (not in Rstudio even if you mannualy change the working directory))

import pandas as pd
import glob
import os
os.chdir( "counts_table_local_cds2019/")
print(os.getcwd())

a=glob.glob("Counts_per_species_per_MGI/*.tsv")

list_table=[]
for i in a:
  df = pd.read_csv(i , sep='\t')
  list_table.append(df)

result=pd.concat(list_table,  axis=1, join='outer' )

result.to_csv("datalocal_cds_counts_fev2022.tsv" ,index=True, sep="\t")

a=glob.glob("Abundances_per_species_per_MGI/*.tsv")

list_table=[]
for i in a:
  list_table.append(pd.read_csv(i , sep='\t'))

result=pd.concat(list_table,  axis=1, join='outer' )

result.to_csv("datalocal_cds_abundances_fev2022.tsv" , sep="\t")
```

# Merge the count tables: local and online

```{r, eval = F}

#Count table
online_cds2019_counts <- read.csv('counts_table_online_cds2019/online_cds_counts_fev2022.tsv', h=T, sep = "\t", row.names=1)
local_cds2019_counts <- read.csv('counts_table_local_cds2019/datalocal_cds_counts_fev2022.tsv', h=T, sep = "\t", row.names=1)

total_counts_cds2022 <- merge(online_cds2019_counts,local_cds2019_counts, by = 0, all = T)
write.table(total_counts_cds2022,file="datalocalandonline_cds_counts_2022.tsv",quote=F,row.names=F, sep="\t")

#Abundance table
online_cds2019_abundance <- read.csv('counts_table_online_cds2019/online_cds_abundances_fev2022.tsv', h=T, sep = "\t", row.names=1)
local_cds2019_abundance <- read.csv('counts_table_local_cds2019/datalocal_cds_abundances_fev2022.tsv', h=T, sep = "\t", row.names=1)

total_abundance_cds2019 <- merge(online_cds2019_abundance,local_cds2019_abundance, by = 0, all = T)
write.table(total_abundance_cds2019,file="datalocalandonline_cds_abundance_2022.tsv",quote=F,row.names=F, sep="\t")

```
