/data/results/tools/rnaseq/subread/subread-1.5.2-source/bin/featureCounts" "-O" "-g" "gene_name" "-t" "exon" "-T" "8" "-a" "/data/images/proton/DKlab/orsalia2/stringtie/merged.gtf3" "-o" "gene_counts_merged3.txt" "/data/images/proton/DKlab/orsalia/sortedwt01pos.sorted.bam" "/data/images/proton/DKlab/orsalia/sortedwt02pos.sorted.bam" "/data/images/proton/DKlab/orsalia/sortedwt21pos.bam" "/data/images/proton/DKlab/orsalia/sortedwt22pos.bam" "/data/images/proton/DKlab/orsalia/sortedwt61pos.bam" "/data/images/proton/DKlab/orsalia/sortedwt62pos.bam" "/data/images/proton/DKlab/orsalia/sortedwtIFN1pos.bam" "/data/images/proton/DKlab/orsalia/sortedwtIFN2pos.bam" "/data/images/proton/DKlab/orsalia/sortedwt01IL4pos.bam" "/data/images/proton/DKlab/orsalia/sortedwt02IL4pos.bam" "/data/images/proton/DKlab/orsalia/sortedKo01pos.sorted.bam" "/data/images/proton/DKlab/orsalia/sortedKo02pos.sorted.bam" "/data/images/proton/DKlab/orsalia/sortedKO21pos.bam" "/data/images/proton/DKlab/orsalia/sortedKO22pos.bam" "/data/images/proton/DKlab/orsalia/sortedKo61pos.bam" "/data/images/proton/DKlab/orsalia/sortedKo62pos.bam" "/data/images/proton/DKlab/orsalia/sortedKOIFN1pos.bam" "/data/images/proton/DKlab/orsalia/sortedKOIFN2pos.bam" "/data/images/proton/DKlab/orsalia/sortedKo1IL4pos.bam" "/data/images/proton/DKlab/orsalia/sortedKo2IL4pos.bam 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 sortedwt01pos.sorted.bam" sortedwt02pos.sorted.bam" sortedwt21pos.bam" sortedwt22pos.bam" sortedwt61pos.bam" sortedwt62pos.bam" sortedwtIFN1pos.bam" sortedwtIFN2pos.bam" sortedwt01IL4pos.bam" sortedwt02IL4pos.bam" sortedKo01pos.sorted.bam" sortedKo02pos.sorted.bam" sortedKO21pos.bam" sortedKO22pos.bam" sortedKo61pos.bam" sortedKo62pos.bam" sortedKOIFN1pos.bam" sortedKOIFN2pos.bam" sortedKo1IL4pos.bam" sortedKo2IL4pos.bam require(edgeR) rnaseqMatrix1= read.table("gene_counts_merged4_featurecounts.txt",header=T,sep="\t",skip=1) rnaseqMatrix0 =rnaseqMatrix1[,7:dim(rnaseqMatrix1)[2]]; rownames(rnaseqMatrix0)=rnaseqMatrix1$Geneid; #rnaseqMatrix = round(rnaseqMatrix) rnaseqMatrix = rnaseqMatrix0[rowSums(rnaseqMatrix0)>=2,] #conditions = factor(c(rep("FAGsA", 2), rep("FAGsB", 2))) # 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 #conditions = factor(c("wt0","wt0","wt2","wt2","wt6","wt6","wtIFN","wtIFN","wt0IL4","wt0IL4","Ko0","Ko0","KO2","KO2","Ko6","Ko6","KOIFN","KOIFN","KoIL4","KoIL4")) conditions = factor(c("wt","wt","wt","wt","wt","wt","wt","wt","wt","wt","KO","KO","KO","KO","KO","KO","KO","KO","KO","KO")) 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)$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)) #detags <- rownames(exp_study)[as.logical(de)] #plotSmear(et, de.tags=detags) #abline(h = c(-1, 1), col = "blue") write.table(tTags, file="edgeRgenes_featurecount.csv") write.table(x, file="edgeRgenes_featurecountExpr.csv") col_ordering = c(1,2,11,12)#0h rnaseqMatrix =rnaseqMatrix0[,col_ordering];rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,];conditions = factor(c("WT","WT","KO","KO")); 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)) #detags <- rownames(exp_study)[as.logical(de)] #plotSmear(et, de.tags=detags) #abline(h = c(-1, 1), col = "blue") write.table(tTags, file="edgeRgenes_featurecount-0h.csv") write.table(x, file="edgeRgenes_featurecountExpr-0h.csv") col_ordering = c(3,4,13,14)#2h rnaseqMatrix =rnaseqMatrix0[,col_ordering];rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,];conditions = factor(c("WT","WT","KO","KO")); 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")) 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) write.table(tTags, file="edgeRgenes_featurecount-2h.csv") write.table(x, file="edgeRgenes_featurecountExpr-2h.csv") col_ordering = c(5,6,15,16)#6h rnaseqMatrix =rnaseqMatrix0[,col_ordering];rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,];conditions = factor(c("WT","WT","KO","KO")); 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) write.table(tTags, file="edgeRgenes_featurecount-6h.csv") write.table(x, file="edgeRgenes_featurecountExpr-6h.csv") col_ordering = c(7,8,17,18)#IFN rnaseqMatrix =rnaseqMatrix0[,col_ordering];rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,];conditions = factor(c("WT","WT","KO","KO")); 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) write.table(tTags, file="edgeRgenes_featurecount-IFN.csv") write.table(x, file="edgeRgenes_featurecountExpr-IFN.csv") col_ordering = c(9,10,19,20)#IL4 rnaseqMatrix =rnaseqMatrix0[,col_ordering];rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,];conditions = factor(c("WT","WT","KO","KO")); 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) write.table(tTags, file="edgeRgenes_featurecount-IL4.csv") write.table(x, file="edgeRgenes_featurecountExpr-IL4.csv") awk -f mergeEdgerDEresults1.awk all_pairs_edger.txt > edgeRgenes_featurecount_all_pairs.csv