Finish the AC-16 analysis.

RIP-seq.Rmd

@@ -1059,3 +1059,300 @@ alk_allsnotx_rna_rpkm$gene_name <- as.character(gencode_sn_mapping[sn_inx,"V3"])
 write.table(alk_allsnotx_rna_rpkm, file="ac16-snRNA_RPKM_RNA.tsv", quote=F, sep="\t", row.names=F, col.names=T)
 
 ```
+
+```{r alignments TNF}
+ALKT_RNA1 <- "/Volumes/project/GCRB/Lee_Lab/s163035/DK_RIP-seq/PARP1_RNAseq_in_AC16_MCF7_withTNFa_E2/raw/Sample_ALKT_RNA1/align-star-se.sh-1.0.0/ALKT_RNA1_GCCAAT_L007_R1.rDNA.filtered.fastq.gz_filtered_chrMAligned.sortedByCoord.out.bam"
+ALKT_RNA2 <- "/Volumes/project/GCRB/Lee_Lab/s163035/DK_RIP-seq/PARP1_RNAseq_in_AC16_MCF7_withTNFa_E2/raw/Sample_ALKT_RNA2/align-star-se.sh-1.0.0/ALKT_RNA2_TAGCTT_L007_R1.rDNA.filtered.fastq.gz_filtered_chrMAligned.sortedByCoord.out.bam"
+
+ALKT_RIP_PI1 <- "/Volumes/project/GCRB/Lee_Lab/s163035/DK_RIP-seq/PARP1_RIPseq_in_AC16_MCF7_withTNFa_E2/raw/Sample_ALKT_RIP1_PI/align-star-se.sh-2.0.0/ALKT_RIP1_PI_ATGTCA_L003_R1.rDNA.filtered.fastq.gz_filtered_chrMAligned.sortedByCoord.out.bam"
+ALKT_RIP_PI2 <- "/Volumes/project/GCRB/Lee_Lab/s163035/DK_RIP-seq/PARP1_RIPseq_in_AC16_MCF7_withTNFa_E2/raw/Sample_ALKT_RIP2_PI/align-star-se.sh-2.0.0/ALKT_RIP2_PI_GAGTGG_L003_R1.rDNA.filtered.fastq.gz_filtered_chrMAligned.sortedByCoord.out.bam"
+
+ALKT_RIP_PARP1 <- "/Volumes/project/GCRB/Lee_Lab/s163035/DK_RIP-seq/PARP1_RIPseq_in_AC16_MCF7_withTNFa_E2/raw/Sample_ALKT_RIP1_PARP1/align-star-se.sh-2.0.0/ALKT_RIP1_PARP1_CCGTCC_L006_R1.rDNA.filtered.fastq.gz_filtered_chrMAligned.sortedByCoord.out.bam"
+ALKT_RIP_PARP2 <- "/Volumes/project/GCRB/Lee_Lab/s163035/DK_RIP-seq/PARP1_RIPseq_in_AC16_MCF7_withTNFa_E2/raw/Sample_ALKT_RIP2_PARP1/align-star-se.sh-2.0.0/ALKT_RIP2_PARP1_GGTAGC_L006_R1.rDNA.filtered.fastq.gz_filtered_chrMAligned.sortedByCoord.out.bam"
+```
+
+```{r counts TNF}
+ALKT_RNA1_counts_ref <- CountRefseq(ref_tx, ALKT_RNA1, all_sn_tx, removeSNORNA=TRUE)
+ALKT_RNA2_counts_ref <- CountRefseq(ref_tx, ALKT_RNA2, all_sn_tx, removeSNORNA=TRUE)
+
+ALKT_RNA1_counts_lnc <- CountRefseq(lnc_tx, ALKT_RNA1, all_sn_tx, removeSNORNA=TRUE)
+ALKT_RNA2_counts_lnc <- CountRefseq(lnc_tx, ALKT_RNA2, all_sn_tx, removeSNORNA=TRUE)
+
+ALKT_RNA1_counts_sn <- CountSNRNA(sn_tx, ALKT_RNA1)
+ALKT_RNA2_counts_sn <- CountSNRNA(sn_tx, ALKT_RNA2)
+
+ALKT_RNA1_counts_sno <- CountSNRNA(sno_tx, ALKT_RNA1)
+ALKT_RNA2_counts_sno <- CountSNRNA(sno_tx, ALKT_RNA2)
+
+
+ALKT_RIP_PI1_counts_ref <- CountRefseq(ref_tx, ALKT_RIP_PI1, all_sn_tx, removeSNORNA=TRUE)
+ALKT_RIP_PI2_counts_ref <- CountRefseq(ref_tx, ALKT_RIP_PI2, all_sn_tx, removeSNORNA=TRUE)
+
+ALKT_RIP_PI1_counts_lnc <- CountRefseq(lnc_tx, ALKT_RIP_PI1, all_sn_tx, removeSNORNA=TRUE)
+ALKT_RIP_PI2_counts_lnc <- CountRefseq(lnc_tx, ALKT_RIP_PI2, all_sn_tx, removeSNORNA=TRUE)
+
+ALKT_RIP_PI1_counts_sn <- CountSNRNA(sn_tx, ALKT_RIP_PI1)
+ALKT_RIP_PI2_counts_sn <- CountSNRNA(sn_tx, ALKT_RIP_PI2)
+
+ALKT_RIP_PI1_counts_sno <- CountSNRNA(sno_tx, ALKT_RIP_PI1)
+ALKT_RIP_PI2_counts_sno <- CountSNRNA(sno_tx, ALKT_RIP_PI2)
+
+ALKT_RIP_PARP1_counts_ref <- CountRefseq(ref_tx, ALKT_RIP_PARP1, all_sn_tx, removeSNORNA=TRUE)
+ALKT_RIP_PARP2_counts_ref <- CountRefseq(ref_tx, ALKT_RIP_PARP2, all_sn_tx, removeSNORNA=TRUE)
+
+ALKT_RIP_PARP1_counts_lnc <- CountRefseq(lnc_tx, ALKT_RIP_PARP1, all_sn_tx, removeSNORNA=TRUE)
+ALKT_RIP_PARP2_counts_lnc <- CountRefseq(lnc_tx, ALKT_RIP_PARP2, all_sn_tx, removeSNORNA=TRUE)
+
+ALKT_RIP_PARP1_counts_sn <- CountSNRNA(sn_tx, ALKT_RIP_PARP1)
+ALKT_RIP_PARP2_counts_sn <- CountSNRNA(sn_tx, ALKT_RIP_PARP2)
+
+ALKT_RIP_PARP1_counts_sno <- CountSNRNA(sno_tx, ALKT_RIP_PARP1)
+ALKT_RIP_PARP2_counts_sno <- CountSNRNA(sno_tx, ALKT_RIP_PARP2)
+
+
+```
+
+# RPKM Cutoff
+``` {r rpkm TNF}
+ALKT_RNA_counts_ref <- ALKT_RNA1_counts_ref$counts +  ALKT_RNA2_counts_ref$counts
+ALKT_RNA_lib_ref <- ALKT_RNA1_counts_ref$lib +  ALKT_RNA2_counts_ref$lib
+ALKT_RNA_ref_d <- DGEList(counts=as.matrix(ALKT_RNA_counts_ref), lib.size=ALKT_RNA_lib_ref)
+ALKT_RNA_ref_rpkm <- data.frame(rpkm(ALKT_RNA_ref_d,width(ref_tx)))
+alkt_pc_rpkm <- ref_tx[ALKT_RNA_ref_rpkm>=1]
+NROW(unique(unlist(mcols(alkt_pc_rpkm)$tx_name)))
+NROW(unique(unlist(mcols(alkt_pc_rpkm)$gene_name)))
+
+ALKT_RIP_PARP_counts_ref <- ALKT_RIP_PARP1_counts_ref$counts +  ALKT_RIP_PARP2_counts_ref$counts
+ALKT_RIP_PARP_lib_ref <- ALKT_RIP_PARP1_counts_ref$lib +  ALKT_RIP_PARP2_counts_ref$lib
+ALKT_RIP_PI_counts_ref <- ALKT_RIP_PI1_counts_ref$counts +  ALKT_RIP_PI2_counts_ref$counts
+ALKT_RIP_PI_lib_ref <- ALKT_RIP_PI1_counts_ref$lib +  ALKT_RIP_PI1_counts_ref$lib
+ALKT_RIP_PI_d <- DGEList(counts=as.matrix(ALKT_RIP_PI_counts_ref), lib.size=ALKT_RIP_PI_lib_ref)
+ALKT_RIP_PI_rpkm <- data.frame(rpkm(ALKT_RIP_PI_d,width(ref_tx)))
+ALKT_RIP_PARP_d <- DGEList(counts=as.matrix(ALKT_RIP_PARP_counts_ref), lib.size=ALKT_RIP_PARP_lib_ref)
+ALKT_RIP_PARP_rpkm <- data.frame(rpkm(ALKT_RIP_PARP_d,width(ref_tx)))
+
+ALKT_RIP_FC <- data.frame((ALKT_RIP_PARP_rpkm+1)/(ALKT_RIP_PI_rpkm+1))
+alkt_pc_rpkm_fc <- ref_tx[ALKT_RNA_ref_rpkm>=1 & ALKT_RIP_FC>=2 ]
+NROW(unique(unlist(mcols(alkt_pc_rpkm_fc)$tx_name)))
+NROW(unique(unlist(mcols(alkt_pc_rpkm_fc)$gene_name)))
+pc_ratio <- NROW(unique(unlist(mcols(alkt_pc_rpkm_fc)$gene_name)))/NROW(unique(unlist(mcols(alkt_pc_rpkm)$gene_name)))
+
+ALKT_RNA_counts_lnc <- ALKT_RNA1_counts_lnc$counts +  ALKT_RNA2_counts_lnc$counts
+ALKT_RNA_lib_lnc <- ALKT_RNA1_counts_lnc$lib +  ALKT_RNA2_counts_lnc$lib
+ALKT_RNA_lnc_d <- DGEList(counts=as.matrix(ALKT_RNA_counts_lnc), lib.size=ALKT_RNA_lib_lnc)
+ALKT_RNA_lnc_rpkm <- data.frame(rpkm(ALKT_RNA_lnc_d,width(lnc_tx)))
+alkt_lnc_rpkm <- lnc_tx[ALKT_RNA_lnc_rpkm>=1]
+NROW(unique(unlist(mcols(alkt_lnc_rpkm)$tx_name)))
+NROW(unique(unlist(mcols(alkt_lnc_rpkm)$gene_id)))
+
+ALKT_RIP_PARP_counts_lnc <- ALKT_RIP_PARP1_counts_lnc$counts +  ALKT_RIP_PARP2_counts_lnc$counts
+ALKT_RIP_PARP_lib_lnc<- ALKT_RIP_PARP1_counts_lnc$lib +  ALKT_RIP_PARP2_counts_lnc$lib
+ALKT_RIP_PI_counts_lnc <- ALKT_RIP_PI1_counts_lnc$counts +  ALKT_RIP_PI2_counts_lnc$counts
+ALKT_RIP_PI_lib_lnc <- ALKT_RIP_PI1_counts_lnc$lib +  ALKT_RIP_PI1_counts_lnc$lib
+ALKT_RIP_PI_d <- DGEList(counts=as.matrix(ALKT_RIP_PI_counts_lnc), lib.size=ALKT_RIP_PI_lib_lnc)
+ALKT_RIP_PI_rpkm <- data.frame(rpkm(ALKT_RIP_PI_d,width(lnc_tx)))
+ALKT_RIP_PARP_d <- DGEList(counts=as.matrix(ALKT_RIP_PARP_counts_lnc), lib.size=ALKT_RIP_PARP_lib_lnc)
+ALKT_RIP_PARP_rpkm <- data.frame(rpkm(ALKT_RIP_PARP_d,width(lnc_tx)))
+
+ALKT_RIP_FC <- data.frame((ALKT_RIP_PARP_rpkm+1)/(ALKT_RIP_PI_rpkm+1))
+alkt_lnc_rpkm_fc <- lnc_tx[ALKT_RNA_lnc_rpkm>=1 & ALKT_RIP_FC>=2 ]
+NROW(unique(unlist(mcols(alkt_lnc_rpkm_fc)$tx_name)))
+NROW(unique(unlist(mcols(alkt_lnc_rpkm_fc)$gene_id)))
+lnc_ratio <- NROW(unique(unlist(mcols(alkt_lnc_rpkm_fc)$gene_id)))/NROW(unique(unlist(mcols(alkt_lnc_rpkm)$gene_id)))
+
+
+ALKT_RNA_counts_sno <- ALKT_RNA1_counts_sno$counts +  ALKT_RNA2_counts_sno$counts
+ALKT_RNA_lib_sno <- ALKT_RNA1_counts_sno$lib +  ALKT_RNA1_counts_sno$lib
+ALKT_RNA_sno_d <- DGEList(counts=as.matrix(ALKT_RNA_counts_sno), lib.size=ALKT_RNA_lib_sno)
+ALKT_RNA_sno_rpkm <- data.frame(rpkm(ALKT_RNA_sno_d,width(sno_tx)))
+alkt_sno_rpkm <- sno_tx[ALKT_RNA_sno_rpkm>=1]
+NROW(unique(unlist(mcols(alkt_sno_rpkm)$tx_name)))
+
+ALKT_RIP_PARP_counts_sno <- ALKT_RIP_PARP1_counts_sno$counts +  ALKT_RIP_PARP2_counts_sno$counts
+ALKT_RIP_PARP_lib_sno<- ALKT_RIP_PARP1_counts_sno$lib +  ALKT_RIP_PARP2_counts_sno$lib
+ALKT_RIP_PI_counts_sno <- ALKT_RIP_PI1_counts_sno$counts +  ALKT_RIP_PI2_counts_sno$counts
+ALKT_RIP_PI_lib_sno <- ALKT_RIP_PI1_counts_sno$lib +  ALKT_RIP_PI1_counts_sno$lib
+ALKT_RIP_PI_d <- DGEList(counts=as.matrix(ALKT_RIP_PI_counts_sno), lib.size=ALKT_RIP_PI_lib_sno)
+ALKT_RIP_PI_rpkm <- data.frame(rpkm(ALKT_RIP_PI_d,width(sno_tx)))
+ALKT_RIP_PARP_d <- DGEList(counts=as.matrix(ALKT_RIP_PARP_counts_sno), lib.size=ALKT_RIP_PARP_lib_sno)
+ALKT_RIP_PARP_rpkm <- data.frame(rpkm(ALKT_RIP_PARP_d,width(sno_tx)))
+
+
+ALKT_RIP_FC <- data.frame((ALKT_RIP_PARP_rpkm+1)/(ALKT_RIP_PI_rpkm+1))
+alkt_sno_rpkm_fc <- sno_tx[ALKT_RNA_sno_rpkm>=1 & ALKT_RIP_FC>=2 ]
+NROW(unique(unlist(mcols(alkt_sno_rpkm_fc)$tx_name)))
+sno_ratio <- NROW(unique(unlist(mcols(alkt_sno_rpkm_fc)$tx_name)))/NROW(unique(unlist(mcols(alkt_sno_rpkm)$tx_name)))
+
+ALKT_RIP_sno <- data.frame(ALKT_RIP_PI_rpkm,ALKT_RIP_PARP_rpkm)
+colnames(ALKT_RIP_sno) <- c('PI', 'PARP1')
+ALKT_RIP_sno_rpkm <- ALKT_RIP_sno[ALKT_RNA_sno_rpkm>=1 & ALKT_RIP_FC>=2,]
+
+
+ALKT_RNA_counts_sn <- ALKT_RNA1_counts_sn$counts +  ALKT_RNA2_counts_sn$counts
+ALKT_RNA_lib_sn <- ALKT_RNA1_counts_sn$lib +  ALKT_RNA1_counts_sn$lib
+ALKT_RNA_sn_d <- DGEList(counts=as.matrix(ALKT_RNA_counts_sn), lib.size=ALKT_RNA_lib_sn)
+ALKT_RNA_sn_rpkm <- data.frame(rpkm(ALKT_RNA_sn_d,width(sn_tx)))
+alkt_sn_rpkm <- sn_tx[ALKT_RNA_sn_rpkm>=1]
+NROW(unique(unlist(mcols(alkt_sn_rpkm)$tx_name)))
+
+ALKT_RIP_PARP_counts_sn <- ALKT_RIP_PARP1_counts_sn$counts +  ALKT_RIP_PARP2_counts_sn$counts
+ALKT_RIP_PARP_lib_sn<- ALKT_RIP_PARP1_counts_sn$lib +  ALKT_RIP_PARP2_counts_sn$lib
+ALKT_RIP_PI_counts_sn <- ALKT_RIP_PI1_counts_sn$counts +  ALKT_RIP_PI2_counts_sn$counts
+ALKT_RIP_PI_lib_sn <- ALKT_RIP_PI1_counts_sn$lib +  ALKT_RIP_PI1_counts_sn$lib
+ALKT_RIP_PI_d <- DGEList(counts=as.matrix(ALKT_RIP_PI_counts_sn), lib.size=ALKT_RIP_PI_lib_sn)
+ALKT_RIP_PI_rpkm <- data.frame(rpkm(ALKT_RIP_PI_d,width(sn_tx)))
+ALKT_RIP_PARP_d <- DGEList(counts=as.matrix(ALKT_RIP_PARP_counts_sn), lib.size=ALKT_RIP_PARP_lib_sn)
+ALKT_RIP_PARP_rpkm <- data.frame(rpkm(ALKT_RIP_PARP_d,width(sn_tx)))
+
+ALKT_RIP_FC <- data.frame((ALKT_RIP_PARP_rpkm+1)/(ALKT_RIP_PI_rpkm+1))
+alkt_sn_rpkm_fc <- sn_tx[ALKT_RNA_sn_rpkm>=1 & ALKT_RIP_FC>=2 ]
+NROW(unique(unlist(mcols(alkt_sn_rpkm_fc)$tx_name)))
+
+sn_ratio <- NROW(unique(unlist(mcols(alkt_sn_rpkm_fc)$tx_name)))/NROW(unique(unlist(mcols(alkt_sn_rpkm)$tx_name)))
+
+ALKT_RIP_sn <- data.frame(ALKT_RIP_PI_rpkm,ALKT_RIP_PARP_rpkm)
+colnames(ALKT_RIP_sn) <- c('PI', 'PARP1')
+ALKT_RIP_sn_rpkm <- ALKT_RIP_sn[ALKT_RNA_sn_rpkm>=1 & ALKT_RIP_FC>=2,]
+
+alkt_allsno <- c(alkt_sno_rpkm_fc,alkt_sn_rpkm_fc)
+alkt_allsno_rpkm <- rbind(ALKT_RIP_sno_rpkm,ALKT_RIP_sn_rpkm)
+allsno_rna <- c(sno_tx,sn_tx)
+alkt_allsno_rna_rpkm <- rbind(ALKT_RNA_sno_rpkm,ALKT_RNA_sn_rpkm)
+
+# bar plot figure (B)
+# Grouped Bar Plot
+jpeg('barchart-ac16-tnf-species.jpg')
+ratios <- c(pc_ratio,lnc_ratio,sno_ratio,sn_ratio)
+barplot(ratios, main="SNORNA AC16 TNF", col=c("red","black", "dimgrey", "lightgrey"),ylim=c(0,.5))
+dev.off()
+
+boundTesting <-matrix(c(1644,5047,225,249),nrow = 2,dimnames = list(Guess = c("bound", "unbound"),Truth = c("mRNA", "snoRNA")))
+fisher.test(boundTesting)
+
+# Boxplot of Preimmune vs PARP-1 (C)
+jpeg('boxplot-ac16-tnf-species.jpg')
+boxplot(ALKT_RIP_sno_rpkm,col=(c("orange","blue")),outline=FALSE,ylim = c(0,400),yaxt="n", cex.axis=1,las=2,lwd=4,lty=1)
+axis(2, at=seq(0,400,100))
+dev.off()
+wilcox.test(ALKT_RIP_sno_rpkm[,1], ALKT_RIP_sno_rpkm[,2],paired=T)
+
+
+# Get Categories of  snoRNA
+alkt_gr_snotx <- data.frame(seqnames=seqnames(alkt_sno_rpkm),
+  starts=start(alkt_sno_rpkm)-1,
+  ends=end(alkt_sno_rpkm),
+  strand=strand(alkt_sno_rpkm),
+  tx_name = elementMetadata(alkt_sno_rpkm)$tx_name,
+  gene_id = unlist(elementMetadata(alkt_sno_rpkm)$gene_id))
+
+sno_annotations <- read.table('RNA_small_nucleolar.txt',header = TRUE, sep = "\t")
+sno_inx <- match(matrix(unlist(strsplit(as.character(alkt_gr_snotx$gene_id),split='.',fixed=T)),ncol=2,byrow=T)[,1], sno_annotations$ensembl_gene_id)
+alkt_gr_snotx$symbol <- as.character(sno_annotations[sno_inx,"symbol"])
+alkt_gr_snotx$gene_family <- as.character(sno_annotations[sno_inx,"gene_family"])
+alkt_gr_snotx$gene_family[alkt_gr_snotx$gene_family == "Small nucleolar RNAs, H/ACA box"] <- "H/ACA box"
+alkt_gr_snotx$gene_family[alkt_gr_snotx$gene_family == "Small nucleolar RNAs, C/D box"] <- "C/D box"
+alkt_gr_snotx$gene_family[alkt_gr_snotx$gene_family == "Small Cajal body-specific RNAs"] <- "scaRNA"
+alkt_gr_snotx$gene_family[alkt_gr_snotx$gene_family == ""] <- "Other"
+alkt_gr_snotx$gene_family[is.na(alkt_gr_snotx$gene_family)] <- "Other"
+
+ha_rna <- dim(alkt_gr_snotx[alkt_gr_snotx$gene_family == "H/ACA box",])[1]
+cd_rna <- dim(alkt_gr_snotx[alkt_gr_snotx$gene_family == "C/D box",])[1]
+sc_rna <- dim(alkt_gr_snotx[alkt_gr_snotx$gene_family == "scaRNA",])[1]
+
+dat = data.frame(count=c(ha_rna, cd_rna, sc_rna), category=c("H/ACA box", "C/D box", "scaRNA"))
+dat$fraction = dat$count / sum(dat$count)
+
+dat$ymax = cumsum(dat$fraction)
+dat$ymin = c(0, head(dat$ymax, n=-1))
+
+dat$category <- factor(dat$category, levels = c("H/ACA box", "C/D box", "scaRNA"))
+
+# (E)
+jpeg('piechart-ac16-tnf-species.jpg')
+ggplot(dat, aes(fill=category, ymax=ymax, ymin=ymin, xmax=4, xmin=2)) +
+  geom_rect(fill=c('darkorchid1', "darkorange1","darkgreen")) +
+   coord_polar(theta="y") +
+   xlim(c(0, 4)) +
+   theme_bw() +
+   theme(panel.grid=element_blank()) +
+   theme(axis.text=element_blank()) +
+   theme(axis.ticks=element_blank())
+dev.off()
+
+
+alkt_gr_snotx <- data.frame(seqnames=seqnames(alkt_sno_rpkm_fc),
+  starts=start(alkt_sno_rpkm_fc)-1,
+  ends=end(alkt_sno_rpkm_fc),
+  strand=strand(alkt_sno_rpkm_fc),
+  tx_name = elementMetadata(alkt_sno_rpkm_fc)$tx_name,
+  gene_id = unlist(elementMetadata(alkt_sno_rpkm_fc)$gene_id))
+
+sno_inx <- match(matrix(unlist(strsplit(as.character(alkt_gr_snotx$gene_id),split='.',fixed=T)),ncol=2,byrow=T)[,1], sno_annotations$ensembl_gene_id)
+alkt_gr_snotx$symbol <- as.character(sno_annotations[sno_inx,"symbol"])
+alkt_gr_snotx$gene_family <- as.character(sno_annotations[sno_inx,"gene_family"])
+alkt_gr_snotx$gene_family[alkt_gr_snotx$gene_family == "Small nucleolar RNAs, H/ACA box"] <- "H/ACA box"
+alkt_gr_snotx$gene_family[alkt_gr_snotx$gene_family == "Small nucleolar RNAs, C/D box"] <- "C/D box"
+alkt_gr_snotx$gene_family[alkt_gr_snotx$gene_family == "Small Cajal body-specific RNAs"] <- "scaRNA"
+alkt_gr_snotx$gene_family[alkt_gr_snotx$gene_family == ""] <- "Other"
+alkt_gr_snotx$gene_family[is.na(alkt_gr_snotx$gene_family)] <- "Other"
+
+ha_rip <- dim(alkt_gr_snotx[alkt_gr_snotx$gene_family == "H/ACA box",])[1]
+cd_rip <- dim(alkt_gr_snotx[alkt_gr_snotx$gene_family == "C/D box",])[1]
+sc_rip<- dim(alkt_gr_snotx[alkt_gr_snotx$gene_family == "scaRNA",])[1]
+
+# bar plot figure (F)
+# Bar plot
+jpeg('barchart-ac16-tnf-species_snospecies.jpg')
+ratios <- c(sn_ratio,ha_rip/ha_rna,cd_rip/cd_rna,sc_rip/sc_rna)
+barplot(ratios, main="snoRNA species AC16 TNF", col=c('pink1','darkorchid1', "darkorange1","darkgreen"),ylim=c(0,1))
+dev.off()
+
+typeTesting <-matrix(c(68,14,61,98),nrow = 2,dimnames = list(Guess = c("bound", "unbound"),Truth = c("HA", "CD")))
+fisher.test(typeTesting)
+
+# Combine snRNA and snoRNA and print out
+allsnotx_tnf <- data.frame(seqnames=seqnames(alkt_allsno),
+  starts=start(alkt_allsno)-1,
+  ends=end(alkt_allsno),
+  strand=strand(alkt_allsno),
+  tx_id = elementMetadata(alkt_allsno)$tx_name,
+  gene_id = unlist(elementMetadata(alkt_allsno)$gene_id))
+
+allsnotx_tnf_rpkm <- cbind(allsnotx_tnf,alkt_allsno_rpkm)
+sno_mapping <- read.csv("./gencode.v19.annotation_snoRNA_mapping.txt", header=F, sep=',')
+sn_mapping <- read.csv("./gencode.v19.annotation_snRNA_mapping.txt", header=F, sep=',')
+gencode_sn_mapping <-rbind(sno_mapping,sn_mapping)
+sn_inx <- match(allsnotx_tnf_rpkm$tx_id, gencode_sn_mapping$V2)
+allsnotx_tnf_rpkm$tx_name <- as.character(gencode_sn_mapping[sn_inx,"V4"])
+allsnotx_tnf_rpkm$gene_name <- as.character(gencode_sn_mapping[sn_inx,"V3"])
+write.table(allsnotx_tnf_rpkm, file="ac16-tnf-snRNA_RPKM_1_FC_2.tsv", quote=F, sep="\t", row.names=F, col.names=T)
+
+allsnotx_rna <- data.frame(seqnames=seqnames(allsno_rna),
+  starts=start(allsno_rna)-1,
+  ends=end(allsno_rna),
+  strand=strand(allsno_rna),
+  tx_id = elementMetadata(allsno_rna)$tx_name,
+  gene_id = unlist(elementMetadata(allsno_rna)$gene_id))
+
+alkt_allsnotx_rna_rpkm <- cbind(allsnotx_rna,alkt_allsno_rna_rpkm)
+sno_mapping <- read.csv("./gencode.v19.annotation_snoRNA_mapping.txt", header=F, sep=',')
+sn_mapping <- read.csv("./gencode.v19.annotation_snRNA_mapping.txt", header=F, sep=',')
+gencode_sn_mapping <-rbind(sno_mapping,sn_mapping)
+sn_inx <- match(alkt_allsnotx_rna_rpkm$tx_id, gencode_sn_mapping$V2)
+alkt_allsnotx_rna_rpkm$tx_name <- as.character(gencode_sn_mapping[sn_inx,"V4"])
+alkt_allsnotx_rna_rpkm$gene_name <- as.character(gencode_sn_mapping[sn_inx,"V3"])
+write.table(alkt_allsnotx_rna_rpkm, file="ac16-tnf-snRNA_RPKM_RNA.tsv", quote=F, sep="\t", row.names=F, col.names=T)
+
+```
+
+# Overlap of snoRNA and snRNA
+``` {r snoRNA Overlap}
+
+# Overlap (A)
+ac16_unique <- unique(alk_gr_snotx$tx_name)
+tnf_unique <- unique(alkt_gr_snotx$tx_name)
+overlap <- calculate.overlap(x = list("Basal" = ac16_unique,"TNF" = tnf_unique))
+jpeg('venndigram_basal_tnf_snoRNA.jpg')
+draw.pairwise.venn(area1 = length(overlap$a1), area2 = length(overlap$a2), cross.area = length(overlap$a3),col = c("royalblue4", "seagreen4"),lwd=4)
+dev.off()
+```