#/data/results/tools/rnaseq/subread/subread-1.5.2-source/bin/featureCounts -O -g "gene_id" -t exon -T 8  -a /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.gtf -o gene_counts_merged_NCBIM37.64_id.txt ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/KO_0h_replicate_I/accepted_hits.bam ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/KO_0h_replicate_II/accepted_hits.bam ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/KO_2h_replicate_I/accepted_hits.bam ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/KO_2h_replicate_II/accepted_hits.bam ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/KO_6h_replicate_I/accepted_hits.bam ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/KO_6h_replicate_II/accepted_hits.bam ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/KO_IFN-gamma_replicate_II/accepted_hits.bam ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/KO_IFN_replicate_I/accepted_hits.bam ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/KO_IL4_replicate_I/accepted_hits.bam ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/KO_IL4_replicate_II/accepted_hits.bam ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/WT_0h_replicate_I/accepted_hits.bam ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/WT_0h_replicate_II/accepted_hits.bam ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/WT_2h_replicate_I/accepted_hits.bam ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/WT_2h_replicate_II/accepted_hits.bam ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/WT_6h_replicate_I/accepted_hits.bam ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/WT_6h_replicate_II/accepted_hits.bam ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/WT_IFN-gamma_replicate_I/accepted_hits.bam ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/WT_IFN-gamma_replicate_II/accepted_hits.bam ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/WT_IL4_replicate_I/accepted_hits.bam ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/WT_IL4_replicate_II/accepted_hits.bam

/data/results/tools/rnaseq/subread/subread-1.5.2-source/bin/featureCounts -O -g "gene_id" -t exon -T 8  -a /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.gtf -o gene_counts_merged_NCBIM37.64_id.txt KO_0h_replicate_I.bam KO_0h_replicate_II.bam KO_2h_replicate_I.bam KO_2h_replicate_II.bam KO_6h_replicate_I.bam KO_6h_replicate_II.bam KO_IFN-gamma_replicate_II.bam KO_IFN_replicate_I.bam KO_IL4_replicate_I.bam KO_IL4_replicate_II.bam WT_0h_replicate_I.bam WT_0h_replicate_II.bam WT_2h_replicate_I.bam WT_2h_replicate_II.bam WT_6h_replicate_I.bam WT_6h_replicate_II.bam WT_IFN-gamma_replicate_I.bam WT_IFN-gamma_replicate_II.bam WT_IL4_replicate_I.bam WT_IL4_replicate_II.bam 

ln -s ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/KO_0h_replicate_I/accepted_hits.bam KO_0h_replicate_I.bam
ln -s ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/KO_0h_replicate_II/accepted_hits.bam KO_0h_replicate_II.bam
ln -s ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/KO_2h_replicate_I/accepted_hits.bam KO_2h_replicate_I.bam
ln -s ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/KO_2h_replicate_II/accepted_hits.bam KO_2h_replicate_II.bam
ln -s ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/KO_6h_replicate_I/accepted_hits.bam KO_6h_replicate_I.bam
ln -s ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/KO_6h_replicate_II/accepted_hits.bam KO_6h_replicate_II.bam
ln -s ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/KO_IFN-gamma_replicate_II/accepted_hits.bam KO_IFN-gamma_replicate_II.bam
ln -s ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/KO_IFN_replicate_I/accepted_hits.bam KO_IFN_replicate_I.bam
ln -s ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/KO_IL4_replicate_I/accepted_hits.bam KO_IL4_replicate_I.bam
ln -s ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/KO_IL4_replicate_II/accepted_hits.bam KO_IL4_replicate_II.bam
ln -s ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/WT_0h_replicate_I/accepted_hits.bam WT_0h_replicate_I.bam
ln -s ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/WT_0h_replicate_II/accepted_hits.bam WT_0h_replicate_II.bam
ln -s ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/WT_2h_replicate_I/accepted_hits.bam WT_2h_replicate_I.bam
ln -s ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/WT_2h_replicate_II/accepted_hits.bam WT_2h_replicate_II.bam
ln -s ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/WT_6h_replicate_I/accepted_hits.bam WT_6h_replicate_I.bam
ln -s ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/WT_6h_replicate_II/accepted_hits.bam WT_6h_replicate_II.bam
ln -s ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/WT_IFN-gamma_replicate_I/accepted_hits.bam WT_IFN-gamma_replicate_I.bam
ln -s ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/WT_IFN-gamma_replicate_II/accepted_hits.bam WT_IFN-gamma_replicate_II.bam
ln -s ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/WT_IL4_replicate_I/accepted_hits.bam WT_IL4_replicate_I.bam
ln -s ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/tophat/WT_IL4_replicate_II/accepted_hits.bam WT_IL4_replicate_II.bam


# Rscript for edgeR in /data/images/proton/DKlab/mr/Polysomes/reads/stringtie/per-gene/README-polysomeSeq-per-gene.txt

require(edgeR)

gl=read.table("~/bak/doc/fleming/kafasla/DKlab/mr/rnaseq/deg/Mus_musculus.NCBIM37.64-toMM9.gene-lenghts.txt",sep="\t")
gl.names=gl$V1;

#rnaseqMatrix1= read.table("gene_counts_merged_k20_NCBIM37.64.txt",header=T,sep="\t",skip=1)
rnaseqMatrix1= read.table("~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/deg/gene_counts_merged_NCBIM37.64_id.txt",header=T,sep="\t",skip=1)
rnaseqMatrix0 =rnaseqMatrix1[,7:dim(rnaseqMatrix1)[2]];
rownames(rnaseqMatrix0)=rnaseqMatrix1$Geneid;
#rnaseqMatrix = rnaseqMatrix0[rowSums(rnaseqMatrix0)>=2,]
 for (i in 1:dim(rnaseqMatrix0)[2]) {
print (colnames(rnaseqMatrix0)[i])
print (summary(rnaseqMatrix0[,i]>0))
}

 colnames(rnaseqMatrix0)
 [1] "KO_0h_replicate_I.bam"         "KO_0h_replicate_II.bam"       
 [3] "KO_2h_replicate_I.bam"         "KO_2h_replicate_II.bam"       
 [5] "KO_6h_replicate_I.bam"         "KO_6h_replicate_II.bam"       
 [7] "KO_IFN.gamma_replicate_II.bam" "KO_IFN_replicate_I.bam"       
 [9] "KO_IL4_replicate_I.bam"        "KO_IL4_replicate_II.bam"      
[11] "WT_0h_replicate_I.bam"         "WT_0h_replicate_II.bam"       
[13] "WT_2h_replicate_I.bam"         "WT_2h_replicate_II.bam"       
[15] "WT_6h_replicate_I.bam"         "WT_6h_replicate_II.bam"       
[17] "WT_IFN.gamma_replicate_I.bam"  "WT_IFN.gamma_replicate_II.bam"
[19] "WT_IL4_replicate_I.bam"        "WT_IL4_replicate_II.bam"

col_ordering = c(1,2,11,12)#0h
rnaseqMatrix =rnaseqMatrix0[,col_ordering];
rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,]
conditions = factor(c("KO","KO","WT","WT"));
exp_study = DGEList(counts=rnaseqMatrix, group=conditions);exp_study = calcNormFactors(exp_study);exp_study = estimateCommonDisp(exp_study);exp_study = estimateTagwiseDisp(exp_study);et = exactTest(exp_study)
tTags = topTags(et,n=NULL,adjust.method = "fdr")
#r <- rownames(topTags(et)$table)
r <- rownames(tTags)
c=cpm(exp_study)[r, order(exp_study$samples$group)]
#x=cbind(as.data.frame(tTags),r)
x=cbind(as.data.frame(tTags), gl[r,2],gl[r,3])
x=cbind(as.data.frame(tTags), "-",1)
x[,5]=lapply(x[,5],as.character)
for (i in 1:length(r)){
    y=gl[(gl.names==r[i]),];
    x[i,5]=as.character(y$V2);x[i,6]=y$V3;
    }
x1=cbind(x,c,c[,1]*1000/x[,6],c[,2]*1000/x[,6],c[,3]*1000/x[,6],c[,4]*1000/x[,6])
x2=cbind(x1,0.5*(x1[,11]+x1[,12]),0.5*(x1[,13]+x1[,14]))
y=x2[ (x2[,15]>1)&(x2[,16]>1),]
plot(log2(y[,15]),log2(y[,16]),pch=".")
#x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(et, p=0.05))
#detags <- rownames(exp_study)[as.logical(de)]
#plotSmear(et, de.tags=detags)
#abline(h = c(-1, 1), col = "blue")
write.table(tTags, file="edgeRgenes-0h.csv")
write.table(x, file="edgeRgenesExpr-0h.csv")
#-1   101
#0  22027
#1    113

col_ordering = c(17,18,7,8)#IFN
rnaseqMatrix =rnaseqMatrix0[,col_ordering];
rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,]
exp_study = DGEList(counts=rnaseqMatrix, group=conditions);exp_study = calcNormFactors(exp_study);exp_study = estimateCommonDisp(exp_study);exp_study = estimateTagwiseDisp(exp_study);et = exactTest(exp_study)
tTags = topTags(et,n=NULL,adjust.method = "fdr")
r <- rownames(tTags)
c=cpm(exp_study)[r, order(exp_study$samples$group)]
x=cbind(as.data.frame(tTags), gl[rownames(tTags),2],gl[rownames(tTags),3],c,0.001*gl[rownames(tTags),3]*c)
x=cbind(x,0.5*(x[,10]+x[,11]),0.5*(x[,12]+x[,13]))
y=x[ (x[,14]>1)&(x[,15]>1),]
write.table(tTags, file="edgeRgenes-IFN.csv")
write.table(x, file="edgeRgenesExpr-IFN.csv")
summary(de <- decideTestsDGE(et, p=0.05))
#-1    73
#0  19924
#1    126

col_ordering = c(19,20,9,10)#IL4
rnaseqMatrix =rnaseqMatrix0[,col_ordering];
rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,]
exp_study = DGEList(counts=rnaseqMatrix, group=conditions);exp_study = calcNormFactors(exp_study);exp_study = estimateCommonDisp(exp_study);exp_study = estimateTagwiseDisp(exp_study);et = exactTest(exp_study)
tTags = topTags(et,n=NULL,adjust.method = "fdr")
r <- rownames(tTags)
c=cpm(exp_study)[r, order(exp_study$samples$group)]
x=cbind(as.data.frame(tTags), gl[rownames(tTags),2],gl[rownames(tTags),3],c,0.001*gl[rownames(tTags),3]*c)
x=cbind(x,0.5*(x[,10]+x[,11]),0.5*(x[,12]+x[,13]))
y=x[ (x[,14]>1)&(x[,15]>1),]
write.table(tTags, file="edgeRgenes-IL4.csv")
write.table(x, file="edgeRgenesExpr-IL4.csv")
summary(de <- decideTestsDGE(et, p=0.05))

# get > 1 RPKM genes for final table

#genes detected
   Mode   FALSE    TRUE    NA's 
logical   14840   22470       0 
[1] "KO_0h_replicate_II.bam"
   Mode   FALSE    TRUE    NA's 
logical   14386   22924       0 
[1] "KO_2h_replicate_I.bam"
   Mode   FALSE    TRUE    NA's 
logical   13672   23638       0 
[1] "KO_2h_replicate_II.bam"
   Mode   FALSE    TRUE    NA's 
logical   14114   23196       0 
[1] "KO_6h_replicate_I.bam"
   Mode   FALSE    TRUE    NA's 
logical   12612   24698       0 
[1] "KO_6h_replicate_II.bam"
   Mode   FALSE    TRUE    NA's 
logical   13236   24074       0 
[1] "KO_IFN.gamma_replicate_II.bam"
   Mode   FALSE    TRUE    NA's 
logical   13455   23855       0 
[1] "KO_IFN_replicate_I.bam"
   Mode   FALSE    TRUE    NA's 
logical   11884   25426       0 
[1] "KO_IL4_replicate_I.bam"
   Mode   FALSE    TRUE    NA's 
logical   14282   23028       0 
[1] "KO_IL4_replicate_II.bam"
   Mode   FALSE    TRUE    NA's 
logical   15362   21948       0 
[1] "WT_0h_replicate_I.bam"
   Mode   FALSE    TRUE    NA's 
logical   14747   22563       0 
[1] "WT_0h_replicate_II.bam"
   Mode   FALSE    TRUE    NA's 
logical   17350   19960       0 
[1] "WT_2h_replicate_I.bam"
   Mode   FALSE    TRUE    NA's 
logical   14596   22714       0 
[1] "WT_2h_replicate_II.bam"
   Mode   FALSE    TRUE    NA's 
logical   15803   21507       0 
[1] "WT_6h_replicate_I.bam"
   Mode   FALSE    TRUE    NA's 
logical   17374   19936       0 
[1] "WT_6h_replicate_II.bam"
   Mode   FALSE    TRUE    NA's 
logical   14727   22583       0 
[1] "WT_IFN.gamma_replicate_I.bam"
   Mode   FALSE    TRUE    NA's 
logical   16354   20956       0 
[1] "WT_IFN.gamma_replicate_II.bam"
   Mode   FALSE    TRUE    NA's 
logical   14696   22614       0 
[1] "WT_IL4_replicate_I.bam"
   Mode   FALSE    TRUE    NA's 
logical   15254   22056       0 
[1] "WT_IL4_replicate_II.bam"
   Mode   FALSE    TRUE    NA's 
logical   17127   20183       0 
>

#genes >1read