/data/results/tools/rnaseq/subread/subread-1.5.2-source/bin/featureCounts -O -g gene_name -t exon -T 8 -a ~/bak/doc/fleming/kafasla/DKlab/mr/Polysomes/reads/stringtie/merged.gtf -o gene_counts_merged3_featurecounts.txt ../../tophat/KO_0h_replicate_I/accepted_hits.bam ../../tophat/KO_0h_replicate_II/accepted_hits.bam ../../tophat/KO_2h_replicate_I/accepted_hits.bam ../../tophat/KO_2h_replicate_II/accepted_hits.bam ../../tophat/KO_6h_replicate_I/accepted_hits.bam ../../tophat/KO_6h_replicate_II/accepted_hits.bam ../../tophat/KO_IFN-gamma_replicate_II/accepted_hits.bam ../../tophat/KO_IFN_replicate_I/accepted_hits.bam ../../tophat/KO_IL4_replicate_I/accepted_hits.bam ../../tophat/KO_IL4_replicate_II/accepted_hits.bam ../../tophat/WT_0h_replicate_I/accepted_hits.bam ../../tophat/WT_0h_replicate_II/accepted_hits.bam ../../tophat/WT_2h_replicate_I/accepted_hits.bam ../../tophat/WT_2h_replicate_II/accepted_hits.bam ../../tophat/WT_6h_replicate_I/accepted_hits.bam ../../tophat/WT_6h_replicate_II/accepted_hits.bam ../../tophat/WT_IFN-gamma_replicate_I/accepted_hits.bam ../../tophat/WT_IFN-gamma_replicate_II/accepted_hits.bam ../../tophat/WT_IL4_replicate_I/accepted_hits.bam ../../tophat/WT_IL4_replicate_II/accepted_hits.bam 

#                1                  2                 3                  4                 5                  6                         7                  8                  9                  10                11                 12                13                 14                15                 16                       17                        18                 19                  20 
#,KO_0h_replicate_I,KO_0h_replicate_II,KO_2h_replicate_I,KO_2h_replicate_II,KO_6h_replicate_I,KO_6h_replicate_II,KO_IFN-gamma_replicate_II,KO_IFN_replicate_I,KO_IL4_replicate_I,KO_IL4_replicate_II,WT_0h_replicate_I,WT_0h_replicate_II,WT_2h_replicate_I,WT_2h_replicate_II,WT_6h_replicate_I,WT_6h_replicate_II,WT_IFN-gamma_replicate_I,WT_IFN-gamma_replicate_II,WT_IL4_replicate_I,WT_IL4_replicate_II

require(edgeR)
# FC sign:  correct is ko vs wt 
rnaseqMatrix1= read.table("gene_counts_merged3_featurecounts.txt",header=T,sep="\t",skip=1)
rnaseqMatrix0 =rnaseqMatrix1[,7:dim(rnaseqMatrix1)[2]];
rownames(rnaseqMatrix0)=rnaseqMatrix1$Geneid;

conditions = factor(c("KO","KO","KO","KO","KO","KO","KO","KO","KO","KO","wt","wt","wt","wt","wt","wt","wt","wt","wt","wt"))
rnaseqMatrix =rnaseqMatrix0;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,pair=c("WT","KO"))
tTags = topTags(et,n=NULL,adjust.method = "fdr")
summary(de <- decideTestsDGE(et, p=0.05))
write.table(tTags, file="edgeRgenes_featurecount-allKO-vs-allWT.csv")
#-1   539
#0  28650
#1    589


col_ordering = c(1,2,11,12)#0h
rnaseqMatrix =rnaseqMatrix0[,col_ordering];
#rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,];
rnaseqMatrix = rnaseqMatrix[(rowSums(rnaseqMatrix[,1:2])>=5)|(rowSums(rnaseqMatrix[,3:4])>=5),];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,pair=c("WT","KO"))
tTags = topTags(et,n=NULL,adjust.method = "fdr");r <- rownames(topTags(et,n=NULL,adjust.method = "fdr")$table)
c=cpm(exp_study)[r, order(exp_study$samples$group)]
x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(et, p=0.05))
write.table(tTags, file="edgeRgenes_featurecount-0h.csv")
write.table(x, file="edgeRgenes_featurecountExpr-0h.csv")
#-1   487
#0  16655
#1    177

col_ordering = c(3,4,13,14)#2h
rnaseqMatrix =rnaseqMatrix0[,col_ordering];
#rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,];
rnaseqMatrix = rnaseqMatrix[(rowSums(rnaseqMatrix[,1:2])>=5)|(rowSums(rnaseqMatrix[,3:4])>=5),];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,pair=c("WT","KO"))
tTags = topTags(et,n=NULL,adjust.method = "fdr");r <- rownames(topTags(et,n=NULL,adjust.method = "fdr")$table)
c=cpm(exp_study)[r, order(exp_study$samples$group)]
x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(et, p=0.05))
write.table(tTags, file="edgeRgenes_featurecount-2h.csv")
write.table(x, file="edgeRgenesExpr_featurecountExpr-2h.csv")
#-1    79
#0  24820
#1     49

col_ordering = c(5,6,15,16)#6h
rnaseqMatrix =rnaseqMatrix0[,col_ordering];
#rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,];conditions = factor(c("KO","KO","WT","WT"));
rnaseqMatrix = rnaseqMatrix[(rowSums(rnaseqMatrix[,1:2])>=5)|(rowSums(rnaseqMatrix[,3:4])>=5),];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,pair=c("WT","KO"))
tTags = topTags(et,n=NULL,adjust.method = "fdr");r <- rownames(topTags(et,n=NULL,adjust.method = "fdr")$table)
c=cpm(exp_study)[r, order(exp_study$samples$group)]
x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(et, p=0.05))
write.table(tTags, file="edgeRgenes_featurecount-6h.csv")
write.table(x, file="edgeRgenes_featurecountExpr-6h.csv")
#-1   113
#0  24993
#1     49


col_ordering = c(7,8,17,18)#IFN
rnaseqMatrix =rnaseqMatrix0[,col_ordering];
#rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,];
rnaseqMatrix = rnaseqMatrix[(rowSums(rnaseqMatrix[,1:2])>=5)|(rowSums(rnaseqMatrix[,3:4])>=5),];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,pair=c("WT","KO"))
tTags = topTags(et,n=NULL,adjust.method = "fdr");r <- rownames(topTags(et,n=NULL,adjust.method = "fdr")$table)
c=cpm(exp_study)[r, order(exp_study$samples$group)]
x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(et, p=0.05))
write.table(tTags, file="edgeRgenes_featurecount-IFN.csv")
write.table(x, file="edgeRgenesExpr_featurecount-IFN.csv")
#-1    30
#0  25364
#1     33

col_ordering = c(9,10,19,20)#IL4
rnaseqMatrix =rnaseqMatrix0[,col_ordering];
#rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,];conditions = factor(c("KO","KO","WT","WT"));
rnaseqMatrix = rnaseqMatrix[(rowSums(rnaseqMatrix[,1:2])>=5)|(rowSums(rnaseqMatrix[,3:4])>=5),];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,pair=c("WT","KO"))
tTags = topTags(et,n=NULL,adjust.method = "fdr");r <- rownames(topTags(et,n=NULL,adjust.method = "fdr")$table)
c=cpm(exp_study)[r, order(exp_study$samples$group)]
x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(et, p=0.05))
write.table(tTags, file="edgeRgenes_featurecount-IL4.csv")
write.table(x, file="edgeRgenes_featurecountExpr-IL4.csv")
#-1   185
#0  23732
#1    119





#@
/data/images/proton/DKlab/mr/Polysomes/reads/stringtie/per-gene/edgeRgenes_featurecountExpr-IL4.csv
"logFC" "logCPM" "PValue" "FDR" "......tophat.KO_IL4_replicate_I.accepted_hits.bam" "......tophat.KO_IL4_replicate_II.accepted_hits.bam" "......tophat.WT_IL4_replicate_I.accepted_hits.bam" "......tophat.WT_IL4_replicate_II.accepted_hits.bam"
"Xist" 0.956230223708131 7.78347077548186 0.166645369656576 0.608522705076508 377.771245244346 191.592721041537 255.305205596653 35.6729749819567

ENSMUST00000153883
ENSMUST00000127786
"logFC" "logCPM" "PValue" "FDR" "KO_IL4_replicate_I" "KO_IL4_replicate_II" "WT_IL4_replicate_I" "WT_IL4_replicate_II"
"ENSMUST00000127786" 0.392687693808114 6.99403139812579 0.755070877485269 1 128.020802025999 84.6737636991242 254.773312979144 23.7041260774964
"ENSMUST00000153883" -5.01334076972655 3.07011562647892 0.0409379450781329 0.601926786734293 11.6382547296363 4.70409798328467 0 0