-
Notifications
You must be signed in to change notification settings - Fork 1
/
39indQC.Rmd
108 lines (70 loc) · 4.85 KB
/
39indQC.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
---
title: "39indQC"
author: "Briana Mittleman "
date: "2018-09-24"
output: workflowr::wflow_html
---
I will use this to look at the map stats and peak stats for the full set of 39 ind.
```{r}
library(tidyverse)
library(workflowr)
library(reshape2)
library(cowplot)
library(tximport)
```
##Map Stats
```{r}
mapstats= read.csv("../data/comb_map_stats_39ind.csv", header = T, stringsAsFactors = F)
mapstats$line=as.factor(mapstats$line)
mapstats$fraction=as.factor(mapstats$fraction)
map_melt=melt(mapstats, id.vars=c("line", "fraction"), measure.vars = c("comb_reads", "comb_mapped", "comb_prop_mapped"))
prop_mapped= map_melt %>% filter(variable=="comb_prop_mapped")
mapped_reads= map_melt %>% filter(variable=="comb_mapped")
mapplot_prop=ggplot(prop_mapped, aes(y=value, x=line, fill=fraction)) + geom_bar(stat="identity",position="dodge") + labs( title="Proportion of reads mapped") + ylab("Proportion mapped")
mapplot_mapped=ggplot(mapped_reads, aes(y=value, x=line, fill=fraction)) + geom_bar(stat="identity",position="dodge") + labs( title="Number of Mapped reads") + ylab("Mapped")
plot_grid(mapplot_prop, mapplot_mapped)
```
Plot boxplots for total vs nuclear.
```{r}
box_mapprop=ggplot(prop_mapped, aes(y=value, x=fraction, fill=fraction)) + geom_boxplot() + geom_jitter(position = position_jitter(.3)) + labs( title="Map Proportion") + ylab("Mapped Proportion")
box_map=ggplot(mapped_reads, aes(y=value, x=fraction, fill=fraction)) + geom_boxplot() + geom_jitter(position = position_jitter(.3)) + labs( title="Number of Mapped reads") + ylab("Mapped")
plot_grid(box_map, box_mapprop)
```
##Genes with multiple peaks
This is similar to the analysis I ran in dataprocfigures.Rmd. I start by overlappping the refseq genes with my peaks. With the script refseq_countdistinct.sh.
```{r}
namesPeak=c("Chr", "Start", "End", "Name", "Score", "Strand", "numPeaks")
Opeakpergene=read.table("../data/peakPerRefSeqGene/filtered_APApeaks_perRefseqGene_oppStrand.txt", stringsAsFactors = F, header = F, col.names = namesPeak) %>% mutate(onePeak=ifelse(numPeaks==1, 1, 0 )) %>% mutate(multPeaks=ifelse(numPeaks > 1, 1, 0 ))
Ogenes1peak=sum(Opeakpergene$onePeak)/nrow(Opeakpergene)
OgenesMultpeak=sum(Opeakpergene$multPeaks)/nrow(Opeakpergene)
Ogenes0peak= 1- Ogenes1peak - OgenesMultpeak
OperPeak= c(round(Ogenes0peak,digits = 3), round(Ogenes1peak,digits = 3),round(OgenesMultpeak, digits = 3))
Category=c("Zero", "One", "Multiple")
OperPeakdf=as.data.frame(cbind(Category,OperPeak))
OperPeakdf$OperPeak=as.numeric(as.character(OperPeakdf$OperPeak))
Olab1=paste("Genes =", Ogenes0peak*nrow(Opeakpergene), sep=" ")
Olab2=paste("Genes =", sum(Opeakpergene$onePeak), sep=" ")
Olab3=paste("Genes =", sum(Opeakpergene$multPeaks), sep=" ")
Ogenepeakplot=ggplot(OperPeakdf, aes(x="", y=OperPeak, by=Category, fill=Category)) + geom_bar(stat="identity")+ labs(title="Characterize Refseq genes by number of PAS- Oppstrand", y="Proportion of Protein Coding gene", x="")+ scale_fill_brewer(palette="Paired") + coord_cartesian(ylim=c(0,1)) + annotate("text", x="", y= .2, label=Olab1) + annotate("text", x="", y= .4, label=Olab2) + annotate("text", x="", y= .9, label=Olab3)
Ogenepeakplot
```
I will now repull in the RNA seq data for one of my lines to look at the expression levels of the genes with at least 1 called peak.
```{r}
tx2gene=read.table("../data/RNAkalisto/ncbiRefSeq.txn2gene.txt" ,header= F, sep="\t", stringsAsFactors = F)
txi.kallisto.tsv <- tximport("../data/RNAkalisto/abundance.tsv", type = "kallisto", tx2gene = tx2gene)
txi.kallisto.tsv$abundance= as.data.frame(txi.kallisto.tsv$abundance) %>% rownames_to_column(var="Name")
colnames(txi.kallisto.tsv$abundance)= c("Name", "TPM")
#genes with >0 TPM and at least 1 peak
refPeakandRNA_withO_TPM=Opeakpergene %>% inner_join(txi.kallisto.tsv$abundance, by="Name") %>% filter(TPM>0, numPeaks>0)
#genes with >0 TPM and 0 peak
refPeakandRNA_noPeakw_withO_TPM=Opeakpergene %>% inner_join(txi.kallisto.tsv$abundance, by="Name") %>% filter(TPM >0, numPeaks==0)
#plot
plot(sort(log10(refPeakandRNA_withO_TPM$TPM), decreasing = T), main="Distribution of RNA expression 18486", ylab="log10 TPM", xlab="Refseq Gene")
points(sort(log10(refPeakandRNA_noPeakw_withO_TPM$TPM), decreasing = T), col="Red")
legend("topright", legend=c("Genes wth Peak", "Genes without Peak"), col=c("black", "red"),pch=19)
```
Plot this as distributions.
```{r}
comp_RNAtpm=ggplot(refPeakandRNA_withO_TPM, aes(x=log10(TPM))) + geom_histogram(binwidth=.5, alpha=.5) +geom_histogram(data = refPeakandRNA_noPeakw_withO_TPM, aes(x=log10(TPM)), fill="Red", alpha=.5, binwidth=.5) + labs(title="Comparing RNA expression for genes with a PAS vs no PAS") + annotate("text", x=-8, y=3000, col="Red", label="Genes without PAS") + annotate("text", x=-8.1, y=2700, col="Black", label="Genes with PAS") + geom_rect(linetype=1, xmin=-10, xmax=-6, ymin=2500, ymax=3200, color="Black", alpha=0)
comp_RNAtpm
```