Skip to content
This repository has been archived by the owner on Oct 23, 2023. It is now read-only.

Commit

Permalink
Skylers changes for centos7 / RNASeq
Browse files Browse the repository at this point in the history
  • Loading branch information
@kopardev
kopardev committed Jun 12, 2018
1 parent 3fa25a9 commit 24784be
Show file tree
Hide file tree
Showing 94 changed files with 540 additions and 382 deletions.
27 changes: 24 additions & 3 deletions Results-template/Scripts/Deseq2Report.Rmd
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@

---
title: "DESeq2 results"
author: "CCBR RNAseq pipeline"
Expand Down Expand Up @@ -143,12 +144,32 @@ Phenotype=sampleinfo$condition
cell_rep=sampleinfo$label
tedf1$group = as.factor(Phenotype)
pc1 = round(pca$sdev[1]^2/sum(pca$sdev^2)*100,2)
pc2 = round(pca$sdev[2]^2/sum(pca$sdev^2)*100,2)
pc3 = round(pca$sdev[3]^2/sum(pca$sdev^2)*100,2)
pcafactor = as.factor(sampleinfo$condition)
library(RColorBrewer)
col <- brewer.pal(nlevels(pcafactor), "Paired")
p <- plot_ly(as.data.frame(pca$x[,1:3]), x = ~PC1, y = ~PC2, z = ~PC3, color = pcafactor, colors = col, hoverinfo="text",
hovertext = ~sampleinfo$label) %>%
add_markers() %>%
layout(title = "PCA plot",
scene = list(xaxis = list(title = paste0("PC1 (",pc1,"%)")),
yaxis = list(title = paste0("PC2 (",pc2,"%)")),
zaxis = list(title = paste0("PC3 (",pc3,"%)"))))
p
# plot(pca,type="lines") #Decide how many PC's are relevant for plotting
#pca$x[,1:3] #look at first 3 PC's
plot3d(pca$x[,1:3],col = as.integer(tedf1$group),type="s",size=2)
group.v<-as.vector(cell_rep)
text3d(pca$x, pca$y, pca$z, group.v, cex=1.0, adj = 1.2)
#plot3d(pca$x[,1:3],col = as.integer(tedf1$group),type="s",size=2)
#group.v<-as.vector(cell_rep)
#text3d(pca$x, pca$y, pca$z, group.v, cex=1.0, adj = 1.2)
# rgl.postscript("pca3d_DESeq2.pdf","pdf")
```
Expand Down
27 changes: 24 additions & 3 deletions Results-template/Scripts/EdgerReport.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -181,12 +181,33 @@ Phenotype=sampleinfo$condition
cell_rep=sampleinfo$label
tedf1$group = as.factor(Phenotype)
pc1 = round(pca$sdev[1]^2/sum(pca$sdev^2)*100,2)
pc2 = round(pca$sdev[2]^2/sum(pca$sdev^2)*100,2)
pc3 = round(pca$sdev[3]^2/sum(pca$sdev^2)*100,2)
pcafactor = as.factor(sampleinfo$condition)
library(RColorBrewer)
col <- brewer.pal(nlevels(pcafactor), "Paired")
p <- plot_ly(as.data.frame(pca$x[,1:3]), x = ~PC1, y = ~PC2, z = ~PC3, color = pcafactor, colors = col, hoverinfo="text",
hovertext = ~sampleinfo$label) %>%
add_markers() %>%
layout(title = "PCA plot",
scene = list(xaxis = list(title = paste0("PC1 (",pc1,"%)")),
yaxis = list(title = paste0("PC2 (",pc2,"%)")),
zaxis = list(title = paste0("PC3 (",pc3,"%)"))))
p
# plot(pca,type="lines") #Decide how many PC's are relevant for plotting
#pca$x[,1:3] #look at first 3 PC's
plot3d(pca$x[,1:3],col = as.integer(tedf1$group),type="s",size=2)
group.v<-as.vector(cell_rep)
text3d(pca$x, pca$y, pca$z, group.v, cex=1.0, adj = 1.2)
#plot3d(pca$x[,1:3],col = as.integer(tedf1$group),type="s",size=2)
#group.v<-as.vector(cell_rep)
#text3d(pca$x, pca$y, pca$z, group.v, cex=1.0, adj = 1.2)
# rgl.postscript("pca3d_edgeR.pdf","pdf")
```
Expand Down
115 changes: 90 additions & 25 deletions Results-template/Scripts/PcaReport.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,25 +116,29 @@ ddslog2= cpm(dds.ndata,log=TRUE,normalized.lib.sizes=FALSE,prior.count=0.5)
### Before Normalization

```{r, echo=FALSE, warning=FALSE,message=FALSE}
print(ggplot(melt(as.data.frame(rawlog2))) + geom_density(aes(x = value,colour = variable)) + labs(x = NULL) + theme(legend.position='right') + scale_x_log10())
beforehist <- ggplotly(ggplot(melt(as.data.frame(rawlog2))) + geom_density(aes(x = value,colour = variable)) + labs(x = NULL) + theme(legend.position='right') + scale_x_log10())
beforehist
```

### TMM

```{r, echo=FALSE, warning=FALSE,message=FALSE}
print(ggplot(melt(as.data.frame(ylog2))) + geom_density(aes(x = value,colour = variable)) + labs(x = NULL) + theme(legend.position='right') + scale_x_log10())
tmmhist <- ggplotly(ggplot(melt(as.data.frame(ylog2))) + geom_density(aes(x = value,colour = variable)) + labs(x = NULL) + theme(legend.position='right') + scale_x_log10())
tmmhist
```

### DESeq2

```{r, echo=FALSE, warning=FALSE,message=FALSE}
print(ggplot(melt(as.data.frame(ddslog2))) + geom_density(aes(x = value,colour = variable)) + labs(x = NULL) + theme(legend.position='right') + scale_x_log10())
deshist <- ggplotly(ggplot(melt(as.data.frame(ddslog2))) + geom_density(aes(x = value,colour = variable)) + labs(x = NULL) + theme(legend.position='right') + scale_x_log10())
deshist
```

### Limma

```{r, echo=FALSE, warning=FALSE,message=FALSE}
print(ggplot(melt(as.data.frame(v1$E))) + geom_density(aes(x = value,colour = variable)) + labs(x = NULL) + theme(legend.position='right') + scale_x_log10())
limmahist <- ggplotly(ggplot(melt(as.data.frame(v1$E))) + geom_density(aes(x = value,colour = variable)) + labs(x = NULL) + theme(legend.position='right') + scale_x_log10())
limmahist
```

## **PCA Plots** {.tabset}
Expand All @@ -156,16 +160,31 @@ before.pc1 = round(before.pca$sdev[1]^2/sum(before.pca$sdev^2)*100,2)
before.pc2 = round(before.pca$sdev[2]^2/sum(before.pca$sdev^2)*100,2)
before.pc3 = round(before.pca$sdev[3]^2/sum(before.pca$sdev^2)*100,2)
pcafactor = as.factor(sampleinfo$condition)
library(RColorBrewer)
col <- brewer.pal(nlevels(pcafactor), "Paired")
p <- plot_ly(as.data.frame(before.pca$x[,1:3]), x = ~PC1, y = ~PC2, z = ~PC3, color = pcafactor, colors = col, hoverinfo="text",
hovertext = ~sampleinfo$label) %>%
add_markers() %>%
layout(title = "Before Normalization PCA plot",
scene = list(xaxis = list(title = paste0("PC1 (",before.pc1,"%)")),
yaxis = list(title = paste0("PC2 (",before.pc2,"%)")),
zaxis = list(title = paste0("PC3 (",before.pc3,"%)"))))
p
# plot(before.pca,type="lines") #Decide how many PC's are relevant for plotting
#before.pca$x[,1:3] #look at first 3 PC's
plot3d(before.pca$x[,1:3],col = as.integer(before.tedf1$group),type="s",size=2,main="PCA before normalization",xlab=paste0("PC1 (",before.pc1,"%)"),ylab=paste0("PC2 (",before.pc2,"%)"),zlab=paste0("PC3 (",before.pc3,"%)"))
group.v<-as.vector(cell_rep)
text3d(before.pca$x, before.pca$y, before.pca$z, group.v, cex=1.0, adj = 1.2)
legend3d("topright", legend = levels(sampleinfo$condition), pch = 16, col = as.numeric(as.factor(levels(sampleinfo$condition))), cex=0.5)
#plot3d(before.pca$x[,1:3],col = as.integer(before.tedf1$group),type="s",size=2,main="PCA before normalization",xlab=paste0("PC1 (",before.pc1,"%)"),ylab=paste0("PC2 (",before.pc2,"%)"),zlab=paste0("PC3 (",before.pc3,"%)"))
#group.v<-as.vector(cell_rep)
#text3d(before.pca$x, before.pca$y, before.pca$z, group.v, cex=1.0, adj = 1.2)
#legend3d("topright", legend = levels(sampleinfo$condition), pch = 16, col = as.numeric(as.factor(levels(sampleinfo$condition))), cex=0.5)
#rgl.postscript("pca3d_raw.pdf","pdf")
rgl.snapshot("pca3d_raw.png","png")
#rgl.snapshot("pca3d_raw.png","png")
```

Expand All @@ -186,16 +205,31 @@ edgeR.pc1 = round(edgeR.pca$sdev[1]^2/sum(edgeR.pca$sdev^2)*100,2)
edgeR.pc2 = round(edgeR.pca$sdev[2]^2/sum(edgeR.pca$sdev^2)*100,2)
edgeR.pc3 = round(edgeR.pca$sdev[3]^2/sum(edgeR.pca$sdev^2)*100,2)
pcafactor = as.factor(sampleinfo$condition)
library(RColorBrewer)
col <- brewer.pal(nlevels(pcafactor), "Paired")
p <- plot_ly(as.data.frame(edgeR.pca$x[,1:3]), x = ~PC1, y = ~PC2, z = ~PC3, color = pcafactor, colors = col, hoverinfo="text",
hovertext = ~sampleinfo$label) %>%
add_markers() %>%
layout(title = "edgeR PCA plot",
scene = list(xaxis = list(title = paste0("PC1 (",edgeR.pc1,"%)")),
yaxis = list(title = paste0("PC2 (",edgeR.pc2,"%)")),
zaxis = list(title = paste0("PC3 (",edgeR.pc3,"%)"))))
p
# plot(edgeR.pca,type="lines") #Decide how many PC's are relevant for plotting
#edgeR.pca$x[,1:3] #look at first 3 PC's
plot3d(edgeR.pca$x[,1:3],col = as.integer(edgeR.tedf1$group),type="s",size=2,main="PCA after TMM normalization",xlab=paste0("PC1 (",edgeR.pc1,"%)"),ylab=paste0("PC2 (",edgeR.pc2,"%)"),zlab=paste0("PC3 (",edgeR.pc3,"%)"))
group.v<-as.vector(cell_rep)
text3d(edgeR.pca$x, edgeR.pca$y, edgeR.pca$z, group.v, cex=1.0, adj = 1.2)
legend3d("topright", legend = levels(sampleinfo$condition), pch = 16, col = as.numeric(as.factor(levels(sampleinfo$condition))), cex=0.5)
#plot3d(edgeR.pca$x[,1:3],col = as.integer(edgeR.tedf1$group),type="s",size=2,main="PCA after TMM normalization",xlab=paste0("PC1 (",edgeR.pc1,"%)"),ylab=paste0("PC2 (",edgeR.pc2,"%)"),zlab=paste0("PC3 (",edgeR.pc3,"%)"))
#group.v<-as.vector(cell_rep)
#text3d(edgeR.pca$x, edgeR.pca$y, edgeR.pca$z, group.v, cex=1.0, adj = 1.2)
#legend3d("topright", legend = levels(sampleinfo$condition), pch = 16, col = as.numeric(as.factor(levels(sampleinfo$condition))), cex=0.5)
#rgl.postscript("pca3d_edgeR.pdf","pdf")
rgl.snapshot("pca3d_edgeR.png","png")
#rgl.snapshot("pca3d_edgeR.png","png")
```

Expand All @@ -221,13 +255,28 @@ deseq2.pc1 = round(deseq2.pca$sdev[1]^2/sum(deseq2.pca$sdev^2)*100,2)
deseq2.pc2 = round(deseq2.pca$sdev[2]^2/sum(deseq2.pca$sdev^2)*100,2)
deseq2.pc3 = round(deseq2.pca$sdev[3]^2/sum(deseq2.pca$sdev^2)*100,2)
pcafactor = as.factor(sampleinfo$condition)
library(RColorBrewer)
plot3d(deseq2.pca$x[,1:3],col = as.integer(deseq2.tedf1$group),type="s",size=2,main="PCA after DESeq2 normalization",xlab=paste0("PC1 (",deseq2.pc1,"%)"),ylab=paste0("PC2 (",deseq2.pc2,"%)"),zlab=paste0("PC3 (",deseq2.pc3,"%)"))
group.v<-as.vector(cell_rep)
text3d(deseq2.pca$x, deseq2.pca$y, deseq2.pca$z, group.v, cex=1.0, adj = 1.2)
legend3d("topright", legend = levels(sampleinfo$condition), pch = 16, col = as.numeric(as.factor(levels(sampleinfo$condition))), cex=0.5)
col <- brewer.pal(nlevels(pcafactor), "Paired")
p <- plot_ly(as.data.frame(deseq2.pca$x[,1:3]), x = ~PC1, y = ~PC2, z = ~PC3, color = pcafactor, colors = col, hoverinfo="text",
hovertext = ~sampleinfo$label) %>%
add_markers() %>%
layout(title = "DESeq2 PCA plot",
scene = list(xaxis = list(title = paste0("PC1 (",deseq2.pc1,"%)")),
yaxis = list(title = paste0("PC2 (",deseq2.pc2,"%)")),
zaxis = list(title = paste0("PC3 (",deseq2.pc3,"%)"))))
p
#plot3d(deseq2.pca$x[,1:3],col = as.integer(deseq2.tedf1$group),type="s",size=2,main="PCA after DESeq2 normalization",xlab=paste0("PC1 (",deseq2.pc1,"%)"),ylab=paste0("PC2 (",deseq2.pc2,"%)"),zlab=paste0("PC3 (",deseq2.pc3,"%)"))
#group.v<-as.vector(cell_rep)
#text3d(deseq2.pca$x, deseq2.pca$y, deseq2.pca$z, group.v, cex=1.0, adj = 1.2)
#legend3d("topright", legend = levels(sampleinfo$condition), pch = 16, col = as.numeric(as.factor(levels(sampleinfo$condition))), cex=0.5)
#rgl.postscript("pca3d_deseq2.pdf","pdf")
rgl.snapshot("pca3d_deseq2.png","png")
#rgl.snapshot("pca3d_deseq2.png","png")
```

Expand All @@ -249,13 +298,29 @@ limma.pc1 = round(limma.pca$sdev[1]^2/sum(limma.pca$sdev^2)*100,2)
limma.pc2 = round(limma.pca$sdev[2]^2/sum(limma.pca$sdev^2)*100,2)
limma.pc3 = round(limma.pca$sdev[3]^2/sum(limma.pca$sdev^2)*100,2)
pcafactor = as.factor(sampleinfo$condition)
library(RColorBrewer)
col <- brewer.pal(nlevels(pcafactor), "Paired")
plot3d(limma.pca$x[,1:3],col = as.integer(limma.tedf1$group),type="s",size=2,main="PCA after Limma normalization",xlab=paste0("PC1 (",limma.pc1,"%)"),ylab=paste0("PC2 (",limma.pc2,"%)"),zlab=paste0("PC3 (",limma.pc3,"%)"))
group.v<-as.vector(cell_rep)
text3d(limma.pca$x, limma.pca$y, limma.pca$z, group.v, cex=1.0, adj = 1.2)
legend3d("topright", legend = levels(sampleinfo$condition), pch = 16, col = as.numeric(as.factor(levels(sampleinfo$condition))), cex=0.5)
p <- plot_ly(as.data.frame(limma.pca$x[,1:3]), x = ~PC1, y = ~PC2, z = ~PC3, color = pcafactor, colors = col, hoverinfo="text",
hovertext = ~sampleinfo$label) %>%
add_markers() %>%
layout(title = "Limma PCA plot",
scene = list(xaxis = list(title = paste0("PC1 (",limma.pc1,"%)")),
yaxis = list(title = paste0("PC2 (",limma.pc2,"%)")),
zaxis = list(title = paste0("PC3 (",limma.pc3,"%)"))))
p
#plot3d(limma.pca$x[,1:3],col = as.integer(limma.tedf1$group),type="s",size=2,main="PCA after Limma normalization",xlab=paste0("PC1 (",limma.pc1,"%)"),ylab=paste0("PC2 (",limma.pc2,"%)"),zlab=paste0("PC3 (",limma.pc3,"%)"))
#group.v<-as.vector(cell_rep)
#text3d(limma.pca$x, limma.pca$y, limma.pca$z, group.v, cex=1.0, adj = 1.2)
#legend3d("topright", legend = levels(sampleinfo$condition), pch = 16, col = as.numeric(as.factor(levels(sampleinfo$condition))), cex=0.5)
#rgl.postscript("pca3d_limma.pdf","pdf")
rgl.snapshot("pca3d_limma.png","png")
#rgl.snapshot("pca3d_limma.png","png")
```

Expand Down Expand Up @@ -363,4 +428,4 @@ for(i in 1:length(sampleinfo$label)){
plotMD(v1$E,column=i, main=paste0("Limma ",sampleinfo$label[i]), xlim=c(-5,15), ylim=c(-15,15))
abline(h=0, col="red", lty=2, lwd=2)
}
```
```
6 changes: 3 additions & 3 deletions Results-template/Scripts/avia.pl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@
$cmd = '';
$time = 300; #amount of time (seconds) to wait between pings to AVIA

$cmd="/usr/local/bin/curl https://avia-abcc.ncifcrf.gov/apps/site/upload_viz -X POST -F user_file=\@$vcf -F user.ver=$species -F user_inputformat=bed -F user_api=cFMdtEdwm34iVzXOZ6 -F 'user_email=" . $email . "|nih.gov' --insecure -F user_id=$id";
$cmd="curl https://avia-abcc.ncifcrf.gov/apps/site/upload_viz -X POST -F user_file=\@$vcf -F user.ver=$species -F user_inputformat=bed -F user_api=cFMdtEdwm34iVzXOZ6 -F 'user_email=" . $email . "|nih.gov' --insecure -F user_id=$id";

print STDERR "Executing command: $cmd\n";

Expand All @@ -39,8 +39,8 @@
while (<G>){
chomp;
last if m/^$/;
if (($_ =~ m/INFO/) || ($_ =~ m/ERROR/)) {
if ($_ =~ m/is ready for download/) {
if (($_ =~ m/INFO/) || ($_ =~ m/completed/)) {
if ($_ =~ m/completed/) {
$status++;
}
elsif ($_ =~ m/is not accessible/) {
Expand Down
7 changes: 5 additions & 2 deletions Results-template/Scripts/genejunctioncounts.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,8 @@ myfiles=as.character(unlist(strsplit(FILES, split=" ")))

res=read.delim(myfiles[1],header=F,stringsAsFactors=F)
stranded=res[,4]

stranded=gsub("0",".",stranded)
stranded=gsub("1","+",stranded)
stranded=gsub("2","-",stranded)
filler1 = rep(".",length(res[,1]))
Expand All @@ -29,7 +31,7 @@ result=cbind(res[,1],res[,2],res[,3],filler1,filler2,stranded,res[,7]);
tmpoutfilename1=paste(substr(myfiles[1], 0, nchar(myfiles[1])-3),"bed.recoded.tab",sep="")
write.table(result,file=tmpoutfilename1,sep="\t",row.names=F,col.names=F,quote=F)
tmpoutfilename2=paste(substr(myfiles[1], 0, nchar(myfiles[1])-14),"_gencodeintersects.txt",sep="")
system(paste('module load bedtools/2.19.1; bedtools intersect -a ',"gencode.gtf.filler.bed",' -b ',tmpoutfilename1,' -wao -loj -s > ',tmpoutfilename2,sep=""))
system(paste('bedtools intersect -a ',"gencode.gtf.filler.bed",' -b ',tmpoutfilename1,' -wao -loj -s > ',tmpoutfilename2,sep=""))
tmpoutfilename3=paste(substr(myfiles[1], 0, nchar(myfiles[1])-14),"_genecounts.txt",sep="")
system(paste('awk \'{a[$7]+=$15}END{for(i in a){print i, a[i], OFS=\"\\t\"}}\' ',tmpoutfilename2," > ",tmpoutfilename3,sep=""))
tmp=read.delim(tmpoutfilename3,header=F,stringsAsFactors=F)
Expand All @@ -41,6 +43,7 @@ for(i in seq(2, length(myfiles), by = 1))
{
res=read.delim(myfiles[i],header=F,stringsAsFactors=F)
stranded=res[,4]
stranded=gsub("0",".",stranded)
stranded=gsub("1","+",stranded)
stranded=gsub("2","-",stranded)
filler1 = rep(".",length(res[,1]))
Expand All @@ -49,7 +52,7 @@ for(i in seq(2, length(myfiles), by = 1))
tmpoutfilename1=paste(substr(myfiles[i], 0, nchar(myfiles[i])-3),"bed.recoded.tab",sep="")
write.table(result,file=tmpoutfilename1,sep="\t",row.names=F,col.names=F,quote=F)
tmpoutfilename2=paste(substr(myfiles[i], 0, nchar(myfiles[i])-14),"_gencodeintersects.txt",sep="")
system(paste('module load bedtools/2.19.1; bedtools intersect -a ',"gencode.gtf.filler.bed",' -b ',tmpoutfilename1,' -wao -loj -s > ',tmpoutfilename2,sep=""))
system(paste('bedtools intersect -a ',"gencode.gtf.filler.bed",' -b ',tmpoutfilename1,' -wao -loj -s > ',tmpoutfilename2,sep=""))
tmpoutfilename3=paste(substr(myfiles[i], 0, nchar(myfiles[i])-14),"_genecounts.txt",sep="")
system(paste('awk \'{a[$7]+=$15}END{for(i in a){print i, a[i], OFS=\"\\t\"}}\' ',tmpoutfilename2," > ",tmpoutfilename3,sep=""))
tmp=read.delim(tmpoutfilename3,header=F,stringsAsFactors=F)
Expand Down
Loading

0 comments on commit 24784be

Please sign in to comment.