---
title: "perform_sexdiff.Rmd"
author: "Sahin Naqvi"
date: "9/15/2018"
output: html_document
editor_options: 
  chunk_output_type: console
---

read in ortholog lists, gene counts, and metadata
```{r}
fivespec_orth = read.table("../intermediates/one2oneorth_emblids.txt",header=TRUE,stringsAsFactors=FALSE,sep="\t",fill=TRUE)

human_counts = as.matrix(read.table(gzfile("../exprvals/gtex.filt.salmon.tximport.adj.counts.txt.gz"),header = TRUE,row.names = 1,sep="\t"))
human_tpm = as.matrix(read.table(gzfile("../exprvals/gtex.filt.salmon.tximport.unadj.tpm.txt.gz"),header = TRUE,row.names = 1,sep="\t"))

cyno_counts = as.matrix(read.table(gzfile("../exprvals/cyno.salmon.tximport.counts.txt.gz"),header = TRUE,row.names = 1,sep="\t"))
cyno_tpm = as.matrix(read.table(gzfile("../exprvals/cyno.salmon.tximport.counts.txt.gz"),header = TRUE,row.names = 1,sep="\t"))

mouse_counts = as.matrix(read.table(gzfile("../exprvals/mouse.salmon.tximport.counts.txt.gz"),header = TRUE,row.names = 1,sep="\t"))
mouse_tpm = as.matrix(read.table(gzfile("../exprvals/mouse.salmon.tximport.counts.txt.gz"),header = TRUE,row.names = 1,sep="\t"))

rat_counts = as.matrix(read.table(gzfile("../exprvals/rat.salmon.tximport.counts.txt.gz"),header = TRUE,row.names = 1,sep="\t"))
rat_tpm = as.matrix(read.table(gzfile("../exprvals/rat.salmon.tximport.counts.txt.gz"),header = TRUE,row.names = 1,sep="\t"))

dog_counts = as.matrix(read.table(gzfile("../exprvals/dog.salmon.tximport.counts.txt.gz"),header = TRUE,row.names = 1,sep="\t"))
dog_tpm = as.matrix(read.table(gzfile("../exprvals/dog.salmon.tximport.counts.txt.gz"),header = TRUE,row.names = 1,sep="\t"))

all_counts = cbind(human_counts[subset(fivespec_orth,Gene.name %in% rownames(human_counts))$Gene.name,],cyno_counts[subset(fivespec_orth,Gene.name %in% rownames(human_counts))$Crab.eating.macaque.gene.stable.ID,],mouse_counts[subset(fivespec_orth,Gene.name %in% rownames(human_counts))$Mouse.gene.stable.ID,],rat_counts[subset(fivespec_orth,Gene.name %in% rownames(human_counts))$Rat.gene.stable.ID,],dog_counts[subset(fivespec_orth,Gene.name %in% rownames(human_counts))$Dog.gene.stable.ID,])
all_tpm = cbind(human_tpm[subset(fivespec_orth,Gene.name %in% rownames(human_tpm))$Gene.name,],cyno_tpm[subset(fivespec_orth,Gene.name %in% rownames(human_tpm))$Crab.eating.macaque.gene.stable.ID,],mouse_tpm[subset(fivespec_orth,Gene.name %in% rownames(human_tpm))$Mouse.gene.stable.ID,],rat_tpm[subset(fivespec_orth,Gene.name %in% rownames(human_tpm))$Rat.gene.stable.ID,],dog_tpm[subset(fivespec_orth,Gene.name %in% rownames(human_tpm))$Dog.gene.stable.ID,])


human_metadata = read.table("../metadata/human.metadata.txt",header = TRUE,stringsAsFactors = FALSE)
nonhuman_metadata = read.table("../metadata/sample.metadata.txt",header = TRUE,stringsAsFactors = FALSE,sep="\t")
nonhuman_metadata = subset(nonhuman_metadata,Tissue != "Testis")
nonhuman_metadata$comb_label = paste(nonhuman_metadata$Species,nonhuman_metadata$Sex,nonhuman_metadata$Tissue,nonhuman_metadata$label,sep="_")

all_metadata = data.frame(rbind(as.matrix(human_metadata[,c("RUN","Species","Tissue","Sex")]),as.matrix(nonhuman_metadata[,c("comb_label","Species","Tissue","Sex")])))
colnames(all_metadata) = c("label","Species","Tissue","Sex")
all_metadata$label = as.character(all_metadata$label)
all_metadata$Species = as.character(all_metadata$Species)
all_metadata$Tissue = as.character(all_metadata$Tissue)
all_metadata$Sex = as.character(all_metadata$Sex)


# calculate median tpms in each tissue-species-sex combination, to be used later for filtering prior to sex DE
require(edgeR)
species_mf_med = list()
for (tissue in unique(all_metadata$Tissue)){
  spec5_tissue_tpm_temp = all_tpm[,subset(all_metadata,(Tissue==tissue))$label]
  tmm = DGEList(counts=spec5_tissue_tpm_temp)
  tmm = calcNormFactors(tmm,doWeighting = FALSE)
  tmm_cpm = cpm(tmm)

  print(tissue)
  species_mf_med[[tissue]] = list()
  species_mf_med[[tissue]][["Male"]] = matrix(nrow=nrow(all_counts),ncol=length(unique(all_metadata$Species)),dimnames=list(rownames(all_counts),unique(all_metadata$Species)))
  species_mf_med[[tissue]][["Female"]] = matrix(nrow=nrow(all_counts),ncol=length(unique(all_metadata$Species)),dimnames=list(rownames(all_counts),unique(all_metadata$Species)))
  for (species in unique(all_metadata$Species)){
    print(species)
    species_mf_med[[tissue]][["Male"]][,species] = apply(tmm_cpm[,subset(all_metadata,(Species==species)&(Tissue==tissue)&(Sex=="Male"))$label,drop=FALSE],1,median)
    species_mf_med[[tissue]][["Female"]][,species] = apply(tmm_cpm[,subset(all_metadata,(Species==species)&(Tissue==tissue)&(Sex=="Female"))$label,drop=FALSE],1,median)
  }
}

```


perform DE for conserved sex bias as well as DE for each tissue-species pair 
```{r}
require(limma)
source('../scripts/getSexBiasStats.R')

toptab_specsep_list = list()
toptab_spec5_list = list()

for (tissue in unique(all_metadata$Tissue)) {
  print(tissue)
  weightopt =  'y'
  pcaopt = 'n'

  spec5_meta = subset(all_metadata,(Tissue==tissue))
  rownames(spec5_meta) = spec5_meta$label
  spec5_tissue_tpm_temp = all_tpm[,rownames(spec5_meta)]
  spec5_tissue_counts_temp = all_counts[,rownames(spec5_meta)]
  keep = (rowSums(species_mf_med[[tissue]][["Male"]]>1) > 3) | (rowSums(species_mf_med[[tissue]][["Female"]]>1) > 3)
  if (tissue =="Thyroid") {keep = keep[-5825]}
  #retlist = getSexBiasStats(spec5_tissue_counts_temp[keep,],spec5_meta,'voom_random',outlierWeight='y')
  #toptab_spec5_list[[tissue]] = retlist[[1]]

  retlist = getSexBiasStats(spec5_tissue_counts_temp[keep,],spec5_meta,'fisher_voom',outlierWeight='n')
  toptab_specsep_list[[tissue]] = retlist[[1]]

}


spec5_pval_mat = matrix(nrow=nrow(all_counts),ncol=length(unique(all_metadata$Tissue)),dimnames = list(rownames(all_counts),unique(all_metadata$Tissue)))
spec5_coef_mat = spec5_pval_mat
for (tissue in unique(all_metadata$Tissue)){
  spec5_pval_mat[rownames(toptab_spec5_list[[tissue]]),tissue] = toptab_spec5_list[[tissue]][,"pval"]
  spec5_coef_mat[rownames(toptab_spec5_list[[tissue]]),tissue] = toptab_spec5_list[[tissue]][,"coef"]
}
spec5_pval_mat = spec5_pval_mat[which(apply(!is.na(spec5_pval_mat),1,sum) > 0),]
spec5_coef_mat = spec5_coef_mat[which(apply(!is.na(spec5_coef_mat),1,sum) > 0),]


```


Apply additional filters on log2(M/F) for conserved sex bias. Produces gene x tissue matrix of conserved sex bias (-1 female, 1 male) 
```{r}

# Take all gene-tissue pairs with qvalue < 0.1
require(qvalue)
q = qvalue(c(spec5_pval_mat))
spec5_qval_mat = matrix(q$qvalues,nrow=nrow(spec5_pval_mat),ncol=ncol(spec5_pval_mat),dimnames = list(rownames(spec5_pval_mat),colnames(spec5_pval_mat)))

bias_matrix_tissep_spec5 = matrix(0,nrow=nrow(spec5_coef_mat),ncol = ncol(spec5_coef_mat),dimnames=list(rownames(spec5_coef_mat),colnames(spec5_coef_mat)))

for (tissue in unique(all_metadata$Tissue)){
  fdr = 0.1
  adj = spec5_qval_mat[,tissue]
  bias_matrix_tissep_spec5[names(adj[which(adj < fdr)]),tissue] = 1
}
bias_matrix_tissep_spec5 = bias_matrix_tissep_spec5 * (spec5_coef_mat/abs(spec5_coef_mat))



# require same direction of bias > 1.05-fold in 4/5 species
lfc_cut = log2(1)
bias_matrix_tissep_spec5_nospecint = bias_matrix_tissep_spec5

for (tissue in colnames(bias_matrix_tissep_spec5)){
  common_rownames = intersect(rownames(spec5_coef_mat),rownames(toptab_specsep_list[[tissue]][["coef"]]))
  diff_rownames = unique(c(setdiff(rownames(toptab_spec5_list[[tissue]]),rownames(toptab_specsep_list[[tissue]][["coef"]])),setdiff(rownames(toptab_specsep_list[[tissue]][["coef"]]),rownames(toptab_spec5_list[[tissue]]))))
  tissue_specdiffsign = common_rownames[which((apply(toptab_specsep_list[[tissue]][["coef"]][common_rownames,]*spec5_coef_mat[common_rownames,tissue]/abs(spec5_coef_mat[common_rownames,tissue]) < log2(1.05),1,sum) %in% c(2,3,4,5)) | (toptab_specsep_list[[tissue]][["coef"]][common_rownames,"Human"]*spec5_coef_mat[common_rownames,tissue] < 0) | (abs(toptab_specsep_list[[tissue]][["coef"]][common_rownames,"Human"]) < log2(1.05)) | (apply(toptab_specsep_list[[tissue]][["pval"]][common_rownames,1:5] > 1,1,sum) > 1))]
  tissue_specdiffsign = unique(c(tissue_specdiffsign,diff_rownames))
  bias_matrix_tissep_spec5_nospecint[intersect(rownames(bias_matrix_tissep_spec5_nospecint),tissue_specdiffsign),tissue] = 0
}

bias_matrix_tissep_spec5_nospecint[is.na(bias_matrix_tissep_spec5_nospecint)] = 0
bias_matrix_tissep_spec5_nospecint = bias_matrix_tissep_spec5_nospecint[which(apply(abs(bias_matrix_tissep_spec5_nospecint),1,sum) > 0),]

bias_matrix_tissep_spec5_nospecint_coef = bias_matrix_tissep_spec5_nospecint*abs(spec5_coef_mat[rownames(bias_matrix_tissep_spec5_nospecint),colnames(bias_matrix_tissep_spec5_nospecint)])
bias_matrix_tissep_spec5_nospecint_coef[which(is.na(bias_matrix_tissep_spec5_nospecint_coef))] = 0

```


use mashr to identify sex biases in each tissue-species combination - WARNING: this takes a long time (several hours), so 
```{r}
require(mashr)

missing_betaval = 0
missing_sebetaval = 1000

# note: alternatively, these two matrices can also be read in from the provided supplementary .txt files
mashr_beta_spectis = matrix(missing_betaval,nrow=nrow(spec5_coef_mat),ncol=60)
rownames(mashr_beta_spectis) = rownames(spec5_coef_mat)
mashr_se_spectis = mashr_beta_spectis + missing_sebetaval

# fill out effect size and SE matrices
i=0
spectis_names = c()
for (tissue in colnames(bias_matrix_tissep_spec5)){
  for (spec in colnames(toptab_specsep_list[[tissue]][["coef"]])){
    i = i + 1
    colname = paste(tissue,spec,sep = "_")
    spectis_names = c(spectis_names,colname)
    mashr_beta_spectis[intersect(rownames(toptab_specsep_list[[tissue]][["coef"]]),rownames(mashr_beta_spectis)),i] = toptab_specsep_list[[tissue]][["coef"]][intersect(rownames(toptab_specsep_list[[tissue]][["coef"]]),rownames(mashr_beta_spectis)),spec]
    mashr_se_spectis[intersect(rownames(toptab_specsep_list[[tissue]][["coef"]]),rownames(mashr_beta_spectis)),i] = toptab_specsep_list[[tissue]][["coef_se"]][intersect(rownames(toptab_specsep_list[[tissue]][["coef"]]),rownames(mashr_beta_spectis)),spec]
  }
}
colnames(mashr_beta_spectis) = spectis_names
colnames(mashr_se_spectis) = spectis_names

#read in sparse factors
sf = as.matrix(read.delim("../intermediates/salmon.starref.tximport.voom.spec5orth.sfa_F.out",header = FALSE))
U.sf = mashr:::cov_from_factors(sf,name = c("SFA_1","SFA_2","SFA_3","SFA_4","SFA_5"))

# run mashr
dat   = set_mash_data(mashr_beta_spectis, mashr_se_spectis)
m.1by1 = mash_1by1(dat)
strong = get_significant_results(m.1by1,0.05)
U.c = cov_canonical(dat)
U.pca = cov_pca(dat,5,strong)
U.ed = cov_ed(dat, c(U.pca,U.sf), strong)

m.c.ed.sf.spectis = mash(dat, c(U.c,U.ed))

m.lsfr.spectis = get_lfsr(m.c.ed.sf.spectis)
m.pm.spectis = get_pm(m.c.ed.sf.spectis)

#alternatively, these two matrices can also be read in from the provided supplementary .txt files
m.lsfr.spectis = as.matrix(read.delim("../intermediates/one2oneorth_60spectis_mashr_lfsr.txt",header = TRUE,stringsAsFactors = FALSE,row.names = 1))
m.pm.spectis = as.matrix(read.delim("../intermediates/one2oneorth_60spectis_mashr_pm.txt",header = TRUE,stringsAsFactors = FALSE,row.names = 1))
```


infer lineage of sex bias via parsimony. Produces gene x tissue matrices of sex bias (-1 female, 1 male) for each type of lineage-specific bias considered
```{r}
lfsr_forcut = 0.05
lfc_cut = log2(1.05)
m.lsfr.spectis.signed = (1*(m.lsfr.spectis < lfsr_forcut) * m.pm.spectis) / abs(m.pm.spectis)
m.lsfr.spectis.signed.lfccut = m.lsfr.spectis.signed*(1*(abs(m.pm.spectis) > lfc_cut))


human_gain_mat = matrix(0,nrow=nrow(m.lsfr.spectis.signed),ncol=length(colnames(bias_matrix_tissep_spec5)),dimnames = list(rownames(m.lsfr.spectis.signed),colnames(bias_matrix_tissep_spec5)))
cyno_gain_mat = human_gain_mat
mouse_gain_mat = human_gain_mat
rat_gain_mat = human_gain_mat
dog_gain_mat = human_gain_mat
human_switch_mat = human_gain_mat
primate_gain_mat = human_gain_mat
primate_loss_mat = human_gain_mat
rodent_gain_mat = human_gain_mat
boroeuth_mat = human_gain_mat
single_loss_mat = human_gain_mat
rodent_loss_mat = human_gain_mat
complex_gain_mat = human_gain_mat # representing multiple gains
complex_anc_mat = human_gain_mat # representing multiple losses
complex_mat = human_gain_mat

allspecbias_genes = c()
allspecbias_bytis = list()
allspecbiaspm_bytis = list()

for (tissue in colnames(bias_matrix_tissep_spec5)){
  print(tissue)
  humanbias_genes = rownames(m.lsfr.spectis.signed.lfccut[which(m.lsfr.spectis.signed.lfccut[,paste(tissue,"Human",sep="_")] != 0),])
  cat(" ")
  cat(length(humanbias_genes))
  cat("\n")
    for (gene in humanbias_genes) {
      spec4_biasvec = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,c("Cyno","Mouse","Rat","Dog"),sep="_")] * m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Human",sep="_")]
      spec4_biasvec[is.nan(spec4_biasvec)] = 0
      #spec4_biasvec[which(spec4_biasvec < 0)] = 0
      if (sum(spec4_biasvec == 0)==4) {
        human_gain_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Human",sep="_")]
      } else if (sum(spec4_biasvec==-1)==4) {
          human_switch_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Human",sep="_")]
      } else if ((sum(spec4_biasvec==0)==3)&(spec4_biasvec[1] == 1)) {
          primate_gain_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Human",sep="_")]
      } else if ((spec4_biasvec[1] == 1)&(spec4_biasvec[2] == 0)&(spec4_biasvec[3] == 0)&(spec4_biasvec[4] == 1)) {
          rodent_loss_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Human",sep="_")]
      }  else if ((sum(spec4_biasvec==1) >= 3)) {
            single_loss_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Human",sep="_")]
      } else if ((sum(spec4_biasvec==0)==3)&((spec4_biasvec[2] == 1)|(spec4_biasvec[3] == 1))) {
            complex_gain_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Human",sep="_")]
      } else if ((sum(spec4_biasvec==0)==2)&(((spec4_biasvec[2] == 1)|(spec4_biasvec[3] == 1)) & (spec4_biasvec[4] == 1) )) {
            complex_anc_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Human",sep="_")]
      } else
            complex_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Human",sep="_")]
    }

  rodentbias_genes = rownames(m.lsfr.spectis.signed.lfccut[which(m.lsfr.spectis.signed.lfccut[,paste(tissue,"Mouse",sep="_")] != 0),])
      for (gene in rodentbias_genes) {
      spec4_biasvec = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,c("Human","Cyno","Rat","Dog"),sep="_")] * m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Mouse",sep="_")]
      #spec4_biasvec[which(spec4_biasvec < 0)] = 0
      spec4_biasvec[is.nan(spec4_biasvec)] = 0
      if ( (sum(spec4_biasvec == 0)==3) & (spec4_biasvec[3] == 1)) {
        rodent_gain_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Mouse",sep="_")]
      } else if ( (sum(spec4_biasvec == 0)==2) & (spec4_biasvec[3] == 1) & (spec4_biasvec[4] == 1)) {
        primate_loss_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Mouse",sep="_")]
      }
      }
  
  cynobias_genes = rownames(m.lsfr.spectis.signed.lfccut[which(m.lsfr.spectis.signed.lfccut[,paste(tissue,"Cyno",sep="_")] != 0),])
      for (gene in cynobias_genes) {
      spec4_biasvec = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,c("Human","Mouse","Rat","Dog"),sep="_")] * m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Cyno",sep="_")]
      spec4_biasvec[is.nan(spec4_biasvec)] = 0
      if (sum(spec4_biasvec == 0)==4) {
        cyno_gain_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Cyno",sep="_")]
      } else if ((sum(spec4_biasvec==1) >= 3)) {
            single_loss_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Cyno",sep="_")]
      } else if ((sum(spec4_biasvec==0)==3)&((spec4_biasvec[2] == 1)|(spec4_biasvec[3] == 1))) {
            complex_gain_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Cyno",sep="_")]
      } else if ((sum(spec4_biasvec==0)==2)&(((spec4_biasvec[2] == 1)|(spec4_biasvec[3] == 1)) & (spec4_biasvec[4] == 1) )) {
            complex_anc_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Cyno",sep="_")]
      } else if ((sum(spec4_biasvec == 0)!=4)&!(gene %in% c(rownames(primate_gain_mat[which(primate_gain_mat[,tissue]!=0),]),rownames(single_loss_mat[which(single_loss_mat[,tissue]!=0),]),rownames(complex_mat[which(complex_mat[,tissue]!=0),]),rownames(rodent_loss_mat[which(rodent_loss_mat[,tissue]!=0),])))){
        if (sum(spec4_biasvec) >= 0) {complex_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Cyno",sep="_")]}
        if (sum(spec4_biasvec) < 0) {complex_mat[gene,tissue] = -m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Cyno",sep="_")]}      
        }
      }
  
  mousebias_genes = rownames(m.lsfr.spectis.signed.lfccut[which(m.lsfr.spectis.signed.lfccut[,paste(tissue,"Mouse",sep="_")] != 0),])
      for (gene in mousebias_genes) {
      spec4_biasvec = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,c("Human","Cyno","Rat","Dog"),sep="_")] * m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Mouse",sep="_")]
      spec4_biasvec[is.nan(spec4_biasvec)] = 0
      if (sum(spec4_biasvec == 0)==4) {
        mouse_gain_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Mouse",sep="_")]
      } else if ((sum(spec4_biasvec==1) >= 3)) {
            single_loss_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Mouse",sep="_")]
      } else if ((sum(spec4_biasvec == 0)!=4)&!(gene %in% c(rownames(primate_loss_mat[which(primate_loss_mat[,tissue]!=0),]),rownames(single_loss_mat[which(single_loss_mat[,tissue]!=0),]),rownames(complex_mat[which(complex_mat[,tissue]!=0),]),rownames(complex_gain_mat[which(complex_gain_mat[,tissue]!=0),]),rownames(complex_anc_mat[which(complex_anc_mat[,tissue]!=0),]),rownames(rodent_gain_mat[which(rodent_gain_mat[,tissue]!=0),])))){
        if (sum(spec4_biasvec) >= 0) {complex_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Mouse",sep="_")]}
        if (sum(spec4_biasvec) < 0) {complex_mat[gene,tissue] = -m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Mouse",sep="_")]}      
      }
      }
  
  ratbias_genes = rownames(m.lsfr.spectis.signed.lfccut[which(m.lsfr.spectis.signed.lfccut[,paste(tissue,"Rat",sep="_")] != 0),])
      for (gene in ratbias_genes) {
      spec4_biasvec = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,c("Human","Cyno","Mouse","Dog"),sep="_")] * m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Rat",sep="_")]
      spec4_biasvec[is.nan(spec4_biasvec)] = 0
      if (sum(spec4_biasvec == 0)==4) {
        rat_gain_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Rat",sep="_")]
      } else if ((sum(spec4_biasvec==1) >= 3)) {
            single_loss_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Rat",sep="_")]
      } else if ((sum(spec4_biasvec == 0)!=4)&!(gene %in% c(rownames(primate_loss_mat[which(primate_loss_mat[,tissue]!=0),]),rownames(single_loss_mat[which(single_loss_mat[,tissue]!=0),]),rownames(complex_mat[which(complex_mat[,tissue]!=0),]),rownames(complex_gain_mat[which(complex_gain_mat[,tissue]!=0),]),rownames(complex_anc_mat[which(complex_anc_mat[,tissue]!=0),]),rownames(rodent_gain_mat[which(rodent_gain_mat[,tissue]!=0),])))){
        if (sum(spec4_biasvec) >= 0) {complex_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Rat",sep="_")]}
        if (sum(spec4_biasvec) < 0) {complex_mat[gene,tissue] = -m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Rat",sep="_")]}      
      }
      }
  
  dogbias_genes = rownames(m.lsfr.spectis.signed.lfccut[which(m.lsfr.spectis.signed.lfccut[,paste(tissue,"Dog",sep="_")] != 0),])
      for (gene in dogbias_genes) {
      spec4_biasvec = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,c("Human","Cyno","Rat","Mouse"),sep="_")] * m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Dog",sep="_")]
      spec4_biasvec[is.nan(spec4_biasvec)] = 0
      if (sum(spec4_biasvec == 0)==4) {
        dog_gain_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Dog",sep="_")]
      } else if ((sum(spec4_biasvec==1) >= 3)) {
            single_loss_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Dog",sep="_")]
      } else if ((sum(spec4_biasvec == 0)!=4)&!(gene %in% c(rownames(primate_loss_mat[which(primate_loss_mat[,tissue]!=0),]),rownames(single_loss_mat[which(single_loss_mat[,tissue]!=0),]),rownames(complex_mat[which(complex_mat[,tissue]!=0),]),rownames(complex_gain_mat[which(complex_gain_mat[,tissue]!=0),]),rownames(complex_anc_mat[which(complex_anc_mat[,tissue]!=0),]),rownames(rodent_loss_mat[which(rodent_loss_mat[,tissue]!=0),])))){
        if (sum(spec4_biasvec) >= 0) {complex_mat[gene,tissue] = m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Dog",sep="_")]}
        if (sum(spec4_biasvec) < 0) {complex_mat[gene,tissue] = -m.lsfr.spectis.signed.lfccut[gene,paste(tissue,"Dog",sep="_")]}      
      }
      }

}

human_gain_mat = human_gain_mat[which(apply(abs(human_gain_mat),1,sum) > 0),]
cyno_gain_mat = cyno_gain_mat[which(apply(abs(cyno_gain_mat),1,sum) > 0),]
mouse_gain_mat = mouse_gain_mat[which(apply(abs(mouse_gain_mat),1,sum) > 0),]
rat_gain_mat = rat_gain_mat[which(apply(abs(rat_gain_mat),1,sum) > 0),]
dog_gain_mat = dog_gain_mat[which(apply(abs(dog_gain_mat),1,sum) > 0),]
human_switch_mat = human_switch_mat[which(apply(abs(human_switch_mat),1,sum) > 0),]
primate_gain_mat = primate_gain_mat[which(apply(abs(primate_gain_mat),1,sum) > 0),]
primate_loss_mat = primate_loss_mat[which(apply(abs(primate_loss_mat),1,sum) > 0),]
rodent_gain_mat = rodent_gain_mat[which(apply(abs(rodent_gain_mat),1,sum) > 0),]
rodent_loss_mat = rodent_loss_mat[which(apply(abs(rodent_loss_mat),1,sum) > 0),]
boroeuth_mat = boroeuth_mat[which(apply(abs(boroeuth_mat),1,sum) > 0),]
single_loss_mat = single_loss_mat[which(apply(abs(single_loss_mat),1,sum) > 0),]
complex_mat = complex_mat[which(apply(abs(complex_mat),1,sum) > 0),]
complex_gain_mat = complex_gain_mat[which(apply(abs(complex_gain_mat),1,sum) > 0),]
complex_anc_mat = complex_anc_mat[which(apply(abs(complex_anc_mat),1,sum) > 0),]

```


Exclude overlaps between lineage-specific sex-biased gene-tissue pairs and conserved sex-biased gene-tissue pairs (few, but simplifies downstream analyses). Produces modified gene x tissue matrices of lineage-specific sex bias ({matrixname}_touse)
```{r}
bias_matrix_touse = bias_matrix_tissep_spec5_nospecint
human_gain_mat_touse = human_gain_mat
for (gene in intersect(rownames(human_gain_mat_touse),rownames(bias_matrix_touse))) {
  for (tis in colnames(human_gain_mat_touse)){
    if (human_gain_mat_touse[gene,tis] * bias_matrix_touse[gene,tis] != 0){
      human_gain_mat_touse[gene,tis] = 0
    }
  }
}
human_gain_mat_touse = human_gain_mat_touse[which(apply(abs(human_gain_mat_touse),1,sum) > 0),]

cyno_gain_mat_touse = cyno_gain_mat
for (gene in intersect(rownames(cyno_gain_mat_touse),rownames(bias_matrix_touse))) {
  for (tis in colnames(cyno_gain_mat_touse)){
    if (cyno_gain_mat_touse[gene,tis] * bias_matrix_touse[gene,tis] != 0){
      cyno_gain_mat_touse[gene,tis] = 0
    }
  }
}
cyno_gain_mat_touse = cyno_gain_mat_touse[which(apply(abs(cyno_gain_mat_touse),1,sum) > 0),]

mouse_gain_mat_touse = mouse_gain_mat
for (gene in intersect(rownames(mouse_gain_mat_touse),rownames(bias_matrix_touse))) {
  for (tis in colnames(mouse_gain_mat_touse)){
    if (mouse_gain_mat_touse[gene,tis] * bias_matrix_touse[gene,tis] != 0){
      mouse_gain_mat_touse[gene,tis] = 0
    }
  }
}
mouse_gain_mat_touse = mouse_gain_mat_touse[which(apply(abs(mouse_gain_mat_touse),1,sum) > 0),]

rat_gain_mat_touse = rat_gain_mat
for (gene in intersect(rownames(rat_gain_mat_touse),rownames(bias_matrix_touse))) {
  for (tis in colnames(rat_gain_mat_touse)){
    if (rat_gain_mat_touse[gene,tis] * bias_matrix_touse[gene,tis] != 0){
      rat_gain_mat_touse[gene,tis] = 0
    }
  }
}
rat_gain_mat_touse = rat_gain_mat_touse[which(apply(abs(rat_gain_mat_touse),1,sum) > 0),]

dog_gain_mat_touse = dog_gain_mat
for (gene in intersect(rownames(dog_gain_mat_touse),rownames(bias_matrix_touse))) {
  for (tis in colnames(dog_gain_mat_touse)){
    if (dog_gain_mat_touse[gene,tis] * bias_matrix_touse[gene,tis] != 0){
      dog_gain_mat_touse[gene,tis] = 0
    }
  }
}
dog_gain_mat_touse = dog_gain_mat_touse[which(apply(abs(dog_gain_mat_touse),1,sum) > 0),]

primate_gain_mat_touse = primate_gain_mat
for (gene in intersect(rownames(primate_gain_mat_touse),rownames(bias_matrix_touse))) {
  for (tis in colnames(primate_gain_mat_touse)){
    if (primate_gain_mat_touse[gene,tis] * bias_matrix_touse[gene,tis] != 0){
      primate_gain_mat_touse[gene,tis] = 0
    }
  }
}
primate_gain_mat_touse = primate_gain_mat_touse[which(apply(abs(primate_gain_mat_touse),1,sum) > 0),]


rodent_gain_mat_touse = rodent_gain_mat
for (gene in intersect(rownames(rodent_gain_mat_touse),rownames(bias_matrix_touse))) {
  for (tis in colnames(rodent_gain_mat_touse)){
    if (rodent_gain_mat_touse[gene,tis] * bias_matrix_touse[gene,tis] != 0){
      rodent_gain_mat_touse[gene,tis] = 0
    }
  }
}
rodent_gain_mat_touse = rodent_gain_mat_touse[which(apply(abs(rodent_gain_mat_touse),1,sum) > 0),]

complex_anc_mat_touse = complex_anc_mat
for (gene in intersect(rownames(complex_anc_mat_touse),rownames(bias_matrix_touse))) {
  for (tis in colnames(complex_anc_mat_touse)){
    if (complex_anc_mat_touse[gene,tis] * bias_matrix_touse[gene,tis] != 0){
      complex_anc_mat_touse[gene,tis] = 0
    }
  }
}
complex_anc_mat_touse = complex_gain_mat_touse[which(apply(abs(complex_gain_mat_touse),1,sum) > 0),]

complex_gain_mat_touse = complex_gain_mat
for (gene in intersect(rownames(complex_gain_mat_touse),rownames(bias_matrix_touse))) {
  for (tis in colnames(complex_gain_mat_touse)){
    if (complex_gain_mat_touse[gene,tis] * bias_matrix_touse[gene,tis] != 0){
      complex_gain_mat_touse[gene,tis] = 0
    }
  }
}
complex_gain_mat_touse = complex_gain_mat_touse[which(apply(abs(complex_gain_mat_touse),1,sum) > 0),]

complex_mat_touse = complex_mat
for (gene in intersect(rownames(complex_mat_touse),rownames(bias_matrix_touse))) {
  for (tis in colnames(complex_mat_touse)){
    if (complex_mat_touse[gene,tis] * bias_matrix_touse[gene,tis] != 0){
      complex_mat_touse[gene,tis] = 0
    }
  }
}
complex_mat_touse = complex_mat_touse[which(apply(abs(complex_mat_touse),1,sum) > 0),]

rodent_loss_mat_touse = rodent_loss_mat
for (gene in intersect(rownames(rodent_loss_mat_touse),rownames(bias_matrix_touse))) {
  for (tis in colnames(rodent_loss_mat_touse)){
    if (rodent_loss_mat_touse[gene,tis] * bias_matrix_touse[gene,tis] != 0){
      rodent_loss_mat_touse[gene,tis] = 0
    }
  }
}
rodent_loss_mat_touse = rodent_loss_mat_touse[which(apply(abs(rodent_loss_mat_touse),1,sum) > 0),]

primate_loss_mat_touse = primate_loss_mat
for (gene in intersect(rownames(primate_loss_mat_touse),rownames(bias_matrix_touse))) {
  for (tis in colnames(primate_loss_mat_touse)){
    if (primate_loss_mat_touse[gene,tis] * bias_matrix_touse[gene,tis] != 0){
      primate_loss_mat_touse[gene,tis] = 0
    }
  }
}
primate_loss_mat_touse = primate_loss_mat_touse[which(apply(abs(primate_loss_mat_touse),1,sum) > 0),]

```

Write gene x tissue matrices of each type of sex bias to file
```{r}
write.table(bias_matrix_tissep_spec5_nospecint,"../sexbias.conserved.matrix.txt",col.names = TRUE,row.names = TRUE,sep="\t",quote=FALSE)
write.table(human_gain_mat_touse,"../sexbias.humangain.matrix.txt",col.names = TRUE,row.names = TRUE,sep="\t",quote=FALSE)
write.table(cyno_gain_mat_touse,"../sexbias.cynogain.matrix.txt",col.names = TRUE,row.names = TRUE,sep="\t",quote=FALSE)
write.table(mouse_gain_mat_touse,"../sexbias.mousegain.matrix.txt",col.names = TRUE,row.names = TRUE,sep="\t",quote=FALSE)
write.table(rat_gain_mat_touse,"../sexbias.ratgain.matrix.txt",col.names = TRUE,row.names = TRUE,sep="\t",quote=FALSE)
write.table(dog_gain_mat_touse,"../sexbias.doggain.matrix.txt",col.names = TRUE,row.names = TRUE,sep="\t",quote=FALSE)
write.table(primate_gain_mat_touse,"../sexbias.primategain.matrix.txt",col.names = TRUE,row.names = TRUE,sep="\t",quote=FALSE)
write.table(rodent_gain_mat_touse,"../sexbias.rodentgain.matrix.txt",col.names = TRUE,row.names = TRUE,sep="\t",quote=FALSE)
write.table(complex_gain_mat_touse,"../sexbias.multiplegain.matrix.txt",col.names = TRUE,row.names = TRUE,sep="\t",quote=FALSE)
write.table(complex_anc_mat_touse,"../sexbias.multipleloss.matrix.txt",col.names = TRUE,row.names = TRUE,sep="\t",quote=FALSE)
write.table(complex_mat_touse,"../sexbias.complex.matrix.txt",col.names = TRUE,row.names = TRUE,sep="t",quote=FALSE)
write.table(rodent_loss_mat_touse,"../sexbias.rodentloss.matrix.txt",col.names = TRUE,row.names = TRUE,sep="\t",quote=FALSE)
write.table(primate_loss_mat_touse,"../sexbias.primateloss.matrix.txt",col.names = TRUE,row.names = TRUE,sep="\t",quote=FALSE)

```

