190819_edgeR_CaT_Seq.Rmd

---
title: "190819 - edgeR Huh7 CRISPRa R analysis notebook"
output: html_notebook
---

####Install packages:
```{r}
#install.packages("calibrate")
#source("https://bioconductor.org/biocLite.R")
#biocLite("edgeR")
#biocLite("biomaRt")
```

####Load packages:
```{r}
library(edgeR)
library(biomaRt)
library(calibrate)
```

```{r}
sessionInfo()
```
R version 3.4.3 (2017-11-30)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.6

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

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

other attached packages:
[1] calibrate_1.7.2 MASS_7.3-51.4   biomaRt_2.34.2  edgeR_3.20.9    limma_3.34.9   

loaded via a namespace (and not attached):
 [1] httr_1.4.0                 Biobase_2.38.0             bit64_0.9-7                splines_3.4.3              Formula_1.2-3              assertthat_0.2.1           stats4_3.4.3              
 [8] latticeExtra_0.6-28        blob_1.1.1                 GenomeInfoDbData_1.0.0     progress_1.2.2             yaml_2.2.0                 pillar_1.4.1               RSQLite_2.1.1             
[15] backports_1.1.4            lattice_0.20-38            glue_1.3.1                 digest_0.6.19              GenomicRanges_1.30.3       RColorBrewer_1.1-2         XVector_0.18.0            
[22] checkmate_1.9.3            colorspace_1.4-1           htmltools_0.3.6            Matrix_1.2-17              DESeq2_1.18.1              XML_3.98-1.20              pkgconfig_2.0.2           
[29] genefilter_1.60.0          zlibbioc_1.24.0            purrr_0.3.2                xtable_1.8-4               scales_1.0.0               BiocParallel_1.12.0        htmlTable_1.13.1          
[36] tibble_2.1.3               annotate_1.56.2            IRanges_2.12.0             ggplot2_3.2.0              SummarizedExperiment_1.8.1 nnet_7.3-12                BiocGenerics_0.24.0       
[43] lazyeval_0.2.2             survival_2.44-1.1          magrittr_1.5               crayon_1.3.4               memoise_1.1.0              foreign_0.8-71             prettyunits_1.0.2         
[50] tools_3.4.3                data.table_1.12.0          hms_0.4.2                  matrixStats_0.54.0         stringr_1.4.0              S4Vectors_0.16.0           munsell_0.5.0             
[57] locfit_1.5-9.1             cluster_2.0.7-1            DelayedArray_0.4.1         AnnotationDbi_1.40.0       compiler_3.4.3             GenomeInfoDb_1.14.0        rlang_0.3.4               
[64] grid_3.4.3                 RCurl_1.95-4.12            rstudioapi_0.10            htmlwidgets_1.3            bitops_1.0-6               base64enc_0.1-3            gtable_0.3.0              
[71] DBI_1.0.0                  R6_2.4.0                   gridExtra_2.3              knitr_1.23                 dplyr_0.8.1                bit_1.1-14                 Hmisc_4.2-0               
[78] stringi_1.4.3              parallel_3.4.3             Rcpp_1.0.1                 geneplotter_1.56.0         rpart_4.1-15               acepack_1.4.1              tidyselect_0.2.5          
[85] xfun_0.7                  

####Load data:
```{r}
counts <- read.delim("~/Desktop/190425_EXP_18_CA6218_nsa_2novel.clean.readCount", row.names = 1)
```

####Select single sample (SS) data, i.e. all NT-Ctrl samples and all from one conditions
```{r}
cols_lncFARED1=c(7,8,13,10,18,20)
lncFARED1_SS=counts[,cols_lncFARED1]

cols_lncFARED2=c(7,9,13,11,18,21)
lncFARED2_SS=counts[,cols_lncFARED2]

cols_lncFARED3=c(3,1,18,17)
lncFARED3_SS=counts[,cols_lncFARED3]

cols_lncFARED4=c(3,4,13,15,18,16)
lncFARED4_SS=counts[,cols_lncFARED4]

cols_lncFARED5=c(3,2,7,6,13,12)
lncFARED5_SS=counts[,cols_lncFARED5]

cols_lncFARED6=c(7,5,13,14,18,19)
lncFARED6_SS=counts[,cols_lncFARED6]
```

####Selection of all Huh7 samples
```{r}
cols_Huh7=c(1:21)
Huh7_counts=counts[,cols_Huh7]
```


####Apply cutoff. Use at least 1 count in 3 or more samples:
```{r}
lncFARED1_counts_wco=lncFARED1_SS[(rowSums(lncFARED1_SS>0)>=3),]
lncFARED2_counts_wco=lncFARED2_SS[(rowSums(lncFARED2_SS>0)>=3),]
lncFARED3_counts_wco=lncFARED3_SS[(rowSums(lncFARED3_SS>0)>=2),]
lncFARED4_counts_wco=lncFARED4_SS[(rowSums(lncFARED4_SS>0)>=3),]
lncFARED5_counts_wco=lncFARED5_SS[(rowSums(lncFARED5_SS>0)>=3),]
lncFARED6_counts_wco=lncFARED6_SS[(rowSums(lncFARED6_SS>0)>=3),]
```

#edgeR

####DE analysis with the edgeR package (3 replicates)
```{r}
#First we put the data into a DGEList object:
y_lncFARED1=DGEList(lncFARED1_counts_wco)
y_lncFARED2=DGEList(lncFARED2_counts_wco)
y_lncFARED4=DGEList(lncFARED4_counts_wco)
y_lncFARED5=DGEList(lncFARED5_counts_wco)
y_lncFARED6=DGEList(lncFARED6_counts_wco)

#Then we conduct TMM normalization to account for compositional difference between the RNA-seq libraries:
y_lncFARED1=calcNormFactors(y_lncFARED1)
y_lncFARED2=calcNormFactors(y_lncFARED2)
y_lncFARED4=calcNormFactors(y_lncFARED4)
y_lncFARED5=calcNormFactors(y_lncFARED5)
y_lncFARED6=calcNormFactors(y_lncFARED6)

#The experimental design is defined:
Experiment_SS = factor(c(rep(c(1,2,3), each=2)))
treatment_SS = factor(c("NT","target","NT","target","NT","target"))
design_SS = model.matrix(~Experiment_SS+treatment_SS)
rownames(design_SS) = factor(c("NT","target","NT","target","NT","target"))
design_SS

#Common dispersions are estimated:
y_lncFARED1 = estimateGLMCommonDisp(y_lncFARED1, design=design_SS, verbose=TRUE)
y_lncFARED2 = estimateGLMCommonDisp(y_lncFARED2, design=design_SS, verbose=TRUE)
y_lncFARED4 = estimateGLMCommonDisp(y_lncFARED4, design=design_SS, verbose=TRUE)
y_lncFARED5 = estimateGLMCommonDisp(y_lncFARED5, design=design_SS, verbose=TRUE)
y_lncFARED6 = estimateGLMCommonDisp(y_lncFARED6, design=design_SS, verbose=TRUE)

#Generalized Linear Models are fitted:
fit_lncFARED1 = glmFit(y_lncFARED1, design_SS)
fit_lncFARED2 = glmFit(y_lncFARED2, design_SS)
fit_lncFARED4 = glmFit(y_lncFARED4, design_SS)
fit_lncFARED5 = glmFit(y_lncFARED5, design_SS)
fit_lncFARED6 = glmFit(y_lncFARED6, design_SS)

#Likelihood ratio tests are conducted:
lrt_lncFARED1 <- glmLRT(fit_lncFARED1, coef=4)
lrt_lncFARED2 <- glmLRT(fit_lncFARED2, coef=4)
lrt_lncFARED4 <- glmLRT(fit_lncFARED4, coef=4)
lrt_lncFARED5 <- glmLRT(fit_lncFARED5, coef=4)
lrt_lncFARED6 <- glmLRT(fit_lncFARED6, coef=4)
```

####DE analysis with the edgeR package (2 replicates)
```{r}
y_lncFARED3=DGEList(lncFARED3_counts_wco)
y_lncFARED3=calcNormFactors(y_lncFARED3)

#special experimental design because only 2 replicates:
Experiment_2rep = factor(c(rep(c(1,2), each=2)))
treatment_2rep = factor(c("NT","target","NT","target"))
design_2rep = model.matrix(~Experiment_2rep+treatment_2rep)
rownames(design_2rep) = factor(c("NT","target","NT","target"))
design_2rep

y_lncFARED3 = estimateGLMCommonDisp(y_lncFARED3, design=design_2rep, verbose=TRUE)
fit_lncFARED3 = glmFit(y_lncFARED3, design_2rep)
lrt_lncFARED3 <- glmLRT(fit_lncFARED3, coef=3)
```

####Top DE genes are:
```{r}
topTags(lrt_lncFARED1)
print(summary(decideTestsDGE(lrt_lncFARED1)))
```

```{r}
topTags(lrt_lncFARED2)
print(summary(decideTestsDGE(lrt_lncFARED2)))
```

```{r}
topTags(lrt_lncFARED3)
print(summary(decideTestsDGE(lrt_lncFARED3)))
```
```{r}
topTags(lrt_lncFARED4)
print(summary(decideTestsDGE(lrt_lncFARED4)))
```
```{r}
topTags(lrt_lncFARED5)
print(summary(decideTestsDGE(lrt_lncFARED5)))
```
```{r}
topTags(lrt_lncFARED6)
print(summary(decideTestsDGE(lrt_lncFARED6)))
```

####Add official gene symbols to lists and write the results of the DE analysis to a table:
```{r}
results_lncFARED1_SS <- topTags(lrt_lncFARED1,n = length(y_lncFARED1$counts[,1]))
results_lncFARED2_SS <- topTags(lrt_lncFARED2,n = length(y_lncFARED2$counts[,1]))
results_lncFARED3_SS <- topTags(lrt_lncFARED3,n = length(y_lncFARED3$counts[,1]))
results_lncFARED4_SS <- topTags(lrt_lncFARED4,n = length(y_lncFARED4$counts[,1]))
results_lncFARED5_SS <- topTags(lrt_lncFARED5,n = length(y_lncFARED5$counts[,1]))
results_lncFARED6_SS <- topTags(lrt_lncFARED6,n = length(y_lncFARED6$counts[,1]))

mart <- useMart("ENSEMBL_MART_ENSEMBL", dataset = "hsapiens_gene_ensembl",host="www.ensembl.org")
ensembl2name <- getBM(attributes=c("ensembl_gene_id","external_gene_name"),mart=mart)

results_lncFARED1_SS <- merge(x=results_lncFARED1_SS$table, y=ensembl2name, by.x=0, by.y=1, all.x=TRUE)
results_lncFARED2_SS <- merge(x=results_lncFARED2_SS$table, y=ensembl2name, by.x=0, by.y=1, all.x=TRUE)
results_lncFARED3_SS <- merge(x=results_lncFARED3_SS$table, y=ensembl2name, by.x=0, by.y=1, all.x=TRUE)
results_lncFARED4_SS <- merge(x=results_lncFARED4_SS$table, y=ensembl2name, by.x=0, by.y=1, all.x=TRUE)
results_lncFARED5_SS <- merge(x=results_lncFARED5_SS$table, y=ensembl2name, by.x=0, by.y=1, all.x=TRUE)
results_lncFARED6_SS <- merge(x=results_lncFARED6_SS$table, y=ensembl2name, by.x=0, by.y=1, all.x=TRUE)

write.table(as.matrix(results_lncFARED1_SS),file="~/Desktop/190426_lncFARED1_Huh7_DEG_w2Novel",sep="\t")
write.table(as.matrix(results_lncFARED2_SS),file="~/Desktop/190426_lncFARED2_Huh7_DEG_w2Novel",sep="\t")
write.table(as.matrix(results_lncFARED3_SS),file="~/Desktop/190426_lncFARED3_Huh7_DEG_w2Novel",sep="\t")
write.table(as.matrix(results_lncFARED4_SS),file="~/Desktop/190426_lncFARED4_Huh7_DEG_w2Novel",sep="\t")
write.table(as.matrix(results_lncFARED5_SS),file="~/Desktop/190426_lncFARED5_Huh7_DEG_w2Novel",sep="\t")
write.table(as.matrix(results_lncFARED6_SS),file="~/Desktop/190426_lncFARED6_Huh7_DEG_w2Novel",sep="\t")
```

####Write the fitted (TMM) values to a table:
```{r}
fit_lncFARED1_name <- merge(x=fit_lncFARED1$fitted.values, y=ensembl2name, by.x=0, by.y=1, all.x=TRUE)
write.table(as.matrix(fit_lncFARED1_name),file="~/Desktop/190426_lncFARED1_Huh7_TMM_w2Novel",sep="\t")

fit_lncFARED2_name <- merge(x=fit_lncFARED2$fitted.values, y=ensembl2name, by.x=0, by.y=1, all.x=TRUE)
write.table(as.matrix(fit_lncFARED2_name),file="~/Desktop/190426_lncFARED2_Huh7_TMM_w2Novel",sep="\t")

fit_lncFARED3_name <- merge(x=fit_lncFARED3$fitted.values, y=ensembl2name, by.x=0, by.y=1, all.x=TRUE)
write.table(as.matrix(fit_lncFARED3_name),file="~/Desktop/190426_lncFARED3_Huh7_TMM_w2Novel",sep="\t")

fit_lncFARED4_name <- merge(x=fit_lncFARED4$fitted.values, y=ensembl2name, by.x=0, by.y=1, all.x=TRUE)
write.table(as.matrix(fit_lncFARED4_name),file="~/Desktop/190426_lncFARED4_Huh7_TMM_w2Novel",sep="\t")

fit_lncFARED5_name <- merge(x=fit_lncFARED5$fitted.values, y=ensembl2name, by.x=0, by.y=1, all.x=TRUE)
write.table(as.matrix(fit_lncFARED5_name),file="~/Desktop/190426_lncFARED5_Huh7_TMM_w2Novel",sep="\t")

fit_lncFARED6_name <- merge(x=fit_lncFARED6$fitted.values, y=ensembl2name, by.x=0, by.y=1, all.x=TRUE)
write.table(as.matrix(fit_lncFARED6_name),file="~/Desktop/190426_lncFARED6_Huh7_TMM_w2Novel",sep="\t")
```

#PCA plots

```{r}
#In order to do PCA analysis on all samples together, we first normalize the samples similar to what was done for the DE analysis (TMM in edgeR) 
Huh7_counts_wco=Huh7_counts[(rowSums(counts>0)>=3),]
y=DGEList(Huh7_counts_wco)
y <- calcNormFactors(y)
Experiment <- factor(c("5","5","5","5","8","8","8","8","8","9","9","9","9","9","9","10","10","10","10","10","10"))
Target <- factor(c("lncFARED3","lncFARED5","NT","lncFARED4","lncFARED6","lncFARED5","NT","lncFARED1","lncFARED2","lncFARED1","lncFARED2","lncFARED5","NT","lncFARED6","LINC001419","lncFARED4","lncFARED3","NT","lncFARED6","lncFARED1","lncFARED2"))
design <- model.matrix(~Experiment+Target)
rownames(design) <- colnames(y)
y <- estimateGLMCommonDisp(y, design=design, verbose=TRUE)
fit <- glmFit(y, design)
```
```{r}
#Next we library normalize the TMM values:
TMM=colSums(fit$fitted.values)

TMM_LibrNorm = sweep(fit$fitted.values, 2, TMM, `/`)*10^6

colSums(TMM_LibrNorm)
```

```{r}
##PCA on normalized values
par(pty="s") 
PCA_norm<- prcomp(t(log2(TMM_LibrNorm+1))
Matrix <- summary(PCA_norm)
col.def<-c("yellow"	,"black"	,"red"	,"grey",	"green",	"black",	"red",	"purple",	"brown",	"purple",	"brown",	"black",	"red",	"green",	"grey",	"grey",	"yellow",	"red",	"green",	"purple",	"brown")
colors <- col.def

plot(PCA_norm$x[,1],PCA_norm$x[,2],col=colors,pch=20,xlab=paste0("PC1 (", round(Matrix$importance[2,1]*100), "%)"),ylab=paste0("PC2 (", round(Matrix$importance[2,2]*100),"%)"),cex=2)
text(PCA_norm$x[,1],PCA_norm$x[,2],labels=colnames(TMM_LibrNorm), pos=2, cex=0.5)

plot(PCA_norm$x[,2],PCA_norm$x[,3],col=colors,pch=20,xlab=paste0("PC2 (", round(Matrix$importance[2,2]*100), "%)"),ylab=paste0("PC3 (", round(Matrix$importance[2,3]*100),"%)"),cex=2)
text(PCA_norm$x[,2],PCA_norm$x[,3],labels=colnames(TMM_LibrNorm), pos=2, cex=0.5)

plot(PCA_norm$x[,3],PCA_norm$x[,4],col=colors,pch=20,xlab=paste0("PC3 (", (Matrix$importance[2,3]*100), "%)"),ylab=paste0("PC4 (", (Matrix$importance[2,4]*100),"%)"),cex=2)
text(PCA_norm$x[,3],PCA_norm$x[,4],labels=colnames(TMM_LibrNorm),pos=list(1), cex=0.5)

plot(PCA_norm$x[,4],PCA_norm$x[,5],col=colors,pch=20,xlab=paste0("PC4 (", (Matrix$importance[2,4]*100), "%)"),ylab=paste0("PC5 (", (Matrix$importance[2,5]*100),"%)"),cex=2)
text(PCA_norm$x[,4],PCA_norm$x[,5],labels=colnames(TMM_LibrNorm),pos=list(1), cex=0.5)

plot(PCA_norm$x[,5],PCA_norm$x[,6],col=colors,pch=20,xlab=paste0("PC5 (", (Matrix$importance[2,5]*100), "%)"),ylab=paste0("PC6 (", (Matrix$importance[2,6]*100),"%)"),cex=2)
text(PCA_norm$x[,5],PCA_norm$x[,6],labels=colnames(TMM_LibrNorm),pos=list(1), cex=0.5)
```

```{r}
#library normalize TMM for lncFARED1 dataset
colsum=colSums(fit_lncFARED1$fitted.values)
colsum=t(as.vector(colsum))

lncFARED1_LibrNormTMM_counts = sweep(fit_lncFARED1$fitted.values, 2, colsum, `/`)*10^6

colSums(lncFARED1_LibrNormTMM_counts)


##PCA on normalized TMM-values for lncFARED1
lncFARED1_LibrNormTMM_counts_priorCount=lncFARED1_LibrNormTMM_counts+1

par(pty="s")
PCA_norm<- prcomp(t(log2(lncFARED1_LibrNormTMM_counts_priorCount)))
Matrix <- summary(PCA_norm)
col.def<-c("darkblue","darkgreen", "darkblue","darkgreen", "darkblue","darkgreen")
colors <- col.def


plot(PCA_norm$x[,1],PCA_norm$x[,2],col=colors,pch=20,xlab=paste0("PC1 (", round(Matrix$importance[2,1]*100), "%)"),ylab=paste0("PC2 (", round(Matrix$importance[2,2]*100),"%)"),cex=3, ylim=c(-40,40),xlim=c(-60,60))
text(PCA_norm$x[,1],PCA_norm$x[,2],labels=colnames(fit_lncFARED1$fitted.values), pos=list(1), cex=0.6)

plot(PCA_norm$x[,2],PCA_norm$x[,3],col=colors,pch=20,xlab=paste0("PC2 (", round(Matrix$importance[2,2]*100), "%)"),ylab=paste0("PC3 (", round(Matrix$importance[2,3]*100),"%)"),cex=3)
text(PCA_norm$x[,2],PCA_norm$x[,3],labels=colnames(fit_lncFARED1$fitted.values), pos=2, cex=0.6)

plot(PCA_norm$x[,3],PCA_norm$x[,4],col=colors,pch=20,xlab=paste0("PC3 ", round(Matrix$importance[2,3]*100), "%"),ylab=paste0("PC4 ", round(Matrix$importance[2,4]*100),"%"),cex=2)
text(PCA_norm$x[,3],PCA_norm$x[,4],labels=colnames(fit_lncFARED1$fitted.values),pos=list(1), cex=0.5)

```

```{r}
#library normalize TMM for lncFARED2 dataset
colsum=colSums(fit_lncFARED2$fitted.values)
colsum=t(as.vector(colsum))

lncFARED2_LibrNormTMM_counts = sweep(fit_lncFARED2$fitted.values, 2, colsum, `/`)*10^6

colSums(lncFARED2_LibrNormTMM_counts)


##PCA on normalized TMM-values for lncFARED2
lncFARED2_LibrNormTMM_counts_priorCount=lncFARED2_LibrNormTMM_counts+1

par(pty="s") #make the plot square
PCA_norm<- prcomp(t(log2(lncFARED2_LibrNormTMM_counts_priorCount)))
Matrix <- summary(PCA_norm)
col.def<-c("darkblue","darkgreen", "darkblue","darkgreen", "darkblue","darkgreen")
colors <- col.def


plot(PCA_norm$x[,1],PCA_norm$x[,2],col=colors,pch=20,xlab=paste0("PC1 (", round(Matrix$importance[2,1]*100), "%)"),ylab=paste0("PC2 (", round(Matrix$importance[2,2]*100),"%)"),cex=3, ylim=c(-40,40),xlim=c(-60,60))
text(PCA_norm$x[,1],PCA_norm$x[,2],labels=colnames(fit_lncFARED2$fitted.values), pos=list(1), cex=0.6)

plot(PCA_norm$x[,2],PCA_norm$x[,3],col=colors,pch=20,xlab=paste0("PC2 (", round(Matrix$importance[2,2]*100), "%)"),ylab=paste0("PC3 (", round(Matrix$importance[2,3]*100),"%)"),cex=3)
text(PCA_norm$x[,2],PCA_norm$x[,3],labels=colnames(fit_lncFARED2$fitted.values), pos=2, cex=0.6)

plot(PCA_norm$x[,3],PCA_norm$x[,4],col=colors,pch=20,xlab=paste0("PC3 ", round(Matrix$importance[2,3]*100), "%"),ylab=paste0("PC4 ", round(Matrix$importance[2,4]*100),"%"),cex=2)
text(PCA_norm$x[,3],PCA_norm$x[,4],labels=colnames(fit_lncFARED2$fitted.values),pos=list(1), cex=0.5)

```

```{r}
#library normalize TMM for lncFARED3 dataset
colsum=colSums(fit_lncFARED3$fitted.values)
colsum=t(as.vector(colsum))

lncFARED3_LibrNormTMM_counts = sweep(fit_lncFARED3$fitted.values, 2, colsum, `/`)*10^6

colSums(lncFARED3_LibrNormTMM_counts)


##PCA on normalized TMM-values for lncFARED3
lncFARED3_LibrNormTMM_counts_priorCount=lncFARED3_LibrNormTMM_counts+1

par(pty="s") #make the plot square
PCA_norm<- prcomp(t(log2(lncFARED3_LibrNormTMM_counts_priorCount)))
Matrix <- summary(PCA_norm)
col.def<-c("darkblue","darkgreen", "darkblue","darkgreen", "darkblue","darkgreen")
colors <- col.def


plot(PCA_norm$x[,1],PCA_norm$x[,2],col=colors,pch=20,xlab=paste0("PC1 (", round(Matrix$importance[2,1]*100), "%)"),ylab=paste0("PC2 (", round(Matrix$importance[2,2]*100),"%)"),cex=3, ylim=c(-40,40),xlim=c(-60,60))
text(PCA_norm$x[,1],PCA_norm$x[,2],labels=colnames(fit_lncFARED3$fitted.values), pos=list(1), cex=0.6)

plot(PCA_norm$x[,2],PCA_norm$x[,3],col=colors,pch=20,xlab=paste0("PC2 (", round(Matrix$importance[2,2]*100), "%)"),ylab=paste0("PC3 (", round(Matrix$importance[2,3]*100),"%)"),cex=3)
text(PCA_norm$x[,2],PCA_norm$x[,3],labels=colnames(fit_lncFARED3$fitted.values), pos=2, cex=0.6)

plot(PCA_norm$x[,3],PCA_norm$x[,4],col=colors,pch=20,xlab=paste0("PC3 ", round(Matrix$importance[2,3]*100), "%"),ylab=paste0("PC4 ", round(Matrix$importance[2,4]*100),"%"),cex=2)
text(PCA_norm$x[,3],PCA_norm$x[,4],labels=colnames(fit_lncFARED3$fitted.values),pos=list(1), cex=0.5)

```

```{r}
#library normalize TMM for lncFARED4 dataset
colsum=colSums(fit_lncFARED4$fitted.values)
colsum=t(as.vector(colsum))

lncFARED4_LibrNormTMM_counts = sweep(fit_lncFARED4$fitted.values, 2, colsum, `/`)*10^6

colSums(lncFARED4_LibrNormTMM_counts)


##PCA on normalized TMM-values for lncFARED4
lncFARED4_LibrNormTMM_counts_priorCount=lncFARED4_LibrNormTMM_counts+1

par(pty="s") #make the plot square
PCA_norm<- prcomp(t(log2(lncFARED4_LibrNormTMM_counts_priorCount)))
Matrix <- summary(PCA_norm)
col.def<-c("darkblue","darkgreen", "darkblue","darkgreen", "darkblue","darkgreen")
colors <- col.def


plot(PCA_norm$x[,1],PCA_norm$x[,2],col=colors,pch=20,xlab=paste0("PC1 (", round(Matrix$importance[2,1]*100), "%)"),ylab=paste0("PC2 (", round(Matrix$importance[2,2]*100),"%)"),cex=3, ylim=c(-40,40),xlim=c(-60,60))
text(PCA_norm$x[,1],PCA_norm$x[,2],labels=colnames(fit_lncFARED4$fitted.values), pos=list(1), cex=0.6)

plot(PCA_norm$x[,2],PCA_norm$x[,3],col=colors,pch=20,xlab=paste0("PC2 (", round(Matrix$importance[2,2]*100), "%)"),ylab=paste0("PC3 (", round(Matrix$importance[2,3]*100),"%)"),cex=3)
text(PCA_norm$x[,2],PCA_norm$x[,3],labels=colnames(fit_lncFARED4$fitted.values), pos=2, cex=0.6)

plot(PCA_norm$x[,3],PCA_norm$x[,4],col=colors,pch=20,xlab=paste0("PC3 ", round(Matrix$importance[2,3]*100), "%"),ylab=paste0("PC4 ", round(Matrix$importance[2,4]*100),"%"),cex=2)
text(PCA_norm$x[,3],PCA_norm$x[,4],labels=colnames(fit_lncFARED4$fitted.values),pos=list(1), cex=0.5)

```

```{r}
#library normalize TMM for lncFARED5 dataset
colsum=colSums(fit_lncFARED5$fitted.values)
colsum=t(as.vector(colsum))

lncFARED5_LibrNormTMM_counts = sweep(fit_lncFARED5$fitted.values, 2, colsum, `/`)*10^6

colSums(lncFARED5_LibrNormTMM_counts)


##PCA on normalized TMM-values for lncFARED5
lncFARED5_LibrNormTMM_counts_priorCount=lncFARED5_LibrNormTMM_counts+1

par(pty="s") #make the plot square
PCA_norm<- prcomp(t(log2(lncFARED5_LibrNormTMM_counts_priorCount)))
Matrix <- summary(PCA_norm)
col.def<-c("darkblue","darkgreen", "darkblue","darkgreen", "darkblue","darkgreen")
colors <- col.def


plot(PCA_norm$x[,1],PCA_norm$x[,2],col=colors,pch=20,xlab=paste0("PC1 (", round(Matrix$importance[2,1]*100), "%)"),ylab=paste0("PC2 (", round(Matrix$importance[2,2]*100),"%)"),cex=3, ylim=c(-40,40),xlim=c(-60,60))
text(PCA_norm$x[,1],PCA_norm$x[,2],labels=colnames(fit_lncFARED5$fitted.values), pos=list(1), cex=0.6)

plot(PCA_norm$x[,2],PCA_norm$x[,3],col=colors,pch=20,xlab=paste0("PC2 (", round(Matrix$importance[2,2]*100), "%)"),ylab=paste0("PC3 (", round(Matrix$importance[2,3]*100),"%)"),cex=3)
text(PCA_norm$x[,2],PCA_norm$x[,3],labels=colnames(fit_lncFARED5$fitted.values), pos=2, cex=0.6)

plot(PCA_norm$x[,3],PCA_norm$x[,4],col=colors,pch=20,xlab=paste0("PC3 ", round(Matrix$importance[2,3]*100), "%"),ylab=paste0("PC4 ", round(Matrix$importance[2,4]*100),"%"),cex=2)
text(PCA_norm$x[,3],PCA_norm$x[,4],labels=colnames(fit_lncFARED5$fitted.values),pos=list(1), cex=0.5)

```

```{r}
#library normalize TMM for lncFARED6 dataset
colsum=colSums(fit_lncFARED6$fitted.values)
colsum=t(as.vector(colsum))

lncFARED6_LibrNormTMM_counts = sweep(fit_lncFARED6$fitted.values, 2, colsum, `/`)*10^6

colSums(lncFARED6_LibrNormTMM_counts)


##PCA on normalized TMM-values for lncFARED6
lncFARED6_LibrNormTMM_counts_priorCount=lncFARED6_LibrNormTMM_counts+1

par(pty="s") #make the plot square
PCA_norm<- prcomp(t(log2(lncFARED6_LibrNormTMM_counts_priorCount)))
Matrix <- summary(PCA_norm)
col.def<-c("darkblue","darkgreen", "darkblue","darkgreen", "darkblue","darkgreen")
colors <- col.def


plot(PCA_norm$x[,1],PCA_norm$x[,2],col=colors,pch=20,xlab=paste0("PC1 (", round(Matrix$importance[2,1]*100), "%)"),ylab=paste0("PC2 (", round(Matrix$importance[2,2]*100),"%)"),cex=3, ylim=c(-40,40),xlim=c(-60,60))
text(PCA_norm$x[,1],PCA_norm$x[,2],labels=colnames(fit_lncFARED6$fitted.values), pos=list(1), cex=0.6)

plot(PCA_norm$x[,2],PCA_norm$x[,3],col=colors,pch=20,xlab=paste0("PC2 (", round(Matrix$importance[2,2]*100), "%)"),ylab=paste0("PC3 (", round(Matrix$importance[2,3]*100),"%)"),cex=3)
text(PCA_norm$x[,2],PCA_norm$x[,3],labels=colnames(fit_lncFARED6$fitted.values), pos=2, cex=0.6)

plot(PCA_norm$x[,3],PCA_norm$x[,4],col=colors,pch=20,xlab=paste0("PC3 ", round(Matrix$importance[2,3]*100), "%"),ylab=paste0("PC4 ", round(Matrix$importance[2,4]*100),"%"),cex=2)
text(PCA_norm$x[,3],PCA_norm$x[,4],labels=colnames(fit_lncFARED6$fitted.values),pos=list(1), cex=0.5)

```

#Volcano plots

####lncFARED1 volcano:
```{r}
#import revious analysis results
results_lncFARED1 = read.delim("~/Desktop/190426_lncFARED1_edgeR_DEG.tsv", row.names = 1, header = TRUE)
```

```{r}
# set values outside axis limits to Inf
#upper_FDR_limit = 10^(-15)
upper_FDR_limit = 10
upper_logFC_limit = 6
lower_logFC_limit = -6

results_lncFARED1_limits = results_lncFARED1
results_lncFARED1_limits$logFDR = -log10(results_lncFARED1_limits$FDR)
results_lncFARED1_limits$logFDR[results_lncFARED1_limits$logFDR > upper_FDR_limit] = Inf

results_lncFARED1_limits$logFC[results_lncFARED1_limits$logFC > upper_logFC_limit] = upper_logFC_limit
results_lncFARED1_limits$logFC[results_lncFARED1_limits$logFC < lower_logFC_limit] = lower_logFC_limit
```

```{r}
lncFARED1_vp=results_lncFARED1_limits

#load list of genes to highlight
highlight = read.delim("~/Desktop/200525_gene_highlight_list.txt", row.names = 1, header=TRUE)

lncFARED1_vp_sub = subset(lncFARED1_vp, rownames(lncFARED1_vp) %in% rownames(highlight_lncFARED1))
lncFARED1_vp_sub = subset(lncFARED1_vp_sub, abs(logFC)>log2(0) & FDR<0.01)
lncFARED1_vp_sub = merge(x=lncFARED1_vp_sub, y=highlight_lncFARED1, by.x=0, by.y=0, all.x=TRUE)
lncFARED1_vp_sub = lncFARED1_vp_sub[,2:15]

lncFARED1_vp_plot = ggplot(lncFARED1_vp) +
  geom_point(
      data = lncFARED1_vp,
      aes(x = logFC, y = logFDR),
      alpha = 0.2,
      color = "grey",
      cex = 1
    ) +
  geom_point(
      data = subset(lncFARED1_vp, abs(logFC)>log2(0) & FDR<0.01),
      aes(x = logFC, y = logFDR),
      alpha = 0.5,
      color = "#20CBF8",
      cex = 1
    ) +
  geom_point(
      data = lncFARED1_vp_sub,
      aes(x = logFC, y = logFDR),
      fill = "#20CBF8",
      color = "black",
      pch = 21,
      cex = 3
    ) +
    ylim(0,10) +
    xlim(-6,6) +
    theme_bw(base_size = 14) +
    theme(axis.text = element_blank()) +
    scale_y_continuous(breaks=c(0,5,10)) +
    ylab("") +
    xlab("")

lncFARED1_vp_plot
```

####lncFARED4 volcano:
```{r}
#import revious analysis results
results_lncFARED4 = read.delim("~/Desktop/190426_lncFARED4_edgeR_DEG.tsv", row.names = 1, header = TRUE)
```

```{r}
# set values outside axis limits to Inf
upper_FDR_limit = 10
upper_logFC_limit = 6
lower_logFC_limit = -6

results_lncFARED4_limits = results_lncFARED4
results_lncFARED4_limits$logFDR = -log10(results_lncFARED4_limits$FDR)
results_lncFARED4_limits$logFDR[results_lncFARED4_limits$logFDR > upper_FDR_limit] = Inf

results_lncFARED4_limits$logFC[results_lncFARED4_limits$logFC > upper_logFC_limit] = upper_logFC_limit
results_lncFARED4_limits$logFC[results_lncFARED4_limits$logFC < lower_logFC_limit] = lower_logFC_limit
```

```{r}
lncFARED4_vp=results_lncFARED4_limits

#load list of genes to highlight
highlight = read.delim("~/Desktop/200525_gene_highlight_list.txt", row.names = 1, header=TRUE)

lncFARED4_vp_sub = subset(lncFARED4_vp, rownames(lncFARED4_vp) %in% rownames(highlight_lncFARED4))
lncFARED4_vp_sub = subset(lncFARED4_vp_sub, abs(logFC)>log2(2) & FDR<0.01)
lncFARED4_vp_sub = merge(x=lncFARED4_vp_sub, y=highlight_lncFARED4, by.x=0, by.y=0, all.x=TRUE)
lncFARED4_vp_sub = lncFARED4_vp_sub[,2:15]

lncFARED4_vp_plot = ggplot(lncFARED4_vp) +
  geom_point(
      data = lncFARED4_vp,
      aes(x = logFC, y = logFDR),
      alpha = 0.2,
      color = "grey",
      cex = 1
    ) +
  geom_point(
      data = subset(lncFARED4_vp, abs(logFC)>log2(0) & FDR<0.01),
      aes(x = logFC, y = logFDR),
      alpha = 0.5,
      color = "#FC2600",
      cex = 1
    ) +
  geom_point(
      data = lncFARED4_vp_sub,
      aes(x = logFC, y = logFDR),
      fill = "#FC2600",
      color = "black",
      pch = 21,
      cex = 3
    ) +
    ylim(0,10) +
    xlim(-6,6) +
    theme_bw(base_size = 14) +
    theme(axis.text = element_blank()) +
    scale_y_continuous(breaks=c(0,5,10)) +
    ylab("") +
    xlab("")

lncFARED4_vp_plot
```

####lncFARED3 volcano:
```{r}
#import revious analysis results
results_lncFARED3 = read.delim("~/Desktop/190426_lncFARED3_edgeR_DEG.tsv", row.names = 1, header = TRUE)
```

```{r}
# set values outside axis limits to Inf
upper_FDR_limit = 10
upper_logFC_limit = 6
lower_logFC_limit = -6

results_lncFARED3_limits = results_lncFARED3
results_lncFARED3_limits$logFDR = -log10(results_lncFARED3_limits$FDR)
results_lncFARED3_limits$logFDR[results_lncFARED3_limits$logFDR > upper_FDR_limit] = Inf

results_lncFARED3_limits$logFC[results_lncFARED3_limits$logFC > upper_logFC_limit] = upper_logFC_limit
results_lncFARED3_limits$logFC[results_lncFARED3_limits$logFC < lower_logFC_limit] = lower_logFC_limit

```

```{r}
lncFARED3_vp=results_lncFARED3_limits

#load list of genes to highlight
highlight = read.delim("~/Desktop/200525_gene_highlight_list.txt", row.names = 1, header=TRUE)

lncFARED3_vp_sub = subset(lncFARED3_vp, rownames(lncFARED3_vp) %in% rownames(highlight_lncFARED3))
lncFARED3_vp_sub = subset(lncFARED3_vp_sub, abs(logFC)>log2(2) & FDR<0.01)
lncFARED3_vp_sub = merge(x=lncFARED3_vp_sub, y=highlight_lncFARED3, by.x=0, by.y=0, all.x=TRUE)
lncFARED3_vp_sub = lncFARED3_vp_sub[,2:15]

lncFARED3_vp_plot = ggplot(lncFARED3_vp) +
  geom_point(
      data = lncFARED3_vp,
      aes(x = logFC, y = logFDR),
      alpha = 0.2,
      color = "grey",
      cex = 1
    ) +
  geom_point(
      data = subset(lncFARED3_vp, abs(logFC)>log2(0) & FDR<0.01),
      aes(x = logFC, y = logFDR),
      alpha = 0.5,
      color = "#20CBF8",
      cex = 1
    ) +
  geom_point(
      data = lncFARED3_vp_sub,
      aes(x = logFC, y = logFDR),
      fill = "#20CBF8",
      color = "black",
      pch = 21,
      cex = 3
    ) +
    ylim(0,10) +
    xlim(-6,6) +
    theme_bw(base_size = 14) +
    theme(axis.text = element_blank()) +
    scale_y_continuous(breaks=c(0,5,10)) +
    ylab("") +
    xlab("")

lncFARED3_vp_plot
```

####lncFARED5 volcano:
```{r}
#import revious analysis results

results_lncFARED5 = read.delim("~/Desktop/190426_lncFARED5_edgeR_DEG.tsv", row.names = 1, header = TRUE)
```

```{r}
# set values outside axis limits to Inf
upper_FDR_limit = 10
upper_logFC_limit = 6
lower_logFC_limit = -6

results_lncFARED5_limits = results_lncFARED5
results_lncFARED5_limits$logFDR = -log10(results_lncFARED5_limits$FDR)
results_lncFARED5_limits$logFDR[results_lncFARED5_limits$logFDR > upper_FDR_limit] = Inf

results_lncFARED5_limits$logFC[results_lncFARED5_limits$logFC > upper_logFC_limit] = upper_logFC_limit
results_lncFARED5_limits$logFC[results_lncFARED5_limits$logFC < lower_logFC_limit] = lower_logFC_limit

```

```{r}
lncFARED5_vp=results_lncFARED5_limits

#load list of genes to highlight
highlight = read.delim("~/Desktop/200525_gene_highlight_list.txt", row.names = 1, header=TRUE)

lncFARED5_vp_sub = subset(lncFARED5_vp, rownames(lncFARED5_vp) %in% rownames(highlight_lncFARED5))
lncFARED5_vp_sub = subset(lncFARED5_vp_sub, abs(logFC)>log2(2) & FDR<0.05)
lncFARED5_vp_sub = merge(x=lncFARED5_vp_sub, y=highlight_lncFARED5, by.x=0, by.y=0, all.x=TRUE)
lncFARED5_vp_sub = lncFARED5_vp_sub[,2:15]

lncFARED5_vp_plot = ggplot(lncFARED5_vp) +
  geom_point(
      data = lncFARED5_vp,
      aes(x = logFC, y = logFDR),
      color = "grey",
      alpha = 0.2,
      cex = 1
    ) +
  geom_point(
      data = subset(lncFARED5_vp, abs(logFC)>log2(0) & FDR<0.05),
      aes(x = logFC, y = logFDR),
      alpha = 0.5,
      color = "#AB46DC",
      cex = 1
    ) +
  geom_point(
      data = lncFARED5_vp_sub,
      aes(x = logFC, y = logFDR),
      fill = "#AB46DC",
      color = "black",
      pch = 21,
      cex = 3
    ) +
    ylim(0,10) +
    xlim(-6,6) +
    theme_bw(base_size = 14) +
    theme(axis.text = element_blank()) +
    scale_y_continuous(breaks=c(0,5,10)) +
    ylab("") +
    xlab("")

lncFARED5_vp_plot
```

```{r}
#print all volcanos to one figure
g = plot_grid(ncol = 2,
          lncFARED1_vp_plot,
          lncFARED4_vp_plot,
          lncFARED3_vp_plot,
          lncFARED5_vp_plot
)
g
```