require(edgeR)
# add gene names
gn=read.table("/data/results/reference/hg19/Homo_sapiens.GRCh38.84.gene-names.txt",sep="\t")
uniq=gn[!duplicated(gn$V2),]
rownames(uniq)=uniq$V2;


sample=factor(c(1,1,2,2,3,3))
condition=factor(c("UntrPMNs","PMNsIRAplasma","UntrPMNs","PMNsIRAplasma","UntrPMNs","PMNsIRAplasma"))
#condition=factor(c("UntrPMNs","PMNspolyP","UntrPMNs","PMNspolyP","UntrPMNs","PMNspolyP"))
design <- model.matrix(~sample+condition)
designUnPaired <- model.matrix(~condition)


#rnaseqMatrix1= read.table("gene_counts.txt",header=T,sep="\t",skip=1)
rnaseqMatrix1= read.table("/data/images/proton2/external/DrStakos/run409/gene_counts.txt",header=T,sep="\t",skip=1)

#Contrast: UntrPMN_vs_PMNsIRAplasma , pairs: (A vs B)
SFI=rnaseqMatrix1[,c(15,28,14,18,30,7)]
colnames(SFI)
y0 <- DGEList(counts=SFI, genes=rnaseqMatrix1[,1])
y0 <- estimateDisp(y0,design)

fit <- glmQLFit(y0, design)
qlf <- glmQLFTest(fit)
topTags(qlf)
tTags = topTags(qlf,n=NULL,adjust.method = "fdr");r <- rownames(topTags(qlf,n=NULL,adjust.method = "fdr")$table)
c=cpm(y0)[r, order(y0$samples$group)]
y=(tTags@.Data[[1]][1])
y <- data.frame(lapply(y, as.character), stringsAsFactors=FALSE) 
x=cbind(as.data.frame(tTags), uniq[y$genes,3],uniq[y$genes,4],c)
colnames(x)[7]="gene_name"
colnames(x)[8]="biotype"
#x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(qlf, p=0.05))
#write.table(tTags, file="edgeRgenes_UntrPMN_vs_PMNsIRAplasma_paired.csv",quote=F)
write.table(x, file="edgeRgenes_AvsB_UntrPMN_vs_PMNsIRAplasma_paired_featurecountExpr.csv",quote=F,row.names=F)


#Contrast: UntrPMN_vs_PMNsIRAplasmaTIC , pairs: (A vs C)

SFI=rnaseqMatrix1[,c(15,36,14,22,30,11)]
colnames(SFI)
y0 <- DGEList(counts=SFI, genes=rnaseqMatrix1[,1])
y0 <- estimateDisp(y0,design)
fit <- glmQLFit(y0, design)
qlf <- glmQLFTest(fit)
topTags(qlf)
tTags = topTags(qlf,n=NULL,adjust.method = "fdr");r <- rownames(topTags(qlf,n=NULL,adjust.method = "fdr")$table)
c=cpm(y0)[r, order(y0$samples$group)]
y=(tTags@.Data[[1]][1])
y <- data.frame(lapply(y, as.character), stringsAsFactors=FALSE) 
x=cbind(as.data.frame(tTags), uniq[y$genes,3],uniq[y$genes,4],c)
colnames(x)[7]="gene_name"
colnames(x)[8]="biotype"
#x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(qlf, p=0.05))
#write.table(tTags, file="edgeRgenes_UntrPMN_vs_PMNsIRAplasmaTIC_paired.csv",quote=FALSE)
write.table(x, file="edgeRgenes_AvsC_UntrPMN_vs_PMNsIRAplasmaTIC_paired_featurecountExpr.csv",quote=FALSE,row.names=F)


#Contrast: UntrPMN_vs_PMNsIRAplasmaIL29 , pairs: (A vs D)
SFI=rnaseqMatrix1[,c(15,10,14,26,30,19)]
colnames(SFI)
y0 <- DGEList(counts=SFI, genes=rnaseqMatrix1[,1])
y0 <- estimateDisp(y0,design)
y1 <- estimateDisp(y0,designUnPaired)
fit1 <- glmQLFit(y1, designUnPaired)
qlf1 <- glmQLFTest(fit1)
summary(de <- decideTestsDGE(qlf1, p=0.05))
topTags(qlf1,n=20)

fit <- glmQLFit(y0, design)
qlf <- glmQLFTest(fit)
topTags(qlf)
tTags = topTags(qlf,n=NULL,adjust.method = "fdr");r <- rownames(topTags(qlf,n=NULL,adjust.method = "fdr")$table)
c=cpm(y0)[r, order(y0$samples$group)]
y=(tTags@.Data[[1]][1])
y <- data.frame(lapply(y, as.character), stringsAsFactors=FALSE) 
x=cbind(as.data.frame(tTags), uniq[y$genes,3],uniq[y$genes,4],c)
colnames(x)[7]="gene_name"
colnames(x)[8]="biotype"
#x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(qlf, p=0.05))
#write.table(tTags, file="edgeRgenes_UntrPMN_vs_PMNsIRAplasmaIL29_paired.csv",quote=F)
write.table(x, file="edgeRgenes_AvsD_UntrPMN_vs_PMNsIRAplasmaIL29_paired_featurecountExpr.csv",quote=F,row.names=F)



#Contrast: PMNsIRAplasma_vs_PMNsIRAplasmaTIC , pairs: (B vs C)

SFI=rnaseqMatrix1[,c(28,36,18,22,7,11)]
colnames(SFI)
y0 <- DGEList(counts=SFI, genes=rnaseqMatrix1[,1])
y0 <- estimateDisp(y0,design)
fit <- glmQLFit(y0, design)
qlf <- glmQLFTest(fit)
topTags(qlf)
tTags = topTags(qlf,n=NULL,adjust.method = "fdr");r <- rownames(topTags(qlf,n=NULL,adjust.method = "fdr")$table)
c=cpm(y0)[r, order(y0$samples$group)]
y=(tTags@.Data[[1]][1])
y <- data.frame(lapply(y, as.character), stringsAsFactors=FALSE) 
x=cbind(as.data.frame(tTags), uniq[y$genes,3],uniq[y$genes,4],c)
colnames(x)[7]="gene_name"
colnames(x)[8]="biotype"
#x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(qlf, p=0.05))
#write.table(tTags, file="edgeRgenes_PMNsIRAplasma_vs_PMNsIRAplasmaTIC_paired.csv",quote=FALSE)
write.table(x, file="edgeRgenes_BvsC_PMNsIRAplasma_vs_PMNsIRAplasmaTIC_paired_featurecountExpr.csv",quote=FALSE,row.names=F)

#Contrast: PMNsIRAplasma_vs_PMNsIRAplasmaIL29 , pairs: (B vs D)
SFI=rnaseqMatrix1[,c(28,10,18,26,7,19)]
colnames(SFI)
y0 <- DGEList(counts=SFI, genes=rnaseqMatrix1[,1])
y0 <- estimateDisp(y0,design)
fit <- glmQLFit(y0, design)
qlf <- glmQLFTest(fit)
topTags(qlf)
tTags = topTags(qlf,n=NULL,adjust.method = "fdr");r <- rownames(topTags(qlf,n=NULL,adjust.method = "fdr")$table)
c=cpm(y0)[r, order(y0$samples$group)]
y=(tTags@.Data[[1]][1])
y <- data.frame(lapply(y, as.character), stringsAsFactors=FALSE) 
x=cbind(as.data.frame(tTags), uniq[y$genes,3],uniq[y$genes,4],c)
colnames(x)[7]="gene_name"
colnames(x)[8]="biotype"
#x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(qlf, p=0.05))
#write.table(tTags, file="edgeRgenes_PMNsIRAplasma_vs_PMNsIRAplasmaIL29_paired.csv",quote=FALSE)
write.table(x, file="edgeRgenes_BvsD_PMNsIRAplasma_vs_PMNsIRAplasmaIL29_paired_featurecountExpr.csv",quote=FALSE,row.names=F)



#EvsF_UntrPMNs_vs_PMNspolyP
SFI=rnaseqMatrix1[,c(23,27,16,20,13,17)]
colnames(SFI)
y0 <- DGEList(counts=SFI, genes=rnaseqMatrix1[,1])
y0 <- estimateDisp(y0,design)
fit <- glmQLFit(y0, design)
qlf <- glmQLFTest(fit)
topTags(qlf)
tTags = topTags(qlf,n=NULL,adjust.method = "fdr");r <- rownames(topTags(qlf,n=NULL,adjust.method = "fdr")$table)
c=cpm(y0)[r, order(y0$samples$group)]
y=(tTags@.Data[[1]][1])
y <- data.frame(lapply(y, as.character), stringsAsFactors=FALSE) 
x=cbind(as.data.frame(tTags), uniq[y$genes,3],uniq[y$genes,4],c)
colnames(x)[7]="gene_name"
colnames(x)[8]="biotype"
#x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(qlf, p=0.05))
#write.table(tTags, file="edgeRgenes_EvsF_UntrPMNs_vs_PMNspolyP_paired.csv",quote=FALSE)
write.table(x, file="edgeRgenes_EvsF_UntrPMNs_vs_PMNspolyP_paired_featurecountExpr.csv",quote=FALSE,row.names=F)



#EvsG_UntrPMNs_vs_PMNspolyPTIC
SFI=rnaseqMatrix1[,c(23,31,16,24,13,21)]
colnames(SFI)
y0 <- DGEList(counts=SFI, genes=rnaseqMatrix1[,1])
y0 <- estimateDisp(y0,design)
fit <- glmQLFit(y0, design)
qlf <- glmQLFTest(fit)
topTags(qlf)
tTags = topTags(qlf,n=NULL,adjust.method = "fdr");r <- rownames(topTags(qlf,n=NULL,adjust.method = "fdr")$table)
c=cpm(y0)[r, order(y0$samples$group)]
y=(tTags@.Data[[1]][1])
y <- data.frame(lapply(y, as.character), stringsAsFactors=FALSE) 
x=cbind(as.data.frame(tTags), uniq[y$genes,3],uniq[y$genes,4],c)
colnames(x)[7]="gene_name"
colnames(x)[8]="biotype"
#x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(qlf, p=0.05))
#write.table(tTags, file="edgeRgenes_EvsG_UntrPMNs_vs_PMNspolyPTIC_paired.csv",quote=FALSE)
write.table(x, file="edgeRgenes_EvsG_UntrPMNs_vs_PMNspolyPTIC_paired_featurecountExpr.csv",quote=FALSE,row.names=F)

#EvsH_UntrPMNs_vs_PMNspolyPIL29
SFI=rnaseqMatrix1[,c(23,34,16,32,13,25)]
colnames(SFI)
y0 <- DGEList(counts=SFI, genes=rnaseqMatrix1[,1])
y0 <- estimateDisp(y0,design)
fit <- glmQLFit(y0, design)
qlf <- glmQLFTest(fit)
topTags(qlf)
tTags = topTags(qlf,n=NULL,adjust.method = "fdr");r <- rownames(topTags(qlf,n=NULL,adjust.method = "fdr")$table)
c=cpm(y0)[r, order(y0$samples$group)]
y=(tTags@.Data[[1]][1])
y <- data.frame(lapply(y, as.character), stringsAsFactors=FALSE) 
x=cbind(as.data.frame(tTags), uniq[y$genes,3],uniq[y$genes,4],c)
colnames(x)[7]="gene_name"
colnames(x)[8]="biotype"
#x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(qlf, p=0.05))
#write.table(tTags, file="edgeRgenes_EvsH_UntrPMNs_vs_PMNspolyPIL29_paired.csv",quote=FALSE)
write.table(x, file="edgeRgenes_EvsH_UntrPMNs_vs_PMNspolyPIL29_paired_featurecountExpr.csv",quote=FALSE,row.names=F)

#EvsI_UntrPMNs_vs_PMNsTIC
SFI=rnaseqMatrix1[,c(23,8,16,35,13,29)]
colnames(SFI)
y0 <- DGEList(counts=SFI, genes=rnaseqMatrix1[,1])
y0 <- estimateDisp(y0,design)
fit <- glmQLFit(y0, design)
qlf <- glmQLFTest(fit)
topTags(qlf)
tTags = topTags(qlf,n=NULL,adjust.method = "fdr");r <- rownames(topTags(qlf,n=NULL,adjust.method = "fdr")$table)
c=cpm(y0)[r, order(y0$samples$group)]
y=(tTags@.Data[[1]][1])
y <- data.frame(lapply(y, as.character), stringsAsFactors=FALSE) 
x=cbind(as.data.frame(tTags), uniq[y$genes,3],uniq[y$genes,4],c)
colnames(x)[7]="gene_name"
colnames(x)[8]="biotype"
#x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(qlf, p=0.05))
#write.table(tTags, file="edgeRgenes_EvsI_UntrPMNs_vs_PMNsTIC_paired.csv",quote=FALSE)
write.table(x, file="edgeRgenes_EvsI_UntrPMNs_vs_PMNsTIC_paired_featurecountExpr.csv",quote=FALSE,row.names=F)

#EvsJ_UntrPMNs_vs_PMNsIL29
SFI=rnaseqMatrix1[,c(23,12,16,9,13,33)]
colnames(SFI)
y0 <- DGEList(counts=SFI, genes=rnaseqMatrix1[,1])
y0 <- estimateDisp(y0,design)
fit <- glmQLFit(y0, design)
qlf <- glmQLFTest(fit)
topTags(qlf)
tTags = topTags(qlf,n=NULL,adjust.method = "fdr");r <- rownames(topTags(qlf,n=NULL,adjust.method = "fdr")$table)
c=cpm(y0)[r, order(y0$samples$group)]
y=(tTags@.Data[[1]][1])
y <- data.frame(lapply(y, as.character), stringsAsFactors=FALSE) 
x=cbind(as.data.frame(tTags), uniq[y$genes,3],uniq[y$genes,4],c)
colnames(x)[7]="gene_name"
colnames(x)[8]="biotype"
#x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(qlf, p=0.05))
#write.table(tTags, file="edgeRgenes_EvsJ_UntrPMNs_vs_PMNsIL29_paired.csv",quote=FALSE)
write.table(x, file="edgeRgenes_EvsJ_UntrPMNs_vs_PMNsIL29_paired_featurecountExpr.csv",quote=FALSE,row.names=F)


#FvsG_PMNspolyP_vs_PMNspolyPTIC
SFI=rnaseqMatrix1[,c(27,31,20,24,17,21)]
colnames(SFI)
y0 <- DGEList(counts=SFI, genes=rnaseqMatrix1[,1])
y0 <- estimateDisp(y0,design)
fit <- glmQLFit(y0, design)
qlf <- glmQLFTest(fit)
topTags(qlf)
tTags = topTags(qlf,n=NULL,adjust.method = "fdr");r <- rownames(topTags(qlf,n=NULL,adjust.method = "fdr")$table)
c=cpm(y0)[r, order(y0$samples$group)]
y=(tTags@.Data[[1]][1])
y <- data.frame(lapply(y, as.character), stringsAsFactors=FALSE) 
x=cbind(as.data.frame(tTags), uniq[y$genes,3],uniq[y$genes,4],c)
colnames(x)[7]="gene_name"
colnames(x)[8]="biotype"
#x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(qlf, p=0.05))
#write.table(tTags, file="edgeRgenes_FvsG_PMNspolyP_vs_PMNspolyPTIC_paired.csv",quote=FALSE)
write.table(x, file="edgeRgenes_FvsG_PMNspolyP_vs_PMNspolyPTIC_paired_featurecountExpr.csv",quote=FALSE,row.names=F)


#FvsH_PMNspolyP_vs_PMNspolyPIL29
SFI=rnaseqMatrix1[,c(27,34,20,32,17,25)]
colnames(SFI)
y0 <- DGEList(counts=SFI, genes=rnaseqMatrix1[,1])
y0 <- estimateDisp(y0,design)
fit <- glmQLFit(y0, design)
qlf <- glmQLFTest(fit)
topTags(qlf)
tTags = topTags(qlf,n=NULL,adjust.method = "fdr");r <- rownames(topTags(qlf,n=NULL,adjust.method = "fdr")$table)
c=cpm(y0)[r, order(y0$samples$group)]
y=(tTags@.Data[[1]][1])
y <- data.frame(lapply(y, as.character), stringsAsFactors=FALSE) 
x=cbind(as.data.frame(tTags), uniq[y$genes,3],uniq[y$genes,4],c)
colnames(x)[7]="gene_name"
colnames(x)[8]="biotype"
#x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(qlf, p=0.05))
#write.table(tTags, file="edgeRgenes_FvsH_PMNspolyP_vs_PMNspolyPIL29_paired.csv",quote=FALSE)
write.table(x, file="edgeRgenes_FvsH_PMNspolyP_vs_PMNspolyPIL29_paired_featurecountExpr.csv",quote=FALSE,row.names=F)

#FvsI_PMNspolyP_vs_PMNsTIC
SFI=rnaseqMatrix1[,c(27,8,20,35,17,29)]
colnames(SFI)
y0 <- DGEList(counts=SFI, genes=rnaseqMatrix1[,1])
y0 <- estimateDisp(y0,design)
fit <- glmQLFit(y0, design)
qlf <- glmQLFTest(fit)
topTags(qlf)
tTags = topTags(qlf,n=NULL,adjust.method = "fdr");r <- rownames(topTags(qlf,n=NULL,adjust.method = "fdr")$table)
c=cpm(y0)[r, order(y0$samples$group)]
y=(tTags@.Data[[1]][1])
y <- data.frame(lapply(y, as.character), stringsAsFactors=FALSE) 
x=cbind(as.data.frame(tTags), uniq[y$genes,3],uniq[y$genes,4],c)
colnames(x)[7]="gene_name"
colnames(x)[8]="biotype"
#x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(qlf, p=0.05))
#write.table(tTags, file="edgeRgenes_FvsI_PMNspolyP_vs_PMNsTIC_paired.csv",quote=FALSE)
write.table(x, file="edgeRgenes_FvsI_PMNspolyP_vs_PMNsTIC_paired_featurecountExpr.csv",quote=FALSE,row.names=F)

#FvsJ_PMNspolyP_vs_PMNsIL29
SFI=rnaseqMatrix1[,c(27,12,20,9,17,33)]
colnames(SFI)
y0 <- DGEList(counts=SFI, genes=rnaseqMatrix1[,1])
y0 <- estimateDisp(y0,design)
fit <- glmQLFit(y0, design)
qlf <- glmQLFTest(fit)
topTags(qlf)
tTags = topTags(qlf,n=NULL,adjust.method = "fdr");r <- rownames(topTags(qlf,n=NULL,adjust.method = "fdr")$table)
c=cpm(y0)[r, order(y0$samples$group)]
y=(tTags@.Data[[1]][1])
y <- data.frame(lapply(y, as.character), stringsAsFactors=FALSE) 
x=cbind(as.data.frame(tTags), uniq[y$genes,3],uniq[y$genes,4],c)
colnames(x)[7]="gene_name"
colnames(x)[8]="biotype"
#x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(qlf, p=0.05))
#write.table(tTags, file="edgeRgenes_FvsJ_PMNspolyP_vs_PMNsIL29_paired.csv",quote=FALSE)
write.table(x, file="edgeRgenes_FvsJ_PMNspolyP_vs_PMNsIL29_paired_featurecountExpr.csv",quote=FALSE,row.names=F)


#GvsI_PMNspolyPTIC_vs_PMNsTIC
SFI=rnaseqMatrix1[,c(31,8,24,35,21,29)]
colnames(SFI)
y0 <- DGEList(counts=SFI, genes=rnaseqMatrix1[,1])
y0 <- estimateDisp(y0,design)
fit <- glmQLFit(y0, design)
qlf <- glmQLFTest(fit)
topTags(qlf)
tTags = topTags(qlf,n=NULL,adjust.method = "fdr");r <- rownames(topTags(qlf,n=NULL,adjust.method = "fdr")$table)
c=cpm(y0)[r, order(y0$samples$group)]
y=(tTags@.Data[[1]][1])
y <- data.frame(lapply(y, as.character), stringsAsFactors=FALSE) 
x=cbind(as.data.frame(tTags), uniq[y$genes,3],uniq[y$genes,4],c)
colnames(x)[7]="gene_name"
colnames(x)[8]="biotype"
#x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(qlf, p=0.05))
#write.table(tTags, file="edgeRgenes_GvsI_PMNspolyPTIC_vs_PMNsTIC_paired.csv",quote=FALSE)
write.table(x, file="edgeRgenes_GvsI_PMNspolyPTIC_vs_PMNsTIC_paired_featurecountExpr.csv",quote=FALSE,row.names=F)

#GvsJ_PMNspolyPTIC_vs_PMNsIL29
SFI=rnaseqMatrix1[,c(31,12,24,9,21,33)]
colnames(SFI)
y0 <- DGEList(counts=SFI, genes=rnaseqMatrix1[,1])
y0 <- estimateDisp(y0,design)
fit <- glmQLFit(y0, design)
qlf <- glmQLFTest(fit)
topTags(qlf)
tTags = topTags(qlf,n=NULL,adjust.method = "fdr");r <- rownames(topTags(qlf,n=NULL,adjust.method = "fdr")$table)
c=cpm(y0)[r, order(y0$samples$group)]
y=(tTags@.Data[[1]][1])
y <- data.frame(lapply(y, as.character), stringsAsFactors=FALSE) 
x=cbind(as.data.frame(tTags), uniq[y$genes,3],uniq[y$genes,4],c)
colnames(x)[7]="gene_name"
colnames(x)[8]="biotype"
#x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(qlf, p=0.05))
#write.table(tTags, file="edgeRgenes_GvsJ_PMNspolyPTIC_vs_PMNsIL29_paired.csv",quote=FALSE)
write.table(x, file="edgeRgenes_GvsJ_PMNspolyPTIC_vs_PMNsIL29_paired_featurecountExpr.csv",quote=FALSE,row.names=F)