<jupyter_start><jupyter_text>find all way hc one to one ortholog<jupyter_code>orthopath="/home/rcf-40/amalthom/panfs/software/genomes/orthology_ensembl94/one2one_hc"
species = c("cjacchus","hsapiens","mmulatta","mmusculus","ocuniculus","rnorvegicus")
species.full = c("Callithrix_jacchus","Homo_sapiens","Macaca_mulatta",
                 "Mus_musculus","Oryctolagus_cuniculus","Rattus_norvegicus")

l = list()
ensembl_unique = vector()
for(s in species){
   if(s == "hsapiens") next
    print(s)
 df = read.table(paste0(orthopath,"/hsapiens-vs-",s,"_one2one_hc.tsv"),
                       sep="\t",header=T,stringsAsFactors = F)
head(df)
l[[s]] = df$ensembl_gene_id
ensembl_unique = c(ensembl_unique,as.vector(df$ensembl_gene_id))    
}
ensembl_unique = unique(ensembl_unique)
sapply(l,length)
sum(sapply(l,length))
length(ensembl_unique)
overlap = as.data.frame(sapply(l,function(x) as.integer(ensembl_unique %in% unlist(x))))
row.names(overlap) = ensembl_unique                               
head(overlap)
commongenes = ensembl_unique[rowSums(overlap) == ncol(overlap)]
length(commongenes)
distribution = vennCounts(as.data.frame(overlap))
distribution
paste(species,collapse="-")
outputdir="carelli2018_analysis"
write.table(as.data.frame(commongenes),
            file=paste0(outputdir,"/common_orthologousgenes_human-",paste(species,collapse="-"),".txt"),
           quote=F,row.names=F)
outputdir="carelli2018_analysis"<jupyter_output><empty_output><jupyter_text>subset countstable based on allway orthologs<jupyter_code>counts_path = "/home/rcf-40/amalthom/staging_work/rnaseq/carelli2018/countstables/"
counts_list = list()
for(i in 1:length(species.full)){
    print(species.full[i])
    if(species.full[i] == "Homo_sapiens"){
       
       df = read.table(paste0(counts_path,"RNAseq_carelli2018_",species.full[i],".tsv"),sep="\t",header=T,row.names=1)
       df = df[commongenes,]
       write.table(df,paste0(outputdir,"/counts_",species.full[i],".tsv"),sep="\t",quote=F,col.names=NA)
        print(dim(df))
        counts_list[[species.full[i]]] = df
    }else{
         df = read.table(paste0(orthopath,"/hsapiens-vs-",species[i],"_one2one_hc.tsv"),
                       sep="\t",header=T,stringsAsFactors = F)
        df = df[match(commongenes,df$ensembl_gene_id),]
        
        t = read.table(paste0(counts_path,"RNAseq_carelli2018_",species.full[i],".tsv"),sep="\t",header=T,row.names=1)
        t = t[df[,3],]
        print(dim(t))
        row.names(t) = commongenes
        write.table(t,paste0(outputdir,"/counts_",species.full[i],".tsv"),sep="\t",quote=F,col.names=NA)
        counts_list[[species.full[i]]] = t

    }
}

combined_df =  do.call(cbind,counts_list)
dim(combined_df)
colnames(combined_df)
write.table(combined_df,file = paste0(outputdir,"/allwayortholgous_counts.tsv"),sep="\t",quote = F,col.names = NA)<jupyter_output><empty_output><jupyter_text>manipulating columns of combined df<jupyter_code>colnames(combined_df)
?tstrsplit
tstrsplit(colnames(combined_df),"\\.")[[2]]
combined_df[1:3,1:3]
colnames(combined_df) = tstrsplit(colnames(combined_df),"\\.")[[2]]
combined_df[1:3,1:3]<jupyter_output><empty_output><jupyter_text>change colnames<jupyter_code>cinfo = "/panfs/cmb-panasas2/amalthom/work_staging/species_analysis/carelli2018_analysis/carelli2018_srpinfo.tsv"
cnames = read.table(cinfo, sep="\t",header=T)
head(cnames)
newcolnames_df = cnames[match(colnames(combined_df),cnames$srr),]
summary(newcolnames_df$srr==colnames(combined_df))
head(newcolnames_df)
species = tstrsplit(newcolnames_df$srr_title,"_")[[2]]
head(species)
organ = tstrsplit(newcolnames_df$srr_title,"_")[[3]]
head(organ)
newcolnames_df$species = species
newcolnames_df$organ = organ
head(newcolnames_df)
write.table(newcolnames_df,file=paste0(outputdir,"/carelli2018_srpinfo_withmeta.tsv"),
            sep="\t",row.names = F,quote=F)
colnames(combined_df) = newcolnames_df$srr_title
combined_df[1:3,1:3]
colnames(combined_df)<jupyter_output><empty_output><jupyter_text>CCA analysis using CCA package<jupyter_code>suppressPackageStartupMessages({
    library(CCA)
    library(ggplot2)
})<jupyter_output><empty_output><jupyter_text>CCA between human & mouse<jupyter_code>colnames(combined_df)[grepl("human",colnames(combined_df))|
                       grepl("mouse",colnames(combined_df))]
dim(combined_df)
s1 = "human"
s2 = "mouse"
df_s1 = combined_df[,grepl(s1,colnames(combined_df))]
df_s2 = combined_df[,grepl(s2,colnames(combined_df))]
dim(df_s1)
dim(df_s2)
cc1 = cc(df_s1,df_s2)
cc1$cor
cc1$xcoef
res.cc = cc1
plot(cc1$cor,type="b")
plt.cc(res.cc)
plt.cc(res.cc,type="v",var.label=T)
cc2 <- comput(df_s1, df_s2, cc1)
# display canonical loadings
cc2[3:6]
names(cc2)
head(cc2$xscores)
names(cc1)
cc1$cor
dim(cc1$xcoef)
cc1$xcoef
names(cc1$scores)
library(geigen)
dim(df_s1)
dim(df_s2)
df_s1 = df_s1[complete.cases(df_s1),]
df_s2 = df_s2[complete.cases(df_s2),]
dim(df_s1)
dim(df_s2)
common = intersect(row.names(df_s1),row.names(df_s2))
df_s1 = df_s1[common,]
df_s2 = df_s2[common,]
dim(df_s1)
dim(df_s2)
t(df_s1)[1:3,1:3]
cc1 = cc(t(df_s1),t(df_s2))
plt.cc(cc1)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16595080046856697, 'fraction_numerical': 0.024990238188207732, 'mean_word_length': 6.350071736011477, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '558667', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 2226, 'n_tokens_neox': 2122, 'n_words': 211, 'path': 'notebooks/CCA_analysis.ipynb', 'repo_name': 'saketkc/re-ribo-smk', 'size': '5122'}