---
title: "annotating larval geoduck transcriptome (for ariana)"
author: "Olivia Cattau"
date: "2/24/2022"
output: pdf_document
---
#P.generosa (Pacific Geoduck) combined early stages and late stages transcriptomes for analyses and annotations 
#load count data from kallisto output(s)
```{r}
countMatrix<-read.table(file="https://raw.githubusercontent.com/sr320/nb-2022/main/P_generosa/analyses/kallisto-0207.isoform.counts.matrix", header=TRUE, sep = '\t')
names(countMatrix)[1]<-"target_id"
head(countMatrix)
```

#use binary data to make Venn Diagram
```{r}
library(tidyr)
library(dplyr)
simple_data<-countMatrix
long_data<- reshape(data=simple_data, idvar="target_id", varying=c("heart", "gonad", "ctenidia", "juv_sl", "juv_amb", "larvae"), v.name=c("tpm"), times=c("heart", "gonad", "ctenidia", "juv_sl", "juv_amb", "larvae"), new.row.names=1:8183754, direction="long")
names(long_data)[2]<-"tissue"
summary(long_data)
long_data$binary=ifelse(long_data$tpm > 0, 1, 0)
long_data_zeros<-long_data[apply(long_data, 1, function(row) all(row !=0)),]
df<-long_data_zeros[-c(3,4)] #tissue type and gene names only
table_no_zeros<-as.data.frame(crossprod(table(df)))
#write.csv(table_no_zeros, file="/Users/oliviacattau/Documents/GitHub/code/characterize_larval_transciptome/pgenerosa_table_no_zeros.tab")
simple_data_cleaned<-simple_data[apply(simple_data, 1, function(row) all(row !=0)),]
genelist<-long_data_zeros %>% group_by(tissue)
output<-genelist %>% summarise("total reads (n)"=sum(binary),  "mean tpm" = mean(tpm), "Q4 tpm" = quantile(tpm, c(0.75)))
#write.csv(output, file="/Users/oliviacattau/Documents/GitHub/code/characterize_larval_transciptome/genelist_02_24_22.tab")
```

#AH's tidy alternative approach: 
```{r}
library(devtools)
install_github("jokergoo/ComplexHeatmap")
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(version = "3.14")
library(BiocManager)
library(tidyr)
library(dplyr)
library(UpSetR)
library(ComplexHeatmap)
tidy_data<-countMatrix
head(tidy_data)
#now add the binary column
tidy_data$h=ifelse(tidy_data$heart > 0, 1, 0)
tidy_data$g=ifelse(tidy_data$gonad >0, 1, 0)
tidy_data$c=ifelse(tidy_data$ctenidia > 0, 1, 0)
tidy_data$l=ifelse(tidy_data$larvae > 0, 1, 0)
tidy_data$ja=ifelse(tidy_data$juv_amb > 0, 1, 0)
tidy_data$jsl=ifelse(tidy_data$juv_sl > 0, 1, 0)
#check code
heart_length<-sum(tidy_data$h)
gonad_length<-sum(tidy_data$g)
ctenidia_length<-sum(tidy_data$c)
larvae_length<-sum(tidy_data$l)
juv_amb_length<-sum(tidy_data$ja)
juv_sl_length<-sum(tidy_data$jsl)
UpSet_test<-data.frame(tissue=c(heart_length, gonad_length, ctenidia_length, larvae_length, juv_amb_length, juv_sl_length))
rownames(UpSet_test)<-c("heart", "gonad", "ctenidia", "larvae", "juv ambient", "juv low oa")
#for making upset plot, must install above packages. Same code works fine on my desktop version but is struggling to knit on the Raven server
m1 = make_comb_mat(heart=tidy_data$heart, gonad=tidy_data$gonad, ctenidia=tidy_data$ctenidia, juv_amb=tidy_data$juv_amb, juv_sl=tidy_data$juv_sl, larvae=tidy_data$larvae) #only distinct intersection 
m2 = make_comb_mat(heart=tidy_data$heart, gonad=tidy_data$gonad, ctenidia=tidy_data$ctenidia, juv_amb=tidy_data$juv_amb, juv_sl=tidy_data$juv_sl, larvae=tidy_data$larvae, mode="intersect") #only overlapping intersections (true Venn diagram) see https://jokergoo.github.io/ComplexHeatmap-reference/book/upset-plot.html for explanation 
set_size(m1)
comb_size(m1)
comb_size(m2)
p1<-UpSet(m1 [comb_degree(m1)<=2], column_title ="All intersections of 2 or less") #intersections 2 or less
p2<-UpSet(m1 [comb_degree(m1)>=3], column_title= "All intersection of 3 or more") #intersections 3 or more
p3<-UpSet(m2, top_annotation = HeatmapAnnotation(
    degree = as.character(comb_degree(m2)),
    "Intersection\nsize" = anno_barplot(comb_size(m2), 
        border = FALSE, 
        gp = gpar(fill = "black"), 
        height = unit(2, "cm")
    ), 
    annotation_name_side = "left", 
    annotation_name_rot = 0), column_title="p. generosa transcriptome intersecting intersections") #intersect mode
p4<-UpSet(m1, top_annotation = HeatmapAnnotation(
    degree = as.character(comb_degree(m1)),
    "Intersection\nsize" = anno_barplot(comb_size(m1), 
        border = FALSE, 
        gp = gpar(fill = "black"), 
        height = unit(2, "cm")
    ), 
    annotation_name_side = "left", 
    annotation_name_rot = 0), column_title="p. generosa transcriptome intersecting intersections")
#####
tidy_data_long<-tidy_data %>%
  gather(key="time", value="tpm", heart:larvae) #name key as time, name value as tpm, gather columns heart:larvae
tidy_data_long$binary=ifelse(tidy_data_long$tpm > 0, 1, 0)
head(tidy_data_long)
#filter out zeros
tidy_data_long_no_zeros<-tidy_data_long%>%
  filter(!binary==0)
#filter zeros out of wide format dataframe 
tidy_data_wide_no_zeros<-tidy_data%>%
  mutate(sum = rowSums(across(where(is.numeric))))%>% #sum across numeric rows
  filter(!sum==0) #remove rows with sum=0 
#generate summary
genelist<-tidy_data_long_no_zeros %>% 
  group_by(time)%>%
  summarise(TotalReadsN=sum(binary), meanTPM = mean(tpm), Q4tpm = quantile(tpm, c(0.75)))
genelist
dataUpset<-tidy_data_long_no_zeros %>% group_by(time) 
dataUpset %>% mutate(OnlyLarvae=sum(larvae), LarvaeGonad=sum(larvae,gonad))
```