DiffIsoAnalysis.Rmd

---
title: "Total v Nuclear"
author: "Briana Mittleman"
date: "5/2/2019"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


```{r}
library(workflowr)
library(tidyverse)
library(reshape2)
```

In this analysis I wil use leafcutter to call PAS with differential ussage between fractions.

##Prepare annotation  

I first filter the annotated peak SAF file for peaks passing the 5% coverage in either fraction.  

```{bash,eval=F}
python makeSAFbothfrac5perc.py
```


##Peak quantification
```{bash,eval=F}
mkdir bothFrac_FC
```

Run feature counts with these peaks with both fractions:  

```{bash,eval=F}
sbatch bothFrac_FC.sh
```


Fix the header:
```{bash,eval=F}
python fixFChead_bothfrac.py ../data/bothFrac_FC/APApeaks.ALLChrom.Filtered.Named.GeneLocAnnoPARSED.5percCov.bothfrac.fc ../data/bothFrac_FC/APApeaks.ALLChrom.Filtered.Named.GeneLocAnnoPARSED.5percCov.bothfrac.fixed.fc
```

Remove location demoniaiton:  


#Prepare leafcutter phenotype  

```{bash,eval=F}
mkdir ../data/DiffIso
```


```{bash,eval=F}
python fc2leafphen.py
```

Fix pheno to remove location:

```{bash,eval=F}
python removeloc_pheno.py ../data/DiffIso/APApeaks.ALLChrom.Filtered.Named.GeneLocAnnoPARSED.5percCov.bothfrac.fixed.forLC.fc ../data/DiffIso/APApeaks.ALLChrom.Filtered.Named.GeneLocAnnoPARSED.5percCov.bothfrac.fixed.forLC_noloc.fc
```


```{bash,eval=F}
python subset_diffisopheno.py 1
python subset_diffisopheno.py 2
python subset_diffisopheno.py 3
python subset_diffisopheno.py 4
python subset_diffisopheno.py 5
python subset_diffisopheno.py 6
python subset_diffisopheno.py 7
python subset_diffisopheno.py 8
python subset_diffisopheno.py 9
python subset_diffisopheno.py 10
python subset_diffisopheno.py 11
python subset_diffisopheno.py 12
python subset_diffisopheno.py 13
python subset_diffisopheno.py 14
python subset_diffisopheno.py 15
python subset_diffisopheno.py 16
python subset_diffisopheno.py 18
python subset_diffisopheno.py 19
python subset_diffisopheno.py 20
python subset_diffisopheno.py 21
python subset_diffisopheno.py 22
```

Make the sample groups file: 

```{bash,eval=F}
python LC_samplegroups.py 
```

##Run leafcutter  
The leafcutter environment is not in the three-prime-seq environment. Make sure leafcutter is installed and working.  
```{bash,eval=F}
sbatch run_leafcutterDiffIso.sh
```

Rscript /project2/gilad/briana/davidaknowles-leafcutter-c3d9474/scripts/leafcutter_ds.R --num_threads 4  ../data/DiffIso/APApeaks.ALLChrom.Filtered.Named.GeneLocAnnoPARSED.5percCov.bothfrac.fixed.forLC.fc_22.txt ../data/bothFrac_FC/sample_groups.txt -o ../data/DiffIso/TN_diff_isoform_chr22.txt


Concatinate results:

```{bash,eval=F}
awk '{if(NR>1)print}' ../data/DiffIso/TN_diff_isoform_chr*.txt_effect_sizes.txt > ../data/DiffIso/TN_diff_isoform_allChrom.txt_effect_sizes.txt


awk '{if(NR>1)print}' ../data/DiffIso/TN_diff_isoform_chr*.txt_cluster_significance.txt > ../data/DiffIso/TN_diff_isoform_AllChrom_cluster_significance.txt
```


##Evaluate results
###Significant clusters

```{r}
sig=read.table("../data/DiffIso/TN_diff_isoform_AllChrom_cluster_significance.txt",sep="\t" ,col.names = c('status','loglr','df','p','cluster','p.adjust'),stringsAsFactors = F) %>% filter(status=="Success") 

sig$p.adjust=as.numeric(as.character(sig$p.adjust))

```

```{r}
qqplot(-log10(runif(nrow(sig))), -log10(sig$p.adjust),ylab="-log10 Total Adjusted Leafcutter pvalue", xlab="-log 10 Uniform expectation", main="Leafcutter differencial isoform analysis between fractions")
abline(0,1)
```
```{r}
tested_genes=nrow(sig)
tested_genes
```



```{r}
sig_genes=sig %>% filter(p.adjust<.05)
number_sig_genes=nrow(sig_genes)
number_sig_genes

sig_genesonly=sig_genes %>% separate(cluster,into=c("chrom", "geneName"), sep = ":") %>% dplyr::select(geneName)

write.table(sig_genesonly, file="../data/sigDiffGenes.txt", col.names = T, row.names = F, quote = F)
```

###Effect sizes 

```{r}

effectsize=read.table("../data/DiffIso/TN_diff_isoform_allChrom.txt_effect_sizes.txt", stringsAsFactors = F, col.names=c('intron',  'logef' ,'Nuclear', 'Total','deltaPAU')) %>% filter(intron != "intron")

write.table(effectsize,file="../data/DiffIso/EffectSizes.txt", quote = F, col.names = T, row.names = F)

effectsize$deltaPAU=as.numeric(as.character(effectsize$deltaPAU))
effectsize$logef=as.numeric(as.character(effectsize$logef))
```

Plot delta PAU:  
```{r}
plot(sort(effectsize$deltaPAU),main="Leafcutter delta PAU", ylab="Delta PAU", xlab="PAS Index")
```

Filter PAU > .2

```{r}
effectsize_deltaPAU= effectsize %>% filter(abs(deltaPAU) > .2) 
nrow(effectsize_deltaPAU)

effectSize_highdiffGenes=effectsize_deltaPAU %>% separate(intron, into=c("chrom", "start", "end", "GeneName"), sep=":") %>% dplyr::select(GeneName) %>% unique()


write.table(effectSize_highdiffGenes, file="../data/highdiffsiggenes.txt", col.names = F, row.names = F, quote = F)
```
Genes in this set:  

```{r}
effectsize_deltaPAU_Genes= effectsize_deltaPAU %>% separate(intron, into=c("chrom", "start", "end","gene"),sep=":") %>% group_by(gene) %>% summarise(nperGene=n()) 

nrow(effectsize_deltaPAU_Genes)
```

Filter >.2 in  

```{r}
effectsize_deltaPAU_nuclear= effectsize_deltaPAU %>% filter(deltaPAU < -0.2)

#write out at bed
#need strand info
PAS=read.table("../data/PAS/APAPAS_GeneLocAnno.5perc.bed", stringsAsFactors = F,col.names = c("chrom", "start", "end", "peak", "score", "strand") )%>% separate(peak, into=c("peaknum","peakID"), sep=":") %>% separate(peakID, into=c("gene", "loc"), sep="_") %>% dplyr::select(gene, strand) %>% unique()
effectsize_deltaPAU_nuclear_bed=effectsize_deltaPAU_nuclear %>% separate(intron, into=c("chr", "peakStart", "peakEnd", "gene"), sep=":") %>% inner_join(PAS, by="gene")  %>% mutate(PASstart=ifelse(strand=="+", as.integer(peakEnd)-1, as.integer(peakStart)+1)) %>% mutate(PASend=ifelse(strand=="+", as.integer(peakEnd), as.integer(peakStart))) %>% mutate(score=".") %>%  dplyr::select(chr, peakStart, peakEnd, gene, score, strand) 

write.table(effectsize_deltaPAU_nuclear_bed, file="../data/PAS/UsedMoreNuclearPAU2.bed", col.names = F, row.names = F, quote = F,sep = "\t")
```

Filter >.2 in Total: 

```{r}
effectsize_deltaPAU_total= effectsize_deltaPAU %>% filter(deltaPAU > 0.2)

effectsize_deltaPAU_total_bed=effectsize_deltaPAU_total %>% separate(intron, into=c("chr", "peakStart", "peakEnd", "gene"), sep=":") %>% inner_join(PAS, by="gene")  %>% mutate(PASstart=ifelse(strand=="+", as.integer(peakEnd)-1, as.integer(peakStart)+1)) %>% mutate(PASend=ifelse(strand=="+", as.integer(peakEnd), as.integer(peakStart))) %>% mutate(score=".") %>%  dplyr::select(chr, peakStart, peakEnd, gene, score, strand) 

write.table(effectsize_deltaPAU_total_bed, file="../data/PAS/UsedMoreTotalPAU2.bed", col.names = F, row.names = F, quote = F,sep="\t")
```

Sort the files:  

```{bash,eval=F}
sort -k1,1 -k2,2n ../data/PAS/UsedMoreTotalPAU2.bed > ../data/PAS/UsedMoreTotalPAU2.sort.bed
sort -k1,1 -k2,2n ../data/PAS/UsedMoreNuclearPAU2.bed > ../data/PAS/UsedMoreNuclearPAU2.sort.bed
```

##Location of high >PAU  

###Total:

Pull in location information for each PAS:  

```{r}
PAS=read.table("../data/peaks_5perc/APApeaks.ALLChrom.Filtered.Named.GeneLocAnnoPARSED.5percCov.bothfrac.SAF",stringsAsFactors = F,header = T) %>% separate(GeneID, into=c("num", "chr", "start", "end", "strand", "geneID"), sep=":") %>% separate(geneID, into=c("gene", "loc"),sep="_") %>%  mutate(intron=paste("chr", Chr, ":", Start, ":", End, ":", gene,sep="")) %>% select(intron, loc)
```



```{r}
effectsize_deltaPAU_total_loc=effectsize_deltaPAU_total %>% inner_join(PAS, by="intron") 


ggplot(effectsize_deltaPAU_total_loc,aes(x=loc)) + geom_histogram(stat="count") + labs(title="Location of Total peaks >.2 PAU") 
```

###Nuclear:


```{r}
effectsize_deltaPAU_nuclear_loc=effectsize_deltaPAU_nuclear %>% inner_join(PAS, by="intron") 


ggplot(effectsize_deltaPAU_nuclear_loc,aes(x=loc)) + geom_histogram(stat="count") + labs(title="Location of Nuclear peaks >.2 PAU")
```

I will want to look at proportions. I need to know how many peaks are in each location: 

```{r}
PAS_loc =PAS%>% group_by(loc) %>% summarise(nloc=n())

```


```{r}
effectsize_deltaPAU_total_locProp=effectsize_deltaPAU_total_loc %>% group_by(loc) %>% summarise(nloctotal=n()) 
effectsize_deltaPAU_nuclear_locProp=effectsize_deltaPAU_nuclear_loc %>% group_by(loc) %>% summarise(nlocnuclear=n()) 

effectsize_deltaPAUProp_tot=effectsize_deltaPAU_total_locProp %>% inner_join(PAS_loc, by="loc") %>% mutate(Proportion_tot=nloctotal/nloc)

effectsize_deltaPAUProp_nuc=effectsize_deltaPAU_nuclear_locProp %>% inner_join(PAS_loc, by="loc") %>% mutate(Proportion_nuc=nlocnuclear/nloc)

```

```{r}
ggplot(effectsize_deltaPAUProp_tot, aes(x=loc, y=Proportion_tot)) + geom_bar(stat="identity") + labs(y="Proportion of all called PAS", title="Location of high Total used PAS")
```

```{r}
ggplot(effectsize_deltaPAUProp_nuc, aes(x=loc, y=Proportion_nuc)) + geom_bar(stat="identity") + labs(y="Proportion of all called PAS", title="Location of high nuclear used PAS")
```

Merge to 1 figure: 

```{r}
effectsize_deltaPAUProp_both= effectsize_deltaPAUProp_nuc %>% inner_join(effectsize_deltaPAUProp_tot, by=c("loc","nloc")) %>% dplyr::rename(Nuclear=Proportion_nuc, Total=Proportion_tot) %>% select(loc, Nuclear, Total) 
effectsize_deltaPAUProp_both_melt= effectsize_deltaPAUProp_both %>% melt(id.vars="loc", variable.name="Fraction", value.name = "Proportion") 
effectsize_deltaPAUProp_both_melt$Fraction=as.character(effectsize_deltaPAUProp_both_melt$Fraction)
```

```{r}
ggplot(effectsize_deltaPAUProp_both_melt, aes(x=loc, y=Proportion, by=Fraction, fill=Fraction)) + geom_bar(stat="identity", position="dodge") + scale_fill_manual(values=c("deepskyblue3","darkviolet")) + labs(title="Proportion of PAS differential used by location",x="") +scale_x_discrete(labels = c('Coding','5kb downstream','Intronic',"3' UTR", "5' UTR")) +theme(axis.text.x = element_text(angle = 90, hjust = 1)) +  theme(legend.position = c(0.1,.9), legend.direction = "horizontal") +  theme(panel.background = element_blank()) 
```

```{r}
effectsize_deltaPAU_total_locProp

sum(effectsize_deltaPAU_total_locProp$nloctotal)

effectsize_deltaPAU_nuclear_locProp
sum(effectsize_deltaPAU_nuclear_locProp$nlocnuclear)
```

```{r}
effectsize_deltaPAUProp_both_melt_sm=effectsize_deltaPAUProp_both_melt %>% filter(loc=="intron" | loc=="utr3")


ggplot(effectsize_deltaPAUProp_both_melt_sm, aes(x=loc, y=Proportion, by=Fraction, fill=Fraction)) + geom_bar(stat="identity", position="dodge") + scale_fill_manual(values=c("deepskyblue3","darkviolet")) + labs(title="Proportion of PAS differential used by location",x="") +scale_x_discrete(labels = c('Intronic',"3' UTR")) +theme(axis.text.x = element_text(angle = 90, hjust = 1)) +  theme(legend.position = c(0.1,.9), legend.direction = "horizontal") +  theme(panel.background = element_blank()) 
```


```{r}
#intronic
prop.test(x=c(473,160), n=c(700,1396),alternative = "greater")

#3' utr
prop.test(x=c(104,1085), n=c(700,1396),alternative = "less")
```


More differentiall used in total. this makes sense because there are more used peaks in the nuclear which evens out the distribution of the ratios.  

###Stratify by different $\Delta$ PAU

I want to create a data frame that has the location proportion distribution based on different $\Delta$ PAU. 0-.1 .1-.2 .2-.3 .3-.4 .4-.5 >.5

First I will seperate the total and nuclear but the sign of the $\Delta$ PAU.  

```{r}
colnames(effectsize)=c("intron", "logef","Nuclear", "Total", "deltaPAU")
Total_dpau= effectsize %>% filter(deltaPAU > 0) %>% inner_join(PAS, by="intron") %>% select(-logef, -Nuclear,-Total) %>%  mutate(fraction="Total", PAU_Cat=ifelse(deltaPAU <.1, "<.1", ifelse(deltaPAU >=.1 & deltaPAU <.2, "<.2", ifelse(deltaPAU >=.2 & deltaPAU <.3, "<.3", ifelse(deltaPAU >=.3 & deltaPAU <.4, "<.4", "<.5"))))) 

Nuclear_dpau= effectsize %>% filter(deltaPAU <0) %>% inner_join(PAS, by="intron") %>% select(-logef,-Nuclear,-Total) %>% mutate(fraction="Nuclear", PAU_Cat=ifelse(deltaPAU >-.1, "<.1", ifelse(deltaPAU <=-.1 & deltaPAU > -.2, "<.2", ifelse(deltaPAU <=-.2 & deltaPAU >-.3, "<.3", ifelse(deltaPAU <=-.3 & deltaPAU >-.4, "<.4", "<.5")))))
```


Merge these together to start grouping: 

```{r}
allPAU=as.data.frame(rbind(Total_dpau, Nuclear_dpau)) %>% group_by(fraction, PAU_Cat, loc ) %>% summarise(nperLoc=n()) %>% full_join(PAS_loc, by ="loc") %>% mutate(Prop=nperLoc/nloc)
```

Plot it:

```{r}
ggplot(allPAU, aes(x=loc,y=Prop, group=fraction, fill=fraction)) + geom_bar(stat="identity", position = "dodge") + facet_wrap(~PAU_Cat)+ scale_fill_manual(values=c("deepskyblue3","darkviolet")) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + labs(title="Proportion of PAS by location and delta PAU")
```
```{r}
allPAU_remove.1= allPAU %>% filter(PAU_Cat != "<.1")

ggplot(allPAU_remove.1, aes(x=loc,y=Prop, group=fraction, fill=fraction)) + geom_bar(stat="identity", position = "dodge") + facet_wrap(~PAU_Cat)+ scale_fill_manual(values=c("deepskyblue3","darkviolet")) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + labs(title="Proportion of PAS by location and delta PAU")
```

Proportion within group:


```{r}
allPAU_ingroup= allPAU %>% mutate(nCat=sum(nperLoc),proppercat=nperLoc/nCat)

ggplot(allPAU_ingroup, aes(x=loc,y=proppercat, group=fraction, fill=fraction)) + geom_bar(stat="identity", position = "dodge") + facet_wrap(~PAU_Cat)+ scale_fill_manual(values=c("deepskyblue3","darkviolet")) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + labs(title="Proportion of PAS by location and delta PAU")
```


##Distance to TSS:

I need to pull in the TSS information so I can look at the distance between the differentially used peaks and by distance .

```{r}
tss=read.table("../../genome_anotation_data/refseq.ProteinCoding.bed",col.names = c("chrom", "start", "end", "gene", "score", "strand") ,stringsAsFactors = F) %>% mutate(TSS= ifelse(strand=="+", start, end)) %>% select(gene, TSS, strand)
```


Seperate effect size introns:

PAS base for + strand is end, PAS for neg stand in -
```{r}
effectsize_TSS= effectsize %>% separate(intron, into=c("chrom", "start", "end", "gene"),sep=":") %>% mutate(fraction=ifelse(deltaPAU < 0, "nuclear", "total")) %>% inner_join(tss, by="gene") %>% mutate(dist2PAS=ifelse(strand=="+", as.numeric(end)-as.numeric(TSS), as.numeric(TSS)-as.numeric(start))) 

effectsize_TSS_tot= effectsize_TSS %>% filter(fraction=="total") %>% mutate( PAU_Cat=ifelse(deltaPAU <.1, "<.1", ifelse(deltaPAU >=.1 & deltaPAU <.2, "<.2", ifelse(deltaPAU >=.2 & deltaPAU <.3, "<.3", ifelse(deltaPAU >=.3 & deltaPAU <.4, "<.4", "<.5"))))) 


effectsize_TSS_nuc=effectsize_TSS %>% filter(fraction=="nuclear") %>% mutate( PAU_Cat=ifelse(deltaPAU >-.1, "<.1", ifelse(deltaPAU <=-.1 & deltaPAU > -.2, "<.2", ifelse(deltaPAU <=-.2 & deltaPAU >-.3, "<.3", ifelse(deltaPAU <=-.3 & deltaPAU >-.4, "<.4", "<.5")))))


effectsize_TSS_cat=as.data.frame(rbind(effectsize_TSS_tot, effectsize_TSS_nuc)) %>% filter(dist2PAS >0)
```


```{r}
ggplot(effectsize_TSS_cat, aes(x=log10(dist2PAS), by=fraction, fill=fraction))+ geom_density(alpha=.4) + facet_grid(~PAU_Cat) + labs(title="Distance to TSS for differentialy used PAS")+scale_fill_manual(values=c("deepskyblue3","darkviolet")) 
```

###By length of gene  


```{r}
length=read.table("../../genome_anotation_data/refseq.ProteinCoding.bed",col.names = c("chrom", "start", "end", "gene", "score", "strand") ,stringsAsFactors = F) %>% mutate(length=abs(end-start)) %>%  mutate(TSS= ifelse(strand=="+", start, end)) %>% select(gene, length,TSS, strand)
```


```{r}
effectsize_length= effectsize %>% separate(intron, into=c("chrom", "start", "end", "gene"),sep=":") %>% mutate(fraction=ifelse(deltaPAU < 0, "nuclear", "total")) %>% inner_join(length, by="gene") %>% mutate(PercLength=ifelse(strand=="+", ((as.numeric(end)-as.numeric(TSS))/as.numeric(length)), (1-(as.numeric(start)-as.numeric(TSS))/as.numeric(length)))) 

effectsize_length_tot= effectsize_length %>% filter(fraction=="total") %>% mutate( PAU_Cat=ifelse(deltaPAU <.1, "<.1", ifelse(deltaPAU >=.1 & deltaPAU <.2, "<.2", ifelse(deltaPAU >=.2 & deltaPAU <.3, "<.3", ifelse(deltaPAU >=.3 & deltaPAU <.4, "<.4", "<.5"))))) 


effectsize_length_nuc=effectsize_length %>% filter(fraction=="nuclear") %>% mutate( PAU_Cat=ifelse(deltaPAU >-.1, "<.1", ifelse(deltaPAU <=-.1 & deltaPAU > -.2, "<.2", ifelse(deltaPAU <=-.2 & deltaPAU >-.3, "<.3", ifelse(deltaPAU <=-.3 & deltaPAU >-.4, "<.4", "<.5")))))


effectsize_length_cat=as.data.frame(rbind(effectsize_length_tot, effectsize_length_nuc)) %>% filter(PercLength<=1 & PercLength >0)

effectsize_length_catall=as.data.frame(rbind(effectsize_length_tot, effectsize_length_nuc)) 
```

```{r}
ggplot(effectsize_length_cat, aes(x=PercLength, by=fraction, fill=fraction))+ geom_histogram(alpha=.4,bins=10) + facet_grid(~PAU_Cat) + labs(title="Location of differentially used PAS within a gene body ")+scale_fill_manual(values=c("deepskyblue3","darkviolet")) 
```
```{r}
summary(effectsize_length_catall$PercLength)
```


```{r}
summary(effectsize$logef)

```

```{r}
ggplot(effectsize_length_cat, aes(x=PercLength, by=fraction, fill=fraction))+ geom_histogram(,bins=50)  + labs(title="Location of differentially used PAS \nwithin a gene body", fill="Fraction", y="Number of PAS", x="Percent of Gene Length")+scale_fill_manual(values=c("deepskyblue3","darkviolet"),labels = c("Nuclear", "Total"))+ theme(legend.position = c(0.1,.9), legend.direction = "horizontal")+  theme(panel.background = element_blank())
```
```{r}
ggplot(effectsize_length_cat, aes(x=PercLength, by=fraction, fill=fraction))+ geom_density(alpha=.5)  + labs(title="Location of differentially used PAS \nwithin a gene body", fill="Fraction", x="Percent of Gene Length")+scale_fill_manual(values=c("deepskyblue3","darkviolet"),labels = c("Nuclear", "Total"))+ theme(legend.position = c(0.1,.9), legend.direction = "horizontal")+  theme(panel.background = element_blank())
```

Diff iso gene proportion:

```{r}
genes_sig=sig %>% separate(cluster,into=c("chr", "gene"), sep=":") %>% group_by(gene) %>% summarise(n=n()) %>% nrow
genes_detlapau= effectSize_highdiffGenes %>% nrow()
testedgenes=read.table("../data/DiffIso/APApeaks.ALLChrom.Filtered.Named.GeneLocAnnoPARSED.5percCov.bothfrac.fixed.forLC.fc",header = T, stringsAsFactors = F) %>% rownames_to_column("ID") %>% select(ID)%>% separate(ID, into=c("chr", "start", "end", "geneID"),sep=":") %>% separate(geneID, into=c("gene", "loc"),sep="_")  %>% group_by(gene) %>% summarise(n=n()) %>% nrow()
notsig=testedgenes-genes_sig
sighothighpau=genes_sig-genes_detlapau

cat=c("NotSig", "SigNotHighPAU", "SigandHighPAU")
values=c(unlist(notsig),unlist(sighothighpau),unlist(genes_detlapau))

difiso_df=as.data.frame(cbind(cat, values)) 
difiso_df$values=as.numeric(as.character(difiso_df$values))
difiso_df=difiso_df%>% mutate(proportion=values/testedgenes)


ggplot(difiso_df, aes(x="",y=proportion, fill=cat)) + geom_bar(stat="identity")+geom_text(aes(label=values))
```


```{r}
slices <- c(notsig, sighothighpau,genes_detlapau)
lbls <- c("No Sig PAS", "At least 1 \nSig PAS", "At least 1 Sig PAS\n High Delta PAU")
pct <- round(slices/sum(slices)*100)
lbls <- paste(lbls, pct, sep="\n   ") # add percents to labels 
lbls <- paste(lbls,"%",sep="") # ad % to labels 
pie(slices, labels = lbls,col=c("Azure2", "Aquamarine1","Darkslateblue"))
```