PAS descriptive plots

Peaks per gene:

I want to plot how many genes have 0, 1, 2 and more than 2 PAS in the set. I need to join my PAS with the annotation to find out how many genes have 0 PAS.

pas=read.table("../data/PAS/APAPAS_GeneLocAnno.5perc.bed", header = F, stringsAsFactors = F, col.names = c("Chr", "start", "end", "PeakID", "score", "strand")) %>% separate(PeakID, into=c("peaknum", "geneAnno"), sep=":") %>% separate(geneAnno, into=c("Gene", "Loc"),sep="_")

pasbygene= pas %>% group_by(Gene) %>% summarise(PAS=n())

annotation=read.table("../../genome_anotation_data/RefSeq_annotations/ncbiRefSeq_FormatedallAnnotation.sort.bed", col.names = c("chr", "start", "end", "anno", "score", "strand")) %>% separate(anno, into=c("Loc", "Gene"),sep=":") %>% group_by(Gene) %>% summarise(annos=n()) %>% dplyr::select(Gene)

PASallgene=annotation %>% full_join(pasbygene, by="Gene") %>% replace_na(list(PAS=0)) 

#group with 0,1,2,more than 2
PASallgene_grouped=PASallgene %>% mutate(Zero=ifelse(PAS==0,1, 0), One=ifelse(PAS==1,1,0), Multiple=ifelse(PAS>1,1,0))

Plot this:

Genes=c(sum(PASallgene_grouped$Zero),sum(PASallgene_grouped$One),sum(PASallgene_grouped$Multiple))
PAS=c("Zero", "One", "Multiple")
AllPAS=c(sum(PASallgene_grouped$Zero),sum(PASallgene_grouped$One),sum(PASallgene_grouped$Multiple))
GenebyPAS=as.data.frame(cbind(PAS,AllPAS))
GenebyPAS$AllPAS=as.numeric(as.character(GenebyPAS$AllPAS))

allPASplot=ggplot(GenebyPAS, aes(x="",y=AllPAS, fill=PAS)) + geom_bar(stat="identity", width=.5) + scale_fill_brewer(palette="GnBu") + labs(title="Identified PAS per Gene", y="Genes",x="All Indentified PAS")
allPASplot

Past versions of unnamed-chunk-3-1.png

ggsave(allPASplot, file="../output/GeneswithAPApotentialAllPAS.png", width=3, height=5)

pasUTR=pas %>% filter(Loc=="utr3") %>% group_by(Gene) %>% summarise(PAS=n())

pasUTR_allgene=annotation %>% full_join(pasUTR, by="Gene") %>% replace_na(list(PAS=0)) 

PASUTRallgene_grouped=pasUTR_allgene %>% mutate(Zero=ifelse(PAS==0,1, 0), One=ifelse(PAS==1,1,0), Multiple=ifelse(PAS>1,1,0))


GenesUTR=c(sum(PASUTRallgene_grouped$Zero),sum(PASUTRallgene_grouped$One),sum(PASUTRallgene_grouped$Multiple))
UTR=c(sum(PASUTRallgene_grouped$Zero),sum(PASUTRallgene_grouped$One),sum(PASUTRallgene_grouped$Multiple))
GenebyPASUTR=as.data.frame(cbind(PAS,UTR))
GenebyPASUTR$UTR=as.numeric(as.character(GenebyPASUTR$UTR))

ggplot(GenebyPASUTR, aes(x="",y=UTR, fill=PAS)) + geom_bar(stat="identity")

pasIntron=pas %>% filter(Loc=="intron" | Loc=='utr3') %>% group_by(Gene) %>% summarise(PAS=n())

pasIntron_allgene=annotation %>% full_join(pasIntron, by="Gene") %>% replace_na(list(PAS=0)) 

pasIntronallgene_grouped=pasIntron_allgene %>% mutate(Zero=ifelse(PAS==0,1, 0), One=ifelse(PAS==1,1,0), Multiple=ifelse(PAS>1,1,0))


UTRandIntron=c(sum(pasIntronallgene_grouped$Zero),sum(pasIntronallgene_grouped$One),sum(pasIntronallgene_grouped$Multiple))
GenebyPASIntron=as.data.frame(cbind(PAS,UTRandIntron))
GenebyPASIntron$UTRandIntron=as.numeric(as.character(GenebyPASIntron$UTRandIntron))

ggplot(GenebyPASIntron, aes(x="",y=UTRandIntron, fill=PAS)) + geom_bar(stat="identity")

GenebyPASUTR_melt=melt(GenebyPASUTR, id.vars = "PAS", value.name = "Genes", variable.name = "Set")

GenebyPAS_melt=melt(GenebyPAS, id.vars = "PAS", value.name = "Genes", variable.name = "Set")

GenebyPASIntron_melt=melt(GenebyPASIntron, id.vars = "PAS", value.name = "Genes", variable.name = "Set")

GenebyPAStoplot=rbind(GenebyPAS_melt,GenebyPASUTR_melt,GenebyPASIntron_melt)

geneswithAPA=ggplot(GenebyPAStoplot, aes(x=Set,y=Genes, fill=PAS, by=Set)) + geom_bar(stat="identity")+ scale_fill_brewer(palette="YlGnBu") + labs(title="Genes with APA poential")

geneswithAPA

ggsave(geneswithAPA, file="../output/GeneswithAPApotential.png")

Saving 7 x 5 in image

GenebyPAStoplot

       PAS          Set Genes
1     Zero       AllPAS 11659
2      One       AllPAS  4111
3 Multiple       AllPAS 11345
4     Zero          UTR 14503
5      One          UTR  6996
6 Multiple          UTR  5616
7     Zero UTRandIntron 13025
8      One UTRandIntron  4328
9 Multiple UTRandIntron  9762

Location of PAS

PAS_loc=pas %>% group_by(Loc) %>% summarise(nPAS=n())
loclabel=c("Coding", "Downstream", "Intronic", "3' UTR", "5' UTR")
PASLocPlot=ggplot(PAS_loc, aes(x=Loc, y=nPAS, fill=Loc)) + geom_bar(stat="identity",width=.5)+ scale_fill_brewer(palette = "YlGnBu") + labs(x="Gene location", y="Number of identified PAS", title="Location distribution for identified PAS") + theme(legend.position = "none")+ scale_x_discrete(labels= loclabel)+theme(axis.text.x = element_text(angle = 90, hjust = 1))
PASLocPlot

Past versions of unnamed-chunk-9-1.png

ggsave(PASLocPlot, file="../output/PASlocation.png")

Saving 7 x 5 in image

Number of genes with apa by cutoff

I want to make a script that takes a cuttoff and tells me how many gens have 0,1, >1 PAS. This way I can put these together to make stacked barplots:

I can make the plot for total then again for nuclear.

The annotaiton is annotation:

totapaanno=read.table("../data/phenotype/APApeak_Phenotype_GeneLocAnno.Total.fc",header = T,stringsAsFactors = F) 
indiv=colnames(totapaanno)[2:55]
totapanum=read.table("../data/phenotype/APApeak_Phenotype_GeneLocAnno.Total.CountsOnlyNumeric",header = F, col.names = indiv) 
totapa_mean=rowMeans(totapanum)
totapaMeananno=as.data.frame(cbind(ID=totapaanno$chrom, meanUsage=totapa_mean))
totapaMeananno$meanUsage=as.numeric(as.character(totapaMeananno$meanUsage))
totapaMeananno$ID=as.character(totapaMeananno$ID)

genesbycuttoff_tot=function(fraction){
  totapaMeananno_filt=totapaMeananno %>% filter(meanUsage >=fraction) %>% separate(ID, into=c("chrom", "start","end", "peakID"),sep=":") %>% separate(peakID, into=c("Gene","loc", "strand", "peak"),sep="_") %>% group_by(Gene) %>% summarise(PAS=n())
  PASallgene=annotation %>% full_join(totapaMeananno_filt, by="Gene") %>% replace_na(list(PAS=0))
  PASallgene_cat=PASallgene %>% mutate(Category=ifelse(PAS==0,"Zero", ifelse(PAS==1, "One", "Multiple"))) %>% group_by(Category) %>% summarise(NPer=n())
return(PASallgene_cat$NPer)
}

#multiple, one, zero 
categories=c("Multiple_PAS", "One_PAS", "Zero_PAS")
FullDF=as.data.frame(cbind(categories))
cutoffs=seq(from=0, to=.5, by=.05)
for (val in cutoffs)
{
FullDF=cbind(FullDF,val=genesbycuttoff_tot(val))
}

Warning: Expected 4 pieces. Additional pieces discarded in 4 rows [14735,
14736, 14737, 14738].

Warning: Expected 4 pieces. Additional pieces discarded in 3 rows [9347,
9348, 9349].

Warning: Expected 4 pieces. Additional pieces discarded in 3 rows [7284,
7285, 7286].

Warning: Expected 4 pieces. Additional pieces discarded in 2 rows [6143,
6144].

Warning: Expected 4 pieces. Additional pieces discarded in 2 rows [5328,
5329].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [4704].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [4169].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [3738].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [3327].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [2987].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [2685].

colnames(FullDF)=c("Category",cutoffs)

Melt:

fullDF_melt=melt(FullDF,id.vars = "Category",variable.name = "Cutoff", value.name = "NGenes") %>% mutate(propGene=NGenes/nrow(annotation))

ggplot(fullDF_melt,aes(x=Cutoff, y=propGene, by=Category, fill=Category)) + geom_bar( stat="identity",width = .5) + scale_fill_brewer(palette="GnBu") + labs(title="Proportion of genes with Multiple PAS by Total Usage filter",y="Proportion of 27115 genes", x="Usage Filter cutoff") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_x_discrete(name="Usage Filter cutoff", breaks=c("0","0.1","0.2", "0.3", "0.4","0.5"))

Past versions of unnamed-chunk-13-1.png

Nuclear:

nucapaanno=read.table("../data/phenotype/APApeak_Phenotype_GeneLocAnno.Nuclear.fc",header = T,stringsAsFactors = F) 

nucapanum=read.table("../data/phenotype/APApeak_Phenotype_GeneLocAnno.Nuclear.CountsOnlyNumeric",header = F, col.names = indiv) 
nucapa_mean=rowMeans(nucapanum)
nucapaMeananno=as.data.frame(cbind(ID=nucapaanno$chrom, meanUsage=nucapa_mean))
nucapaMeananno$meanUsage=as.numeric(as.character(nucapaMeananno$meanUsage))
nucapaMeananno$ID=as.character(nucapaMeananno$ID)

genesbycuttoff_nuc=function(fraction){
  nucapaMeananno_filt=nucapaMeananno %>% filter(meanUsage >=fraction) %>% separate(ID, into=c("chrom", "start","end", "peakID"),sep=":") %>% separate(peakID, into=c("Gene","loc", "strand", "peak"),sep="_") %>% group_by(Gene) %>% summarise(PAS=n())
  PASallgene=annotation %>% full_join(nucapaMeananno_filt, by="Gene") %>% replace_na(list(PAS=0))
  PASallgene_cat=PASallgene %>% mutate(Category=ifelse(PAS==0,"Zero", ifelse(PAS==1, "One", "Multiple"))) %>% group_by(Category) %>% summarise(NPer=n())
return(PASallgene_cat$NPer)
}

#multiple, one, zero 
FullDFNuc=as.data.frame(cbind(categories))
for (val in cutoffs)
{
FullDFNuc=cbind(FullDFNuc,val=genesbycuttoff_nuc(val))
}

Warning: Expected 4 pieces. Additional pieces discarded in 4 rows [14735,
14736, 14737, 14738].

Warning: Expected 4 pieces. Additional pieces discarded in 4 rows [9492,
9493, 9494, 9495].

Warning: Expected 4 pieces. Additional pieces discarded in 3 rows [7512,
7513, 7514].

Warning: Expected 4 pieces. Additional pieces discarded in 2 rows [6329,
6330].

Warning: Expected 4 pieces. Additional pieces discarded in 2 rows [5636,
5637].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [5028].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [4561].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [4146].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [3782].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [3457].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [3178].

colnames(FullDFNuc)=c("Category",cutoffs)

FullDFNuc_melt=melt(FullDFNuc,id.vars = "Category",variable.name = "Cutoff", value.name = "NGenes") %>% mutate(propGene=NGenes/nrow(annotation))

ggplot(FullDFNuc_melt,aes(x=Cutoff, y=propGene, by=Category, fill=Category)) + geom_bar( stat="identity") + scale_fill_brewer(palette="GnBu") + labs(title="Proportion of genes with Multiple PAS by Nuclear Usage filter",y="Proportion of 27115 genes", x="Usage Filter cutoff")

Past versions of unnamed-chunk-17-1.png

Location of PAS by filter

I will make these plots but the categories will be location of the PAS.

locbycutoff_tot=function(fraction){
  totapaMeananno_filt=totapaMeananno %>% filter(meanUsage >=fraction) %>% separate(ID, into=c("chrom", "start","end", "peakID"),sep=":") %>% separate(peakID, into=c("Gene","loc", "strand", "peak"),sep="_") %>% group_by(loc) %>% summarise(PerLoc=n()) %>%filter(loc!= "008559")
return(totapaMeananno_filt$PerLoc)
}

locations=c("cds", "end", "intron", "utr3", "utr5")
FullDF_loc=as.data.frame(cbind(locations))
cutoffs=seq(from=0, to=.5, by=.05)
for (val in cutoffs)
{
FullDF_loc=cbind(FullDF_loc,val=locbycutoff_tot(val))
}

Warning: Expected 4 pieces. Additional pieces discarded in 4 rows [14735,
14736, 14737, 14738].

Warning: Expected 4 pieces. Additional pieces discarded in 3 rows [9347,
9348, 9349].

Warning: Expected 4 pieces. Additional pieces discarded in 3 rows [7284,
7285, 7286].

Warning: Expected 4 pieces. Additional pieces discarded in 2 rows [6143,
6144].

Warning: Expected 4 pieces. Additional pieces discarded in 2 rows [5328,
5329].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [4704].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [4169].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [3738].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [3327].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [2987].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [2685].

colnames(FullDF_loc)=c("Location",cutoffs)

Melt:

FullDF_loc_melt=melt(FullDF_loc,id.vars = "Location",variable.name = "Cutoff", value.name = "NPas") %>% group_by(Cutoff) %>%  mutate(propPAS=NPas/sum(NPas))

totplotloc=ggplot(FullDF_loc_melt,aes(x=Cutoff, y=propPAS, by=Location, fill=Location)) + geom_bar(width=.5, stat="identity") + scale_fill_brewer(palette="GnBu") + labs(title="PAS location\n Total Fraction",y="Proportion of PAS", x="Usage Filter cutoff")+ theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_x_discrete(name="Usage Filter cutoff", breaks=c("0","0.1","0.2", "0.3", "0.4","0.5"))

Nuclear

locbycutoff_nuc=function(fraction){
nucapaMeananno_filt=nucapaMeananno %>% filter(meanUsage >=fraction) %>% separate(ID, into=c("chrom", "start","end", "peakID"),sep=":") %>% separate(peakID, into=c("Gene","loc", "strand", "peak"),sep="_") %>% group_by(loc) %>% summarise(PerLoc=n()) %>%filter(loc!= "008559")
return(nucapaMeananno_filt$PerLoc)
}

NucFullDF_loc=as.data.frame(cbind(locations))
cutoffs=seq(from=0, to=.5, by=.05)
for (val in cutoffs)
{
NucFullDF_loc=cbind(NucFullDF_loc,val=locbycutoff_nuc(val))
}

Warning: Expected 4 pieces. Additional pieces discarded in 4 rows [14735,
14736, 14737, 14738].

Warning: Expected 4 pieces. Additional pieces discarded in 4 rows [9492,
9493, 9494, 9495].

Warning: Expected 4 pieces. Additional pieces discarded in 3 rows [7512,
7513, 7514].

Warning: Expected 4 pieces. Additional pieces discarded in 2 rows [6329,
6330].

Warning: Expected 4 pieces. Additional pieces discarded in 2 rows [5636,
5637].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [5028].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [4561].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [4146].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [3782].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [3457].

Warning: Expected 4 pieces. Additional pieces discarded in 1 rows [3178].

colnames(NucFullDF_loc)=c("Location",cutoffs)

Melt:

NucFullDF_locMelt=melt(NucFullDF_loc,id.vars = "Location",variable.name = "Cutoff", value.name = "NPas") %>% group_by(Cutoff) %>%  mutate(propPAS=NPas/sum(NPas))

nucplotloc=ggplot(NucFullDF_locMelt,aes(x=Cutoff, y=propPAS, by=Location, fill=Location)) + geom_bar(width = .5, stat="identity") + scale_fill_brewer(palette="GnBu") + labs(title="PAS location\n Nuclear Fraction",y="Proportion of PAS", x="Usage Filter cutoff")+ theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_x_discrete(name="Usage Filter cutoff", breaks=c("0","0.1","0.2", "0.3", "0.4","0.5"))

Plot next to eachother

plot_grid(nucplotloc, totplotloc)

Past versions of unnamed-chunk-24-1.png

totplotloc

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] cowplot_0.9.4   reshape2_1.4.3  forcats_0.3.0   stringr_1.3.1  
 [5] dplyr_0.8.0.1   purrr_0.3.2     readr_1.3.1     tidyr_0.8.3    
 [9] tibble_2.1.1    ggplot2_3.1.1   tidyverse_1.2.1 workflowr_1.3.0

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.0         RColorBrewer_1.1-2 cellranger_1.1.0  
 [4] pillar_1.3.1       compiler_3.5.1     git2r_0.23.0      
 [7] plyr_1.8.4         tools_3.5.1        digest_0.6.18     
[10] lubridate_1.7.4    jsonlite_1.6       evaluate_0.12     
[13] nlme_3.1-137       gtable_0.2.0       lattice_0.20-38   
[16] pkgconfig_2.0.2    rlang_0.3.1        cli_1.0.1         
[19] rstudioapi_0.10    yaml_2.2.0         haven_1.1.2       
[22] withr_2.1.2        xml2_1.2.0         httr_1.3.1        
[25] knitr_1.20         hms_0.4.2          generics_0.0.2    
[28] fs_1.2.6           rprojroot_1.3-2    grid_3.5.1        
[31] tidyselect_0.2.5   glue_1.3.0         R6_2.3.0          
[34] readxl_1.1.0       rmarkdown_1.10     modelr_0.1.2      
[37] magrittr_1.5       whisker_0.3-2      backports_1.1.2   
[40] scales_1.0.0       htmltools_0.3.6    rvest_0.3.2       
[43] assertthat_0.2.0   colorspace_1.3-2   labeling_0.3      
[46] stringi_1.2.4      lazyeval_0.2.1     munsell_0.5.0     
[49] broom_0.5.1        crayon_1.3.4