Dual RNA-seq of Lupin infected with P. Cinnamomi
Mark Ziemann
2021-01-22
Source: https://github.com/markziemann/phytophthora-lupin
Vew the reports: http://118.138.234.73/public/barry/
Intro
Phytophthora infected Lupinus angustifolius root.
Reference genome info
Here is the reference informtion
Lupinus angustifolius from Ensembl
Lupinus_angustifolius.LupAngTanjil_v1.0.49.gtf.gz
Lupinus_angustifolius.LupAngTanjil_v1.0.dna_sm.toplevel.fa.gz
Phytophthora cinnamomi from FungiDB 21st Jan 2021
Pc gff file was converted to gtf with gffread. fasta files were combined before indexing gtf files were combined before indexing
Data processing
Read trimming with Skewer v0.2.2 followed by STAR v2.7.7a. STAR gene level counts (reverse strand) were imported into R. Lupin and phytophthora counts were separated. DESeq2 was run to quantify differences between treatment groups.
Read in data
#conda activate R40
library("DESeq2")
library("mitch")
library("reshape2")
library("Biostrings")
library("gplots")Load counts into R
tmp <- read.table("3col.tsv.gz")
y<-as.matrix(acast(tmp, V2~V1, value.var="V3"))
colnames(y) <- sapply(strsplit(colnames(y),"_"),"[[",1)
head(y)## BS01 BS02 BS03 BS04 BS05 BS06 BS07 BS08 BS09 BS10 BS11 BS12
## ENSRNA049758817 0 0 0 0 0 0 0 0 0 0 0 0
## ENSRNA049758844 0 0 0 0 0 0 0 0 0 0 0 0
## ENSRNA049758846 0 0 0 0 0 0 0 0 0 0 0 1
## ENSRNA049758848 0 0 0 1 0 0 0 0 0 0 0 0
## ENSRNA049758851 3 1 1 1 2 1 0 0 2 2 2 3
## ENSRNA049758853 0 0 0 0 0 0 0 0 0 0 0 0
## BS13 BS14 BS15 BS16 BS17 BS18 BS19 BS20 BS21 BS22 BS23 BS24
## ENSRNA049758817 0 0 0 0 0 0 0 0 0 0 0 0
## ENSRNA049758844 0 0 0 0 0 0 0 0 0 0 0 0
## ENSRNA049758846 0 0 0 0 1 0 0 0 0 0 0 0
## ENSRNA049758848 0 0 0 0 0 0 0 0 0 0 0 0
## ENSRNA049758851 0 0 0 0 0 2 0 0 1 0 1 0
## ENSRNA049758853 0 0 0 0 0 0 0 0 0 0 0 0Load samplesheet and change the colnames.
ss <- read.table("samplesheet.tsv")
ss## V1 V2
## 1 BS01 LaH-0_1
## 2 BS02 LaH-0_2
## 3 BS03 LaH-0_3
## 4 BS04 LaH-18_1
## 5 BS05 LaH-18_2
## 6 BS06 LaH-18_3
## 7 BS07 LaH-30_1
## 8 BS08 LaH-30_2
## 9 BS09 LaH-30_3
## 10 BS10 LaH-48_1
## 11 BS11 LaH-48_2
## 12 BS12 LaH-48_3
## 13 BS13 PcHY-0_1
## 14 BS14 PcHY-0_2
## 15 BS15 PcHY-0_3
## 16 BS16 LaP-18_1
## 17 BS17 LaP-18_2
## 18 BS18 LaP-18_3
## 19 BS19 LaP-30_1
## 20 BS20 LaP-30_2
## 21 BS21 LaP-30_3
## 22 BS22 LaP-48_1
## 23 BS23 LaP-48_2
## 24 BS24 LaP-48_3ss$V1 == colnames(y)## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUEcolnames(y) <- ss$V2QC analysis of counts and split by species.
par(mar=c(5,10,4,2))
dim(y)## [1] 61301 24pc <- y[grep("PHYCI",rownames(y)),]
dim(pc)## [1] 26131 24la <- y[grep("PHYCI",rownames(y),invert=TRUE),]
dim(la)## [1] 35170 24barplot(colSums(y),horiz=TRUE,las=1,xlab="total number of assigned reads")
barplot(colSums(pc),horiz=TRUE,las=1,xlab="number of PC assigned reads")
barplot(colSums(la),horiz=TRUE,las=1,xlab="number of LA assigned reads")
barplot(log10(colSums(pc)),horiz=TRUE,las=1,xlab="log10 number of PC assigned reads")
barplot(log10(colSums(la)),horiz=TRUE,las=1,xlab="log10 number of LA assigned reads")
barplot( colSums(pc)/colSums(y), horiz=TRUE,las=1,xlab="proportion of PC/total reads" )
# revert to defaults
par(mar=c(5.1,4.1,4.1,2.1))Multidimensional scaling analysis
mds <- cmdscale(dist(t(y)))
plot(mds*1.05, xlab="Coordinate 1", ylab="Coordinate 2", type = "p",pch=19,col="lightblue",bty="n")
text(mds, labels=colnames(y))
mtext("overall MDS plot")
mds <- cmdscale(dist(t(pc)))
plot(mds*1.05, xlab="Coordinate 1", ylab="Coordinate 2", type = "p",pch=19,col="lightblue",bty="n")
text(mds, labels=colnames(pc))
mtext("Pc MDS plot")
mds <- cmdscale(dist(t(la)))
plot(mds*1.05, xlab="Coordinate 1", ylab="Coordinate 2", type = "p",pch=19,col="lightblue",bty="n")
text(mds, labels=colnames(la))
mtext("La MDS plot")
Contrasts
P cinnamomi. These will be called Pc1 to Pc6
PcHY vs LaP18. “DE early”
PcHY vs LaP30. “DE mid”
PcHY vs LaP48. “DE late”
LaPc18 vs LaPc30 “Change from early to mid”
LaPc30 vs LaPc48 “Change from mid to late”
LaPc18 vs LaPc48 “Difference between early and late”
ss$grp <- sapply(strsplit(ss$V2,"_"),"[[",1)
head(pc)## LaH-0_1 LaH-0_2 LaH-0_3 LaH-18_1 LaH-18_2 LaH-18_3 LaH-30_1
## PHYCI_100000 0 0 0 0 0 0 0
## PHYCI_100005 0 0 0 0 0 0 0
## PHYCI_100006 0 0 0 0 0 0 0
## PHYCI_100014 0 0 0 0 0 0 0
## PHYCI_100017 0 0 0 0 0 0 0
## PHYCI_100037 0 0 0 0 0 0 0
## LaH-30_2 LaH-30_3 LaH-48_1 LaH-48_2 LaH-48_3 PcHY-0_1 PcHY-0_2
## PHYCI_100000 0 0 0 0 0 714 637
## PHYCI_100005 0 0 0 0 0 81 56
## PHYCI_100006 0 0 0 0 0 54 58
## PHYCI_100014 0 0 0 0 0 55779 55655
## PHYCI_100017 0 0 0 0 0 227 267
## PHYCI_100037 0 0 0 0 0 2 0
## PcHY-0_3 LaP-18_1 LaP-18_2 LaP-18_3 LaP-30_1 LaP-30_2 LaP-30_3
## PHYCI_100000 741 28 39 12 118 139 54
## PHYCI_100005 73 1520 2497 1110 2774 1932 2063
## PHYCI_100006 66 0 0 0 0 2 1
## PHYCI_100014 52086 2102 3226 1099 15044 13468 6225
## PHYCI_100017 189 1 1 1 21 14 7
## PHYCI_100037 3 0 0 0 0 1 0
## LaP-48_1 LaP-48_2 LaP-48_3
## PHYCI_100000 158 277 123
## PHYCI_100005 1526 1219 1558
## PHYCI_100006 4 14 6
## PHYCI_100014 21762 31309 21788
## PHYCI_100017 27 52 37
## PHYCI_100037 0 0 1sspc1 <- subset(ss,grp=="PcHY-0" | grp=="LaP-18")
sspc1$trt <- as.numeric(sspc1$grp=="LaP-18")
pc1 <- pc[,which(colnames(pc) %in% sspc1$V2)]
head(pc1)## PcHY-0_1 PcHY-0_2 PcHY-0_3 LaP-18_1 LaP-18_2 LaP-18_3
## PHYCI_100000 714 637 741 28 39 12
## PHYCI_100005 81 56 73 1520 2497 1110
## PHYCI_100006 54 58 66 0 0 0
## PHYCI_100014 55779 55655 52086 2102 3226 1099
## PHYCI_100017 227 267 189 1 1 1
## PHYCI_100037 2 0 3 0 0 0head(sspc1)## V1 V2 grp trt
## 13 BS13 PcHY-0_1 PcHY-0 0
## 14 BS14 PcHY-0_2 PcHY-0 0
## 15 BS15 PcHY-0_3 PcHY-0 0
## 16 BS16 LaP-18_1 LaP-18 1
## 17 BS17 LaP-18_2 LaP-18 1
## 18 BS18 LaP-18_3 LaP-18 1sspc2 <- subset(ss,grp=="PcHY-0" | grp=="LaP-30")
sspc2$trt <- as.numeric(sspc2$grp=="LaP-30")
pc2 <- pc[,which(colnames(pc) %in% sspc2$V2)]
head(pc2)## PcHY-0_1 PcHY-0_2 PcHY-0_3 LaP-30_1 LaP-30_2 LaP-30_3
## PHYCI_100000 714 637 741 118 139 54
## PHYCI_100005 81 56 73 2774 1932 2063
## PHYCI_100006 54 58 66 0 2 1
## PHYCI_100014 55779 55655 52086 15044 13468 6225
## PHYCI_100017 227 267 189 21 14 7
## PHYCI_100037 2 0 3 0 1 0head(sspc2)## V1 V2 grp trt
## 13 BS13 PcHY-0_1 PcHY-0 0
## 14 BS14 PcHY-0_2 PcHY-0 0
## 15 BS15 PcHY-0_3 PcHY-0 0
## 19 BS19 LaP-30_1 LaP-30 1
## 20 BS20 LaP-30_2 LaP-30 1
## 21 BS21 LaP-30_3 LaP-30 1sspc3 <- subset(ss,grp=="PcHY-0" | grp=="LaP-48")
sspc3$trt <- as.numeric(sspc3$grp=="LaP-48")
pc3 <- pc[,which(colnames(pc) %in% sspc3$V2)]
head(pc3)## PcHY-0_1 PcHY-0_2 PcHY-0_3 LaP-48_1 LaP-48_2 LaP-48_3
## PHYCI_100000 714 637 741 158 277 123
## PHYCI_100005 81 56 73 1526 1219 1558
## PHYCI_100006 54 58 66 4 14 6
## PHYCI_100014 55779 55655 52086 21762 31309 21788
## PHYCI_100017 227 267 189 27 52 37
## PHYCI_100037 2 0 3 0 0 1head(sspc3)## V1 V2 grp trt
## 13 BS13 PcHY-0_1 PcHY-0 0
## 14 BS14 PcHY-0_2 PcHY-0 0
## 15 BS15 PcHY-0_3 PcHY-0 0
## 22 BS22 LaP-48_1 LaP-48 1
## 23 BS23 LaP-48_2 LaP-48 1
## 24 BS24 LaP-48_3 LaP-48 1sspc4 <- subset(ss,grp=="LaP-18" | grp=="LaP-30")
sspc4$trt <- as.numeric(sspc4$grp=="LaP-30")
pc4 <- pc[,which(colnames(pc) %in% sspc4$V2)]
head(pc4)## LaP-18_1 LaP-18_2 LaP-18_3 LaP-30_1 LaP-30_2 LaP-30_3
## PHYCI_100000 28 39 12 118 139 54
## PHYCI_100005 1520 2497 1110 2774 1932 2063
## PHYCI_100006 0 0 0 0 2 1
## PHYCI_100014 2102 3226 1099 15044 13468 6225
## PHYCI_100017 1 1 1 21 14 7
## PHYCI_100037 0 0 0 0 1 0head(sspc4)## V1 V2 grp trt
## 16 BS16 LaP-18_1 LaP-18 0
## 17 BS17 LaP-18_2 LaP-18 0
## 18 BS18 LaP-18_3 LaP-18 0
## 19 BS19 LaP-30_1 LaP-30 1
## 20 BS20 LaP-30_2 LaP-30 1
## 21 BS21 LaP-30_3 LaP-30 1sspc5 <- subset(ss,grp=="LaP-30" | grp=="LaP-48")
sspc5$trt <- as.numeric(sspc5$grp=="LaP-48")
pc5 <- pc[,which(colnames(pc) %in% sspc5$V2)]
head(pc5)## LaP-30_1 LaP-30_2 LaP-30_3 LaP-48_1 LaP-48_2 LaP-48_3
## PHYCI_100000 118 139 54 158 277 123
## PHYCI_100005 2774 1932 2063 1526 1219 1558
## PHYCI_100006 0 2 1 4 14 6
## PHYCI_100014 15044 13468 6225 21762 31309 21788
## PHYCI_100017 21 14 7 27 52 37
## PHYCI_100037 0 1 0 0 0 1head(sspc5)## V1 V2 grp trt
## 19 BS19 LaP-30_1 LaP-30 0
## 20 BS20 LaP-30_2 LaP-30 0
## 21 BS21 LaP-30_3 LaP-30 0
## 22 BS22 LaP-48_1 LaP-48 1
## 23 BS23 LaP-48_2 LaP-48 1
## 24 BS24 LaP-48_3 LaP-48 1sspc6 <- subset(ss,grp=="LaP-18" | grp=="LaP-48")
sspc6$trt <- as.numeric(sspc6$grp=="LaP-48")
pc6 <- pc[,which(colnames(pc) %in% sspc6$V2)]
head(pc6)## LaP-18_1 LaP-18_2 LaP-18_3 LaP-48_1 LaP-48_2 LaP-48_3
## PHYCI_100000 28 39 12 158 277 123
## PHYCI_100005 1520 2497 1110 1526 1219 1558
## PHYCI_100006 0 0 0 4 14 6
## PHYCI_100014 2102 3226 1099 21762 31309 21788
## PHYCI_100017 1 1 1 27 52 37
## PHYCI_100037 0 0 0 0 0 1head(sspc6)## V1 V2 grp trt
## 16 BS16 LaP-18_1 LaP-18 0
## 17 BS17 LaP-18_2 LaP-18 0
## 18 BS18 LaP-18_3 LaP-18 0
## 22 BS22 LaP-48_1 LaP-48 1
## 23 BS23 LaP-48_2 LaP-48 1
## 24 BS24 LaP-48_3 LaP-48 1Lupin. These will be called La1 to La5
LaH-0 vs LaP18 “DE early”
LaH-0 vs LaP30 “DE mid”
LaH-0 vs LaP48 “DE late”
LaPc18 vs LaPc30 “Change from early to mid”
LaPc30 vs LaPc48 “Change from mid to late”
LaPc18 vs LaPc48 “Difference between early and late”
head(la)## LaH-0_1 LaH-0_2 LaH-0_3 LaH-18_1 LaH-18_2 LaH-18_3 LaH-30_1
## ENSRNA049758817 0 0 0 0 0 0 0
## ENSRNA049758844 0 0 0 0 0 0 0
## ENSRNA049758846 0 0 0 0 0 0 0
## ENSRNA049758848 0 0 0 1 0 0 0
## ENSRNA049758851 3 1 1 1 2 1 0
## ENSRNA049758853 0 0 0 0 0 0 0
## LaH-30_2 LaH-30_3 LaH-48_1 LaH-48_2 LaH-48_3 PcHY-0_1 PcHY-0_2
## ENSRNA049758817 0 0 0 0 0 0 0
## ENSRNA049758844 0 0 0 0 0 0 0
## ENSRNA049758846 0 0 0 0 1 0 0
## ENSRNA049758848 0 0 0 0 0 0 0
## ENSRNA049758851 0 2 2 2 3 0 0
## ENSRNA049758853 0 0 0 0 0 0 0
## PcHY-0_3 LaP-18_1 LaP-18_2 LaP-18_3 LaP-30_1 LaP-30_2 LaP-30_3
## ENSRNA049758817 0 0 0 0 0 0 0
## ENSRNA049758844 0 0 0 0 0 0 0
## ENSRNA049758846 0 0 1 0 0 0 0
## ENSRNA049758848 0 0 0 0 0 0 0
## ENSRNA049758851 0 0 0 2 0 0 1
## ENSRNA049758853 0 0 0 0 0 0 0
## LaP-48_1 LaP-48_2 LaP-48_3
## ENSRNA049758817 0 0 0
## ENSRNA049758844 0 0 0
## ENSRNA049758846 0 0 0
## ENSRNA049758848 0 0 0
## ENSRNA049758851 0 1 0
## ENSRNA049758853 0 0 0ssla1 <- subset(ss,grp=="LaH-0" | grp=="LaP-18")
ssla1$trt <- as.numeric(ssla1$grp=="LaP-18")
la1 <- la[,which(colnames(pc) %in% ssla1$V2)]
head(la1)## LaH-0_1 LaH-0_2 LaH-0_3 LaP-18_1 LaP-18_2 LaP-18_3
## ENSRNA049758817 0 0 0 0 0 0
## ENSRNA049758844 0 0 0 0 0 0
## ENSRNA049758846 0 0 0 0 1 0
## ENSRNA049758848 0 0 0 0 0 0
## ENSRNA049758851 3 1 1 0 0 2
## ENSRNA049758853 0 0 0 0 0 0head(ssla1)## V1 V2 grp trt
## 1 BS01 LaH-0_1 LaH-0 0
## 2 BS02 LaH-0_2 LaH-0 0
## 3 BS03 LaH-0_3 LaH-0 0
## 16 BS16 LaP-18_1 LaP-18 1
## 17 BS17 LaP-18_2 LaP-18 1
## 18 BS18 LaP-18_3 LaP-18 1ssla2 <- subset(ss,grp=="LaH-0" | grp=="LaP-30")
ssla2$trt <- as.numeric(ssla2$grp=="LaP-30")
la2 <- la[,which(colnames(pc) %in% ssla2$V2)]
head(la2)## LaH-0_1 LaH-0_2 LaH-0_3 LaP-30_1 LaP-30_2 LaP-30_3
## ENSRNA049758817 0 0 0 0 0 0
## ENSRNA049758844 0 0 0 0 0 0
## ENSRNA049758846 0 0 0 0 0 0
## ENSRNA049758848 0 0 0 0 0 0
## ENSRNA049758851 3 1 1 0 0 1
## ENSRNA049758853 0 0 0 0 0 0head(ssla2)## V1 V2 grp trt
## 1 BS01 LaH-0_1 LaH-0 0
## 2 BS02 LaH-0_2 LaH-0 0
## 3 BS03 LaH-0_3 LaH-0 0
## 19 BS19 LaP-30_1 LaP-30 1
## 20 BS20 LaP-30_2 LaP-30 1
## 21 BS21 LaP-30_3 LaP-30 1ssla3 <- subset(ss,grp=="LaH-0" | grp=="LaP-48")
ssla3$trt <- as.numeric(ssla3$grp=="LaP-48")
la3 <- la[,which(colnames(pc) %in% ssla3$V2)]
head(la3)## LaH-0_1 LaH-0_2 LaH-0_3 LaP-48_1 LaP-48_2 LaP-48_3
## ENSRNA049758817 0 0 0 0 0 0
## ENSRNA049758844 0 0 0 0 0 0
## ENSRNA049758846 0 0 0 0 0 0
## ENSRNA049758848 0 0 0 0 0 0
## ENSRNA049758851 3 1 1 0 1 0
## ENSRNA049758853 0 0 0 0 0 0head(ssla3)## V1 V2 grp trt
## 1 BS01 LaH-0_1 LaH-0 0
## 2 BS02 LaH-0_2 LaH-0 0
## 3 BS03 LaH-0_3 LaH-0 0
## 22 BS22 LaP-48_1 LaP-48 1
## 23 BS23 LaP-48_2 LaP-48 1
## 24 BS24 LaP-48_3 LaP-48 1ssla4 <- subset(ss,grp=="LaP-18" | grp=="LaP-30")
ssla4$trt <- as.numeric(ssla4$grp=="LaP-30")
la4 <- la[,which(colnames(pc) %in% ssla4$V2)]
head(la4)## LaP-18_1 LaP-18_2 LaP-18_3 LaP-30_1 LaP-30_2 LaP-30_3
## ENSRNA049758817 0 0 0 0 0 0
## ENSRNA049758844 0 0 0 0 0 0
## ENSRNA049758846 0 1 0 0 0 0
## ENSRNA049758848 0 0 0 0 0 0
## ENSRNA049758851 0 0 2 0 0 1
## ENSRNA049758853 0 0 0 0 0 0head(ssla4)## V1 V2 grp trt
## 16 BS16 LaP-18_1 LaP-18 0
## 17 BS17 LaP-18_2 LaP-18 0
## 18 BS18 LaP-18_3 LaP-18 0
## 19 BS19 LaP-30_1 LaP-30 1
## 20 BS20 LaP-30_2 LaP-30 1
## 21 BS21 LaP-30_3 LaP-30 1ssla5 <- subset(ss,grp=="LaH-30" | grp=="LaP-48")
ssla5$trt <- as.numeric(ssla5$grp=="LaP-48")
la5 <- la[,which(colnames(pc) %in% ssla4$V2)]
head(la5)## LaP-18_1 LaP-18_2 LaP-18_3 LaP-30_1 LaP-30_2 LaP-30_3
## ENSRNA049758817 0 0 0 0 0 0
## ENSRNA049758844 0 0 0 0 0 0
## ENSRNA049758846 0 1 0 0 0 0
## ENSRNA049758848 0 0 0 0 0 0
## ENSRNA049758851 0 0 2 0 0 1
## ENSRNA049758853 0 0 0 0 0 0head(ssla5)## V1 V2 grp trt
## 7 BS07 LaH-30_1 LaH-30 0
## 8 BS08 LaH-30_2 LaH-30 0
## 9 BS09 LaH-30_3 LaH-30 0
## 22 BS22 LaP-48_1 LaP-48 1
## 23 BS23 LaP-48_2 LaP-48 1
## 24 BS24 LaP-48_3 LaP-48 1ssla6 <- subset(ss,grp=="LaH-18" | grp=="LaP-48")
ssla6$trt <- as.numeric(ssla6$grp=="LaP-48")
la6 <- la[,which(colnames(pc) %in% ssla6$V2)]
head(la6)## LaH-18_1 LaH-18_2 LaH-18_3 LaP-48_1 LaP-48_2 LaP-48_3
## ENSRNA049758817 0 0 0 0 0 0
## ENSRNA049758844 0 0 0 0 0 0
## ENSRNA049758846 0 0 0 0 0 0
## ENSRNA049758848 1 0 0 0 0 0
## ENSRNA049758851 1 2 1 0 1 0
## ENSRNA049758853 0 0 0 0 0 0head(ssla6)## V1 V2 grp trt
## 4 BS04 LaH-18_1 LaH-18 0
## 5 BS05 LaH-18_2 LaH-18 0
## 6 BS06 LaH-18_3 LaH-18 0
## 22 BS22 LaP-48_1 LaP-48 1
## 23 BS23 LaP-48_2 LaP-48 1
## 24 BS24 LaP-48_3 LaP-48 1DE function
Here is a function for differential expression.
run_de <- function(ss,xx){
y <- round(xx)
y <- y[rowMeans(y)>10,]
# MDS
mds <- cmdscale(dist(t(y)))
XMAX=max(mds[,1])*1.1
XMIN=min(mds[,1])*1.1
plot( mds , xlab="Coordinate 1", ylab="Coordinate 2",
type = "n" , xlim=c(XMIN,XMAX),main="MDS plot",bty="n")
text(mds, labels=colnames(y) )
# DE
dds <- DESeqDataSetFromMatrix(countData=y, colData = ss, design = ~ trt)
dds <- DESeq(dds)
de <- DESeq2::results(dds)
de <- de[order(de$pvalue),]
up <- rownames(subset(de, log2FoldChange>0 & padj<0.05 ))
dn <- rownames(subset(de, log2FoldChange<0 & padj<0.05 ))
str(up)
str(dn)
# MA plot
sig <-subset(de, padj < 0.05 )
GENESUP <- length(up)
GENESDN <- length(dn)
SUBHEADER = paste(GENESUP, "up, ", GENESDN, "down")
ns <-subset(de, padj > 0.05 )
plot(log2(de$baseMean),de$log2FoldChange,
xlab="log2 basemean", ylab="log2 foldchange",
pch=19, cex=0.5, col="dark gray",
main="smear plot")
points(log2(sig$baseMean),sig$log2FoldChange,
pch=19, cex=0.5, col="red")
mtext(SUBHEADER)
# heatmap
yn <- y/colSums(y)*1000000
yf <- yn[which(rownames(yn) %in% rownames(de)[1:50]),]
mycols <- gsub("0","yellow",ss$trt)
mycols <- gsub("1","orange",mycols)
colfunc <- colorRampPalette(c("blue", "white", "red"))
heatmap.2( as.matrix(yf), col=colfunc(25),scale="row",
ColSideColors =mycols ,trace="none",
margin = c(10,10), cexRow=0.6, cexCol=0.8 , main="Top 50 genes by p-val")
mtext("yellow=ctrl, orange=trt")
return(de)
}Run differential expression
Here execute the analysis.
pc1de <- run_de(sspc1,pc1)## converting counts to integer mode## the design formula contains one or more numeric variables with integer values,
## specifying a model with increasing fold change for higher values.
## did you mean for this to be a factor? if so, first convert
## this variable to a factor using the factor() function## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing
## chr [1:3097] "PHYCI_100005" "PHYCI_148089" "PHYCI_178772" "PHYCI_261186" ...
## chr [1:5631] "PHYCI_121119" "PHYCI_29764" "PHYCI_273086" "PHYCI_105396" ...
head(as.data.frame(pc1de),20)## baseMean log2FoldChange lfcSE stat pvalue
## PHYCI_100005 7458.5689 10.997984 0.2230196 49.31397 0.000000e+00
## PHYCI_148089 4240.2584 8.645469 0.1582754 54.62296 0.000000e+00
## PHYCI_178772 1454.0794 10.493558 0.2608642 40.22613 0.000000e+00
## PHYCI_261186 1638.8110 9.831829 0.2595520 37.88000 0.000000e+00
## PHYCI_304292 2032.5337 5.963319 0.1285101 46.40349 0.000000e+00
## PHYCI_311759 502.6617 7.419734 0.1824379 40.66992 0.000000e+00
## PHYCI_312248 946.1004 8.014527 0.1838346 43.59641 0.000000e+00
## PHYCI_314066 5082.9059 9.221590 0.2080240 44.32946 0.000000e+00
## PHYCI_94709 2090.7740 3.304664 0.0790668 41.79585 0.000000e+00
## PHYCI_99054 1491.6292 10.570734 0.2782320 37.99251 0.000000e+00
## PHYCI_462950 1079.1206 4.358389 0.1173852 37.12895 9.587183e-302
## PHYCI_226119 991.9203 8.536650 0.2343888 36.42090 1.987581e-290
## PHYCI_99743 1269.6615 7.457800 0.2052497 36.33525 4.493054e-289
## PHYCI_82350 432.5156 6.879233 0.1901251 36.18267 1.140430e-286
## PHYCI_96449 317.9925 7.602397 0.2137732 35.56291 5.247620e-277
## PHYCI_129006 1232.7239 7.581010 0.2159319 35.10833 5.030335e-270
## PHYCI_90225 496.1246 4.510982 0.1288738 35.00310 2.018007e-268
## PHYCI_66467 545.5239 4.614214 0.1323926 34.85251 3.901331e-266
## PHYCI_90480 495.5453 8.490696 0.2460764 34.50431 6.911825e-261
## PHYCI_32109 1870.3312 9.162297 0.2740070 33.43818 3.822824e-245
## padj
## PHYCI_100005 0.000000e+00
## PHYCI_148089 0.000000e+00
## PHYCI_178772 0.000000e+00
## PHYCI_261186 0.000000e+00
## PHYCI_304292 0.000000e+00
## PHYCI_311759 0.000000e+00
## PHYCI_312248 0.000000e+00
## PHYCI_314066 0.000000e+00
## PHYCI_94709 0.000000e+00
## PHYCI_99054 0.000000e+00
## PHYCI_462950 1.184801e-298
## PHYCI_226119 2.251598e-287
## PHYCI_99743 4.698352e-286
## PHYCI_82350 1.107358e-283
## PHYCI_96449 4.755743e-274
## PHYCI_129006 4.273898e-267
## PHYCI_90225 1.613694e-265
## PHYCI_66467 2.946372e-263
## PHYCI_90480 4.945229e-258
## PHYCI_32109 2.598374e-242write.table(pc1de,file="pc1de.tsv",quote=FALSE,sep="\t")
pc2de <- run_de(sspc2,pc2)## converting counts to integer mode## the design formula contains one or more numeric variables with integer values,
## specifying a model with increasing fold change for higher values.
## did you mean for this to be a factor? if so, first convert
## this variable to a factor using the factor() function## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing
## chr [1:4431] "PHYCI_148089" "PHYCI_304292" "PHYCI_562950" "PHYCI_83306" ...
## chr [1:5982] "PHYCI_75824" "PHYCI_22350" "PHYCI_98777" "PHYCI_98676" ...
head(as.data.frame(pc2de),20)## baseMean log2FoldChange lfcSE stat pvalue
## PHYCI_148089 7311.4238 8.098697 0.1975795 40.98957 0.000000e+00
## PHYCI_304292 17192.2704 7.712978 0.1581552 48.76843 0.000000e+00
## PHYCI_562950 1146.5213 5.193429 0.1364593 38.05846 0.000000e+00
## PHYCI_83306 1240.3340 5.722536 0.1299755 44.02781 0.000000e+00
## PHYCI_92107 990.1029 5.065277 0.1299843 38.96837 0.000000e+00
## PHYCI_96416 906.9505 7.480746 0.1911710 39.13117 0.000000e+00
## PHYCI_314066 9749.7982 8.827690 0.2467176 35.78055 2.216680e-280
## PHYCI_178772 2806.5556 10.111498 0.2978683 33.94620 1.387911e-252
## PHYCI_82350 1020.2754 6.739978 0.2028019 33.23429 3.442681e-242
## PHYCI_94850 2009.4592 8.648859 0.2603423 33.22111 5.336546e-242
## PHYCI_314209 579.6973 6.314685 0.1908979 33.07885 5.987175e-240
## PHYCI_95536 1043.6801 8.607659 0.2616344 32.89957 2.229401e-237
## PHYCI_117702 4753.9360 11.341791 0.3464195 32.74005 4.206719e-235
## PHYCI_226119 3844.5930 9.152973 0.2839611 32.23319 6.051493e-228
## PHYCI_195336 2104.5977 5.724698 0.1777601 32.20463 1.520367e-227
## PHYCI_312248 2354.9656 8.002739 0.2496927 32.05036 2.170423e-225
## PHYCI_93493 2747.5803 10.417656 0.3251041 32.04406 2.655837e-225
## PHYCI_261186 3537.7954 9.610484 0.3022085 31.80083 6.303262e-222
## PHYCI_77017 405.0280 6.536352 0.2080040 31.42417 9.463696e-217
## PHYCI_129006 4165.6141 8.003471 0.2555795 31.31500 2.916188e-215
## padj
## PHYCI_148089 0.000000e+00
## PHYCI_304292 0.000000e+00
## PHYCI_562950 0.000000e+00
## PHYCI_83306 0.000000e+00
## PHYCI_92107 0.000000e+00
## PHYCI_96416 0.000000e+00
## PHYCI_314066 4.404227e-277
## PHYCI_178772 2.412884e-249
## PHYCI_82350 5.320090e-239
## PHYCI_94850 7.422068e-239
## PHYCI_314209 7.569966e-237
## PHYCI_95536 2.583876e-234
## PHYCI_117702 4.500543e-232
## PHYCI_226119 6.011726e-225
## PHYCI_195336 1.409684e-224
## PHYCI_312248 1.886640e-222
## PHYCI_93493 2.172787e-222
## PHYCI_261186 4.870320e-219
## PHYCI_77017 6.927426e-214
## PHYCI_129006 2.027917e-212write.table(pc2de,file="pc2de.tsv",quote=FALSE,sep="\t")
pc3de <- run_de(sspc3,pc3)## converting counts to integer mode## the design formula contains one or more numeric variables with integer values,
## specifying a model with increasing fold change for higher values.
## did you mean for this to be a factor? if so, first convert
## this variable to a factor using the factor() function## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing
## chr [1:4418] "PHYCI_117702" "PHYCI_175455" "PHYCI_211022" "PHYCI_211711" ...
## chr [1:4319] "PHYCI_93190" "PHYCI_93788" "PHYCI_120595" "PHYCI_75824" ...
head(as.data.frame(pc3de),20)## baseMean log2FoldChange lfcSE stat pvalue
## PHYCI_117702 7180.3208 11.196418 0.2908850 38.49088 0.000000e+00
## PHYCI_175455 1550.1181 7.181042 0.1751215 41.00605 0.000000e+00
## PHYCI_211022 21589.4507 6.186448 0.1490174 41.51493 0.000000e+00
## PHYCI_211711 17551.6039 12.018585 0.3173072 37.87681 0.000000e+00
## PHYCI_304292 22534.1708 7.365443 0.1249478 58.94816 0.000000e+00
## PHYCI_91208 4977.9701 8.423477 0.2100068 40.11049 0.000000e+00
## PHYCI_270176 2889.9084 5.744200 0.1596402 35.98216 1.590633e-283
## PHYCI_93271 14277.8039 11.740598 0.3282673 35.76536 3.818083e-280
## PHYCI_221789 28790.4124 8.219553 0.2351920 34.94826 1.375939e-267
## PHYCI_93493 2840.4293 9.725936 0.2900585 33.53095 1.706616e-246
## PHYCI_98096 1173.9039 3.649188 0.1094984 33.32642 1.600183e-243
## PHYCI_94539 962.4735 6.449695 0.1948639 33.09845 3.128135e-240
## PHYCI_314066 6649.5784 7.534062 0.2327583 32.36860 7.594337e-230
## PHYCI_127817 7207.8895 8.141552 0.2569893 31.68051 2.883282e-220
## PHYCI_12045 6647.9654 12.566207 0.3972484 31.63312 1.294193e-219
## PHYCI_105524 593.0126 6.417022 0.2040211 31.45274 3.850613e-217
## PHYCI_99949 1230.6371 8.509814 0.2792963 30.46877 6.758654e-204
## PHYCI_20242 1610.8834 4.671332 0.1536513 30.40216 5.143525e-203
## PHYCI_95536 1413.7804 8.305252 0.2737516 30.33864 3.547774e-202
## PHYCI_312248 1607.3979 6.709817 0.2252223 29.79198 4.962244e-195
## padj
## PHYCI_117702 0.000000e+00
## PHYCI_175455 0.000000e+00
## PHYCI_211022 0.000000e+00
## PHYCI_211711 0.000000e+00
## PHYCI_304292 0.000000e+00
## PHYCI_91208 0.000000e+00
## PHYCI_270176 3.207397e-280
## PHYCI_93271 6.736531e-277
## PHYCI_221789 2.157931e-264
## PHYCI_93493 2.408889e-243
## PHYCI_98096 2.053325e-240
## PHYCI_94539 3.679468e-237
## PHYCI_314066 8.245698e-227
## PHYCI_127817 2.906966e-217
## PHYCI_12045 1.217835e-216
## PHYCI_105524 3.396962e-214
## PHYCI_99949 5.611671e-201
## PHYCI_20242 4.033381e-200
## PHYCI_95536 2.635623e-199
## PHYCI_312248 3.502104e-192write.table(pc3de,file="pc3de.tsv",quote=FALSE,sep="\t")
pc4de <- run_de(sspc4,pc4)## converting counts to integer mode## the design formula contains one or more numeric variables with integer values,
## specifying a model with increasing fold change for higher values.
## did you mean for this to be a factor? if so, first convert
## this variable to a factor using the factor() function## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing
## chr [1:1326] "PHYCI_12045" "PHYCI_117702" "PHYCI_92137" "PHYCI_90358" ...
## chr [1:1086] "PHYCI_223483" "PHYCI_15784" "PHYCI_96550" "PHYCI_18201" ...
head(as.data.frame(pc4de),20)## baseMean log2FoldChange lfcSE stat pvalue
## PHYCI_12045 493.01890 3.467028 0.2503337 13.849625 1.278781e-43
## PHYCI_117702 500.40890 3.505902 0.2932693 11.954549 6.146648e-33
## PHYCI_92137 381.19590 2.179718 0.1824115 11.949459 6.534912e-33
## PHYCI_90358 306.12263 3.125956 0.2619359 11.934053 7.864928e-33
## PHYCI_90348 241.21394 3.341102 0.2824812 11.827696 2.807428e-32
## PHYCI_10098 177.57117 4.038288 0.3422626 11.798802 3.959073e-32
## PHYCI_273086 214.14227 3.597291 0.3305310 10.883369 1.383536e-27
## PHYCI_89525 451.95737 3.208252 0.2973471 10.789585 3.855282e-27
## PHYCI_93841 478.20968 3.708064 0.3464435 10.703229 9.829174e-27
## PHYCI_108488 1378.76622 8.297350 0.7860479 10.555781 4.776520e-26
## PHYCI_158482 224.21351 5.143427 0.4882836 10.533687 6.042062e-26
## PHYCI_223483 441.68792 -1.509241 0.1447522 -10.426375 1.879210e-25
## PHYCI_308281 1235.11526 1.599409 0.1538290 10.397319 2.550083e-25
## PHYCI_15784 47.82216 -5.195938 0.5224053 -9.946180 2.620499e-23
## PHYCI_89007 159.31340 3.820137 0.3916107 9.754936 1.757118e-22
## PHYCI_96257 178.93061 3.504188 0.3646497 9.609738 7.273353e-22
## PHYCI_93271 385.17698 2.990896 0.3133651 9.544447 1.368351e-21
## PHYCI_95308 427.23652 1.824879 0.1931011 9.450380 3.376014e-21
## PHYCI_96550 97.64992 -2.436131 0.2626229 -9.276154 1.757061e-20
## PHYCI_304292 2257.62318 1.484546 0.1603177 9.260029 2.043765e-20
## padj
## PHYCI_12045 1.173793e-39
## PHYCI_117702 1.804804e-29
## PHYCI_92137 1.804804e-29
## PHYCI_90358 1.804804e-29
## PHYCI_90348 5.153877e-29
## PHYCI_10098 6.056722e-29
## PHYCI_273086 1.814211e-24
## PHYCI_89525 4.423455e-24
## PHYCI_93841 1.002467e-23
## PHYCI_108488 4.384368e-23
## PHYCI_158482 5.041827e-23
## PHYCI_223483 1.437439e-22
## PHYCI_308281 1.800555e-22
## PHYCI_15784 1.718111e-20
## PHYCI_89007 1.075239e-19
## PHYCI_96257 4.172632e-19
## PHYCI_93271 7.388290e-19
## PHYCI_95308 1.721579e-18
## PHYCI_96550 8.488454e-18
## PHYCI_304292 9.379862e-18write.table(pc4de,file="pc4de.tsv",quote=FALSE,sep="\t")
pc5de <- run_de(sspc5,pc5)## converting counts to integer mode## the design formula contains one or more numeric variables with integer values,
## specifying a model with increasing fold change for higher values.
## did you mean for this to be a factor? if so, first convert
## this variable to a factor using the factor() function## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing
## chr [1:3053] "PHYCI_8031" "PHYCI_233579" "PHYCI_93203" "PHYCI_118265" ...
## chr [1:2648] "PHYCI_12583" "PHYCI_475096" "PHYCI_92428" "PHYCI_120131" ...
head(as.data.frame(pc5de),20)## baseMean log2FoldChange lfcSE stat pvalue
## PHYCI_12583 415.99068 -2.935618 0.2025090 -14.49623 1.279799e-47
## PHYCI_8031 964.85469 5.040191 0.3584301 14.06185 6.515650e-45
## PHYCI_475096 751.84974 -2.927662 0.2107154 -13.89392 6.895967e-44
## PHYCI_233579 415.16751 3.451813 0.2535534 13.61375 3.317605e-42
## PHYCI_93203 12975.14164 4.812676 0.3807587 12.63970 1.275349e-36
## PHYCI_118265 489.31697 2.989477 0.2427597 12.31455 7.563757e-35
## PHYCI_96111 326.61098 3.719613 0.3054987 12.17555 4.195803e-34
## PHYCI_92428 4555.33019 -2.678746 0.2244013 -11.93730 7.563682e-33
## PHYCI_208360 1154.52945 2.613164 0.2191556 11.92378 8.897718e-33
## PHYCI_230688 252.71981 6.081632 0.5114518 11.89092 1.319452e-32
## PHYCI_99495 385.28074 3.508029 0.2964809 11.83222 2.659999e-32
## PHYCI_77249 1126.84057 4.407950 0.3754663 11.73993 7.954898e-32
## PHYCI_120131 351.52736 -2.467330 0.2121501 -11.63011 2.897086e-31
## PHYCI_248382 226.35334 4.587882 0.3949973 11.61497 3.459104e-31
## PHYCI_397273 231.12373 4.266243 0.3692748 11.55303 7.126240e-31
## PHYCI_91722 303.81923 3.083876 0.2676523 11.52195 1.022692e-30
## PHYCI_261804 2708.65621 -2.631025 0.2353801 -11.17777 5.238886e-29
## PHYCI_437341 72.28638 -5.520404 0.4942563 -11.16911 5.775604e-29
## PHYCI_148089 3709.16170 -3.073780 0.2754378 -11.15961 6.427071e-29
## PHYCI_22350 162.73714 4.803891 0.4377361 10.97440 5.074204e-28
## padj
## PHYCI_12583 1.511059e-43
## PHYCI_8031 3.846514e-41
## PHYCI_475096 2.714023e-40
## PHYCI_233579 9.792741e-39
## PHYCI_93203 3.011610e-33
## PHYCI_118265 1.488421e-31
## PHYCI_96111 7.077121e-31
## PHYCI_92428 1.116305e-29
## PHYCI_208360 1.167282e-29
## PHYCI_230688 1.557877e-29
## PHYCI_99495 2.855147e-29
## PHYCI_77249 7.826956e-29
## PHYCI_120131 2.631223e-28
## PHYCI_248382 2.917260e-28
## PHYCI_397273 5.609301e-28
## PHYCI_91722 7.546825e-28
## PHYCI_261804 3.638560e-26
## PHYCI_437341 3.788475e-26
## PHYCI_148089 3.993917e-26
## PHYCI_22350 2.995556e-25write.table(pc5de,file="pc5de.tsv",quote=FALSE,sep="\t")
la1de <- run_de(ssla1,la1)## converting counts to integer mode## the design formula contains one or more numeric variables with integer values,
## specifying a model with increasing fold change for higher values.
## did you mean for this to be a factor? if so, first convert
## this variable to a factor using the factor() function## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing
## chr [1:6518] "TanjilG_25668" "TanjilG_13054" "TanjilG_31530" ...
## chr [1:7034] "TanjilG_21738" "TanjilG_30005" "TanjilG_03097" ...
head(as.data.frame(la1de),20)## baseMean log2FoldChange lfcSE stat pvalue
## TanjilG_21738 6084.0947 -5.652528 0.13930749 -40.57591 0.000000e+00
## TanjilG_25668 3933.4631 3.767639 0.10573880 35.63157 4.547324e-278
## TanjilG_30005 1200.9058 -6.397750 0.19248660 -33.23738 3.106683e-242
## TanjilG_03097 1621.5440 -5.792440 0.17643279 -32.83086 2.137015e-236
## TanjilG_08363 1845.9082 -3.806889 0.12168259 -31.28540 7.371131e-215
## TanjilG_03248 21389.4434 -7.396194 0.24433783 -30.27036 2.815846e-201
## TanjilG_10617 1120.1031 -5.145379 0.17660238 -29.13539 1.279558e-186
## TanjilG_26762 2768.7053 -8.161854 0.28742482 -28.39648 2.234732e-177
## TanjilG_13054 9628.9486 5.093092 0.17972001 28.33903 1.142542e-176
## TanjilG_06910 1323.2013 -3.372833 0.12116335 -27.83708 1.544363e-170
## TanjilG_25621 4335.8706 -2.137132 0.07884346 -27.10601 8.363966e-162
## TanjilG_31530 8651.7545 4.561844 0.16840806 27.08804 1.362170e-161
## TanjilG_01633 1052.7061 -5.487506 0.20437896 -26.84966 8.511093e-159
## TanjilG_00341 867.8256 5.108753 0.19029955 26.84585 9.429831e-159
## TanjilG_30454 1151.2735 5.869116 0.21960312 26.72601 2.346905e-157
## TanjilG_23160 2955.5919 -4.705178 0.17618628 -26.70570 4.041517e-157
## TanjilG_14202 2118.6838 -5.539758 0.21267236 -26.04832 1.405517e-149
## TanjilG_06104 575.4518 -5.522615 0.21272651 -25.96110 1.362457e-148
## TanjilG_27770 2787.5733 -3.449582 0.13382908 -25.77602 1.647014e-146
## TanjilG_31035 12137.9734 -4.078743 0.15880544 -25.68390 1.768797e-145
## padj
## TanjilG_21738 0.000000e+00
## TanjilG_25668 5.791927e-274
## TanjilG_30005 2.637988e-238
## TanjilG_03097 1.360958e-232
## TanjilG_08363 3.755444e-211
## TanjilG_03248 1.195514e-197
## TanjilG_10617 4.656495e-183
## TanjilG_26762 7.115945e-174
## TanjilG_13054 3.233902e-173
## TanjilG_06910 3.934109e-167
## TanjilG_25621 1.936942e-158
## TanjilG_31530 2.891661e-158
## TanjilG_01633 1.667781e-155
## TanjilG_00341 1.715825e-155
## TanjilG_30454 3.985671e-154
## TanjilG_23160 6.434600e-154
## TanjilG_14202 2.106126e-146
## TanjilG_06104 1.928180e-145
## TanjilG_27770 2.208212e-143
## TanjilG_31035 2.252916e-142write.table(la1de,file="la1de.tsv",quote=FALSE,sep="\t")
la2de <- run_de(ssla2,la2)## converting counts to integer mode## the design formula contains one or more numeric variables with integer values,
## specifying a model with increasing fold change for higher values.
## did you mean for this to be a factor? if so, first convert
## this variable to a factor using the factor() function## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing
## chr [1:7804] "TanjilG_31530" "TanjilG_32339" "TanjilG_22068" ...
## chr [1:7822] "TanjilG_18821" "TanjilG_00128" "TanjilG_09966" ...
head(as.data.frame(la2de),20)## baseMean log2FoldChange lfcSE stat pvalue
## TanjilG_31530 18831.049 5.970301 0.1183801 50.43331 0.000000e+00
## TanjilG_32339 3137.729 5.036788 0.1503985 33.48960 6.829453e-246
## TanjilG_22068 3672.534 8.102453 0.2528472 32.04487 2.588392e-225
## TanjilG_13054 24324.690 6.705454 0.2100289 31.92633 1.151295e-223
## TanjilG_13015 4361.562 4.665514 0.1475629 31.61712 2.148236e-219
## TanjilG_17787 1569.649 7.148794 0.2267476 31.52753 3.644347e-218
## TanjilG_12890 2251.805 6.526765 0.2113502 30.88128 2.131034e-209
## TanjilG_28363 7500.044 6.952005 0.2328994 29.84981 8.827638e-196
## TanjilG_00065 2935.895 3.472683 0.1190245 29.17621 3.886634e-187
## TanjilG_21470 4248.520 4.295556 0.1492404 28.78279 3.522390e-182
## TanjilG_18821 2291.981 -4.777481 0.1662683 -28.73356 1.453791e-181
## TanjilG_17822 15370.904 7.411157 0.2609336 28.40246 1.885294e-177
## TanjilG_00128 2405.552 -7.394093 0.2606243 -28.37070 4.649857e-177
## TanjilG_14085 4546.297 4.318613 0.1541708 28.01187 1.164606e-172
## TanjilG_09966 6826.986 -4.020823 0.1442504 -27.87391 5.528228e-171
## TanjilG_17902 26138.068 3.882039 0.1401352 27.70209 6.588497e-169
## TanjilG_07217 6522.758 3.879310 0.1407028 27.57094 2.482738e-167
## TanjilG_01075 1772.947 5.692016 0.2067276 27.53390 6.899291e-167
## TanjilG_20161 1733.517 7.952862 0.2897401 27.44827 7.286875e-166
## TanjilG_06910 1053.586 -4.714447 0.1735662 -27.16224 1.815165e-162
## padj
## TanjilG_31530 0.000000e+00
## TanjilG_32339 8.644038e-242
## TanjilG_22068 2.184085e-221
## TanjilG_13054 7.285969e-220
## TanjilG_13015 1.087609e-215
## TanjilG_17787 1.537550e-214
## TanjilG_12890 7.706429e-206
## TanjilG_28363 2.793285e-192
## TanjilG_00065 1.093181e-183
## TanjilG_21470 8.916579e-179
## TanjilG_18821 3.345569e-178
## TanjilG_17822 3.977028e-174
## TanjilG_00128 9.054345e-174
## TanjilG_14085 2.105774e-169
## TanjilG_09966 9.329437e-168
## TanjilG_17902 1.042383e-165
## TanjilG_07217 3.696944e-164
## TanjilG_01075 9.702703e-164
## TanjilG_20161 9.708419e-163
## TanjilG_06910 2.297454e-159write.table(la2de,file="la2de.tsv",quote=FALSE,sep="\t")
la3de <- run_de(ssla3,la3)## converting counts to integer mode## the design formula contains one or more numeric variables with integer values,
## specifying a model with increasing fold change for higher values.
## did you mean for this to be a factor? if so, first convert
## this variable to a factor using the factor() function## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing
## chr [1:9144] "TanjilG_14085" "TanjilG_22068" "TanjilG_32339" ...
## chr [1:9822] "TanjilG_05765" "TanjilG_03248" "TanjilG_13344" ...
head(as.data.frame(la3de),20)## baseMean log2FoldChange lfcSE stat pvalue
## TanjilG_14085 10544.451 6.012654 0.1506101 39.92199 0.000000e+00
## TanjilG_22068 11968.768 10.242133 0.2684224 38.15677 0.000000e+00
## TanjilG_32339 5098.222 6.191774 0.1657555 37.35486 2.114549e-305
## TanjilG_04385 12639.966 7.233356 0.1951585 37.06400 1.068503e-300
## TanjilG_15363 17074.265 10.933488 0.2950412 37.05750 1.359810e-300
## TanjilG_05765 5217.183 -5.403651 0.1458608 -37.04663 2.034852e-300
## TanjilG_10317 2383.926 6.416327 0.1754468 36.57135 8.165695e-293
## TanjilG_13823 19785.918 5.507165 0.1549989 35.53036 1.670849e-276
## TanjilG_13054 34525.517 7.649860 0.2153720 35.51929 2.476631e-276
## TanjilG_11778 3028.944 8.226423 0.2357819 34.88996 1.055538e-266
## TanjilG_07217 7116.667 4.467228 0.1313454 34.01130 1.516531e-253
## TanjilG_17902 37223.668 4.869821 0.1438508 33.85329 3.246496e-251
## TanjilG_00940 2665.445 8.374035 0.2475449 33.82834 7.557868e-251
## TanjilG_31530 40299.460 7.513685 0.2227741 33.72783 2.260379e-249
## TanjilG_25612 11144.982 6.413451 0.1913927 33.50938 3.519062e-246
## TanjilG_03248 13282.255 -10.356194 0.3128179 -33.10614 2.424699e-240
## TanjilG_28363 13185.665 8.203310 0.2489496 32.95169 4.000586e-238
## TanjilG_13842 6944.686 4.756833 0.1454892 32.69545 1.812743e-234
## TanjilG_32922 6552.569 5.109301 0.1564465 32.65845 6.079320e-234
## TanjilG_13015 7087.619 5.828564 0.1787085 32.61493 2.519471e-233
## padj
## TanjilG_14085 0.000000e+00
## TanjilG_22068 0.000000e+00
## TanjilG_32339 1.766917e-301
## TanjilG_04385 6.696309e-297
## TanjilG_15363 6.817546e-297
## TanjilG_05765 8.501613e-297
## TanjilG_10317 2.924252e-289
## TanjilG_13823 5.235606e-273
## TanjilG_13054 6.898242e-273
## TanjilG_11778 2.646022e-263
## TanjilG_07217 3.456036e-250
## TanjilG_17902 6.781931e-248
## TanjilG_00940 1.457389e-247
## TanjilG_31530 4.047370e-246
## TanjilG_25612 5.881056e-243
## TanjilG_03248 3.798897e-237
## TanjilG_28363 5.899217e-235
## TanjilG_13842 2.524547e-231
## TanjilG_32922 8.020863e-231
## TanjilG_13015 3.157906e-230write.table(la3de,file="la3de.tsv",quote=FALSE,sep="\t")
la4de <- run_de(ssla4,la4)## converting counts to integer mode## the design formula contains one or more numeric variables with integer values,
## specifying a model with increasing fold change for higher values.
## did you mean for this to be a factor? if so, first convert
## this variable to a factor using the factor() function## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing
## chr [1:6150] "TanjilG_24871" "TanjilG_32339" "TanjilG_21471" ...
## chr [1:5356] "TanjilG_18821" "TanjilG_29255" "TanjilG_28251" ...
head(as.data.frame(la4de),20)## baseMean log2FoldChange lfcSE stat pvalue
## TanjilG_18821 958.6247 -3.368914 0.13902077 -24.23317 9.951092e-130
## TanjilG_24871 1605.7445 1.872117 0.08262968 22.65672 1.197838e-113
## TanjilG_29255 2061.9297 -3.455047 0.15263411 -22.63614 1.910723e-113
## TanjilG_32339 3883.3694 2.217731 0.10549650 21.02185 4.139881e-98
## TanjilG_21471 308.1423 2.863579 0.15231568 18.80029 7.510458e-79
## TanjilG_28744 3799.6755 1.635618 0.08742142 18.70958 4.136511e-78
## TanjilG_25070 1107.9430 4.176298 0.22322301 18.70908 4.175010e-78
## TanjilG_13779 346.0633 5.484429 0.30645959 17.89609 1.264899e-71
## TanjilG_05531 1692.0481 3.534705 0.20144793 17.54649 6.326910e-69
## TanjilG_28251 3274.8345 -1.853097 0.10710549 -17.30161 4.574352e-67
## TanjilG_31962 2778.0230 2.516025 0.14721632 17.09066 1.741798e-65
## TanjilG_02661 3432.4887 3.347197 0.19953694 16.77482 3.729585e-63
## TanjilG_01449 888.5795 3.454625 0.21214527 16.28424 1.277034e-59
## TanjilG_06979 6173.4875 -1.259014 0.07792704 -16.15631 1.025124e-58
## TanjilG_09968 3110.7802 -1.244185 0.07747173 -16.05986 4.876232e-58
## TanjilG_03391 7387.1028 1.938485 0.12166806 15.93257 3.765546e-57
## TanjilG_06485 1715.8228 4.857889 0.30921939 15.71017 1.288362e-55
## TanjilG_05088 1458.4024 -4.076460 0.25997392 -15.68027 2.063853e-55
## TanjilG_21567 909.0565 -2.522357 0.16296146 -15.47824 4.866072e-54
## TanjilG_07959 968.4502 1.627990 0.10633481 15.31004 6.552226e-53
## padj
## TanjilG_18821 2.485484e-125
## TanjilG_24871 1.495920e-109
## TanjilG_29255 1.590804e-109
## TanjilG_32339 2.585045e-94
## TanjilG_21471 3.751774e-75
## TanjilG_28744 1.489703e-74
## TanjilG_25070 1.489703e-74
## TanjilG_13779 3.949173e-68
## TanjilG_05531 1.755858e-65
## TanjilG_28251 1.142536e-63
## TanjilG_31962 3.954989e-62
## TanjilG_02661 7.762821e-60
## TanjilG_01449 2.453576e-56
## TanjilG_06979 1.828895e-55
## TanjilG_09968 8.119577e-55
## TanjilG_03391 5.878253e-54
## TanjilG_06485 1.892907e-52
## TanjilG_05088 2.863825e-52
## TanjilG_21567 6.396836e-51
## TanjilG_07959 8.182748e-50write.table(la4de,file="la4de.tsv",quote=FALSE,sep="\t")
la5de <- run_de(ssla5,la5)## converting counts to integer mode## the design formula contains one or more numeric variables with integer values,
## specifying a model with increasing fold change for higher values.
## did you mean for this to be a factor? if so, first convert
## this variable to a factor using the factor() function## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing
## chr [1:6150] "TanjilG_24871" "TanjilG_32339" "TanjilG_21471" ...
## chr [1:5356] "TanjilG_18821" "TanjilG_29255" "TanjilG_28251" ...
head(as.data.frame(la5de),20)## baseMean log2FoldChange lfcSE stat pvalue
## TanjilG_18821 958.6247 -3.368914 0.13902077 -24.23317 9.951092e-130
## TanjilG_24871 1605.7445 1.872117 0.08262968 22.65672 1.197838e-113
## TanjilG_29255 2061.9297 -3.455047 0.15263411 -22.63614 1.910723e-113
## TanjilG_32339 3883.3694 2.217731 0.10549650 21.02185 4.139881e-98
## TanjilG_21471 308.1423 2.863579 0.15231568 18.80029 7.510458e-79
## TanjilG_28744 3799.6755 1.635618 0.08742142 18.70958 4.136511e-78
## TanjilG_25070 1107.9430 4.176298 0.22322301 18.70908 4.175010e-78
## TanjilG_13779 346.0633 5.484429 0.30645959 17.89609 1.264899e-71
## TanjilG_05531 1692.0481 3.534705 0.20144793 17.54649 6.326910e-69
## TanjilG_28251 3274.8345 -1.853097 0.10710549 -17.30161 4.574352e-67
## TanjilG_31962 2778.0230 2.516025 0.14721632 17.09066 1.741798e-65
## TanjilG_02661 3432.4887 3.347197 0.19953694 16.77482 3.729585e-63
## TanjilG_01449 888.5795 3.454625 0.21214527 16.28424 1.277034e-59
## TanjilG_06979 6173.4875 -1.259014 0.07792704 -16.15631 1.025124e-58
## TanjilG_09968 3110.7802 -1.244185 0.07747173 -16.05986 4.876232e-58
## TanjilG_03391 7387.1028 1.938485 0.12166806 15.93257 3.765546e-57
## TanjilG_06485 1715.8228 4.857889 0.30921939 15.71017 1.288362e-55
## TanjilG_05088 1458.4024 -4.076460 0.25997392 -15.68027 2.063853e-55
## TanjilG_21567 909.0565 -2.522357 0.16296146 -15.47824 4.866072e-54
## TanjilG_07959 968.4502 1.627990 0.10633481 15.31004 6.552226e-53
## padj
## TanjilG_18821 2.485484e-125
## TanjilG_24871 1.495920e-109
## TanjilG_29255 1.590804e-109
## TanjilG_32339 2.585045e-94
## TanjilG_21471 3.751774e-75
## TanjilG_28744 1.489703e-74
## TanjilG_25070 1.489703e-74
## TanjilG_13779 3.949173e-68
## TanjilG_05531 1.755858e-65
## TanjilG_28251 1.142536e-63
## TanjilG_31962 3.954989e-62
## TanjilG_02661 7.762821e-60
## TanjilG_01449 2.453576e-56
## TanjilG_06979 1.828895e-55
## TanjilG_09968 8.119577e-55
## TanjilG_03391 5.878253e-54
## TanjilG_06485 1.892907e-52
## TanjilG_05088 2.863825e-52
## TanjilG_21567 6.396836e-51
## TanjilG_07959 8.182748e-50write.table(la5de,file="la5de.tsv",quote=FALSE,sep="\t")Session Information
sessionInfo()## R version 4.0.1 (2020-06-06)
## Platform: x86_64-conda_cos6-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.1 LTS
##
## Matrix products: default
## BLAS/LAPACK: /ceph-g/opt/miniconda3/envs/R40/lib/libopenblasp-r0.3.9.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] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] gplots_3.1.0 Biostrings_2.56.0
## [3] XVector_0.28.0 reshape2_1.4.4
## [5] mitch_1.0.10 DESeq2_1.28.1
## [7] SummarizedExperiment_1.18.2 DelayedArray_0.14.1
## [9] matrixStats_0.57.0 Biobase_2.48.0
## [11] GenomicRanges_1.40.0 GenomeInfoDb_1.24.2
## [13] IRanges_2.22.2 S4Vectors_0.26.1
## [15] BiocGenerics_0.34.0
##
## loaded via a namespace (and not attached):
## [1] bit64_4.0.5 splines_4.0.1 gtools_3.8.2
## [4] shiny_1.4.0.2 blob_1.2.1 GenomeInfoDbData_1.2.3
## [7] yaml_2.2.1 pillar_1.4.4 RSQLite_2.2.1
## [10] lattice_0.20-41 glue_1.4.1 digest_0.6.25
## [13] RColorBrewer_1.1-2 promises_1.1.1 colorspace_1.4-1
## [16] htmltools_0.5.0 httpuv_1.5.4 Matrix_1.2-18
## [19] plyr_1.8.6 XML_3.99-0.5 pkgconfig_2.0.3
## [22] genefilter_1.70.0 zlibbioc_1.34.0 purrr_0.3.4
## [25] xtable_1.8-4 scales_1.1.1 later_1.1.0.1
## [28] BiocParallel_1.22.0 tibble_3.0.1 annotate_1.66.0
## [31] echarts4r_0.3.3 generics_0.0.2 ggplot2_3.3.2
## [34] ellipsis_0.3.1 survival_3.2-3 magrittr_1.5
## [37] crayon_1.3.4 mime_0.9 evaluate_0.14
## [40] memoise_1.1.0 GGally_2.0.0 MASS_7.3-51.6
## [43] beeswarm_0.2.3 tools_4.0.1 lifecycle_0.2.0
## [46] stringr_1.4.0 munsell_0.5.0 locfit_1.5-9.4
## [49] AnnotationDbi_1.50.3 compiler_4.0.1 caTools_1.18.0
## [52] rlang_0.4.6 grid_4.0.1 RCurl_1.98-1.2
## [55] htmlwidgets_1.5.1 rmarkdown_2.3 bitops_1.0-6
## [58] gtable_0.3.0 DBI_1.1.0 reshape_0.8.8
## [61] R6_2.4.1 gridExtra_2.3 knitr_1.28
## [64] dplyr_1.0.0 fastmap_1.0.1 bit_4.0.4
## [67] KernSmooth_2.23-17 stringi_1.4.6 Rcpp_1.0.4.6
## [70] vctrs_0.3.1 geneplotter_1.66.0 tidyselect_1.1.0
## [73] xfun_0.15