reczko@fix:~/bak/doc/fleming/kafasla/DKlab/mr/rnaseq-vs-polysomeseq$ R R version 3.2.0 (2015-04-16) -- "Full of Ingredients" Copyright (C) 2015 The R Foundation for Statistical Computing Platform: x86_64-unknown-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > require(edgeR) gl=read.table("/mnt/max/b/genomics_facility/DKlab/mr/tzanos/mm10-GRCm38.p6/genome/mm10b2.gene-lenghts.txt",sep="\t") gl.names=gl$V1; rownames(gl)=gl$V1 colnames(gl)=c("ID","geneID","genelength"); #rnaseq data rnaseqMatrix1= read.table("/mnt/max/b/genomics_facility/DKlab/mr/tzanos/mm10-GRCm38.p6/rnaseq/featurecounts/gene_counts_merged_k20_mm10.txt",header=T,sep="\t",skip=1) rnaseqMatrix0 =rnaseqMatrix1[,7:dim(rnaseqMatrix1)[2]]; rownames(rnaseqMatrix0)=as.character(rnaseqMatrix1$Geneid); colnames(rnaseqMatrix0) # [1] "wt01_trm_k20.bam" "wt02_trm_k20.bam" "wt21_trm_k20.bam" # [4] "wt22_trm_k20.bam" "wt61_trm_k20.bam" "wt62_trm_k20.bam" # [7] "wtIFN_1_trm_k20.bam" "wtIFN_2_trm_k20.bam" "wtIL4_1_trm_k20.bam" #[10] "wtIL4_2_trm_k20.bam" "KO01_trm_k20.bam" "KO02_trm_k20.bam" #[13] "KO21_trm_k20.bam" "KO22_trm_k20.bam" "KO61_trm_k20.bam" #[16] "KO62_trm_k20.bam" "KOIFN_1_trm_k20.bam" "KOIFN_2_trm_k20.bam" #[19] "KOIL4_1_trm_k20.bam" "KOIL4_2_trm_k20.bam" col_ordering = c(9,10,7,8)#WT_M1vsM2 rnaseqMatrix =rnaseqMatrix0[,col_ordering]; rnaseqMatrix = rnaseqMatrix[(rowSums(rnaseqMatrix[,1:2])>=2)&(rowSums(rnaseqMatrix[,3:4])>=2),] conditions = factor(c("IFN","IFN","IL4","IL4")); 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)] rx=merge(x=as.data.frame(tTags),y=gl,by="row.names",all.x=TRUE) rx=rx[,!(names(rx) %in% "ID")] rownames(rx)=rx$Row.names; rx=rx[,!(names(rx) %in% "Row.names")] rx2=merge(x=rx,y=c,by="row.names",all.x=TRUE) rx=rx2[,!(names(rx2) %in% "ID")] rx=cbind(rx,1000*rx[,8]/rx$genelength,1000*rx[,9]/rx$genelength,1000*rx[,10]/rx$genelength,1000*rx[,11]/rx$genelength) rx=cbind(rx,0.5*(rx[,12]+rx[,13]),0.5*(rx[,14]+rx[,15])) write.table(rx, file="edgeRgenesExpr-WT_M1vsM2.csv",quote=FALSE,row.names=FALSE) rx_m1=rx for (i in 8:11){ print (colnames(rx)[i]) tt=summary(rx[,i]>0) print(tt) } Loading required package: edgeR Loading required package: limma > > > > > > > > > > > [1] "wt01_trm_k20.bam" "wt02_trm_k20.bam" "wt21_trm_k20.bam" [4] "wt22_trm_k20.bam" "wt61_trm_k20.bam" "wt62_trm_k20.bam" [7] "wtIFN_1_trm_k20.bam" "wtIFN_2_trm_k20.bam" "wtIL4_1_trm_k20.bam" [10] "wtIL4_2_trm_k20.bam" "KO01_trm_k20.bam" "KO02_trm_k20.bam" [13] "KO21_trm_k20.bam" "KO22_trm_k20.bam" "KO61_trm_k20.bam" [16] "KO62_trm_k20.bam" "KOIFN_1_trm_k20.bam" "KOIFN_2_trm_k20.bam" [19] "KOIL4_1_trm_k20.bam" "KOIL4_2_trm_k20.bam" > > > > > > > > > > > > > > > > > > > > > > > > > > > > + + + + [1] "wtIL4_1_trm_k20.bam" Mode FALSE TRUE NA's logical 469 20510 0 [1] "wtIL4_2_trm_k20.bam" Mode FALSE TRUE NA's logical 709 20270 0 [1] "wtIFN_1_trm_k20.bam" Mode FALSE TRUE NA's logical 1037 19942 0 [1] "wtIFN_2_trm_k20.bam" Mode FALSE TRUE NA's logical 668 20311 0 > > col_ordering = c(19,20,17,18)#KO_M1vsM2 rnaseqMatrix =rnaseqMatrix0[,col_ordering]; rnaseqMatrix = rnaseqMatrix[(rowSums(rnaseqMatrix[,1:2])>=2)&(rowSums(rnaseqMatrix[,3:4])>=2),] conditions = factor(c("IFN","IFN","IL4","IL4")); 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)] rx=merge(x=as.data.frame(tTags),y=gl,by="row.names",all.x=TRUE) rx=rx[,!(names(rx) %in% "ID")] rownames(rx)=rx$Row.names; rx=rx[,!(names(rx) %in% "Row.names")] rx2=merge(x=rx,y=c,by="row.names",all.x=TRUE) rx=rx2[,!(names(rx2) %in% "ID")] rx=cbind(rx,1000*rx[,8]/rx$genelength,1000*rx[,9]/rx$genelength,1000*rx[,10]/rx$genelength,1000*rx[,11]/rx$genelength) rx=cbind(rx,0.5*(rx[,12]+rx[,13]),0.5*(rx[,14]+rx[,15])) write.table(rx, file="edgeRgenesExpr-KO_M1vsM2.csv",quote=FALSE,row.names=FALSE) rx_m2=rx for (i in 8:11){ print (colnames(rx)[i]) tt=summary(rx[,i]>0) print(tt) } #[1] "KOIL4_1_trm_k20.bam" #Mode FALSE TRUE #logical 514 20577 #[1] "KOIL4_2_trm_k20.bam" #Mode FALSE TRUE #logical 518 20573 #[1] "KOIFN_1_trm_k20.bam" #Mode FALSE TRUE #logical 757 20334 #[1] "KOIFN_2_trm_k20.bam" #Mode FALSE TRUE #logical 922 20169 for (i in 12:15){ print (colnames(rx)[i]) tt=summary(rx[,i]>1) print(tt) } #1] "1000 * rx[, 8]/rx$genelength" #Mode FALSE TRUE #logical 8092 12999 #[1] "1000 * rx[, 9]/rx$genelength" #Mode FALSE TRUE #logical 8086 13005 #[1] "1000 * rx[, 10]/rx$genelength" #Mode FALSE TRUE #logical 8706 12385 #[1] "1000 * rx[, 11]/rx$genelength" #Mode FALSE TRUE #logical 8690 12401 #polysomeseq poly1= read.table("/mnt/max/b/genomics_facility/DKlab/mr/tzanos/mm10-GRCm38.p6/polysomeseq/featurecounts/gene_counts_merged_mm10.txt",header=T,sep="\t",skip=1) poly0 =poly1[,7:dim(poly1)[2]]; rownames(poly0)=poly1$Geneid; colnames(poly0) #[1] "KO_0h_replicate_IIb" "KO_0h_replicate_Ib" "KO_2h_replicate_IIb" "KO_2h_replicate_Ib" #[5] "KO_6h_replicate_IIb" "KO_6h_replicate_Ib" "KO_IFN.gamma_replicate_IIb" "KO_IFN_replicate_Ib" #[9] "KO_IL4_replicate_IIb" "KO_IL4_replicate_Ib" "WT_0h_replicate_IIb" "WT_0h_replicate_Ib" #[13] "WT_2h_replicate_IIb" "WT_2h_replicate_Ib" "WT_6h_replicate_IIb" "WT_6h_replicate_Ib" #[17] "WT_IFN.gamma_replicate_IIb" "WT_IFN.gamma_replicate_Ib" "WT_IL4_replicate_IIb" "WT_IL4_replicate_Ib" col_ordering = c(19,20,17,18)#WT_M1vsM2 poly =poly0[,col_ordering]; poly = poly[(rowSums(poly[,1:2])>=2)&(rowSums(poly[,3:4])>=2),] conditions = factor(c("IFN","IFN","IL4","IL4")); exp_study = DGEList(counts=poly, 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=merge(x=as.data.frame(tTags),y=gl,by="row.names",all.x=TRUE) x=x[,!(names(x) %in% "ID")] rownames(x)=x$Row.names; x=x[,!(names(x) %in% "Row.names")] x2=merge(x=x,y=c,by="row.names",all.x=TRUE) x=x2[,!(names(x2) %in% "ID")] x=cbind(x,1000*x[,8]/x$genelength,1000*x[,9]/x$genelength,1000*x[,10]/x$genelength,1000*x[,11]/x$genelength) x=cbind(x,0.5*(x[,12]+x[,13]),0.5*(x[,14]+x[,15])) write.table(x, file="edgeRgenesPolyExpr-WT_M1vsM2.csv") x_m1=x for (i in 8:11){ print (colnames(x)[i]) tt=summary(x[,i]>0) print(tt) } #[1] "WT_IL4_replicate_IIb" #Mode FALSE TRUE #logical 1379 17851 #[1] "WT_IL4_replicate_Ib" #Mode FALSE TRUE #logical 514 18716 #[1] "WT_IFN.gamma_replicate_IIb" #Mode FALSE TRUE #logical 616 18614 #[1] "WT_IFN.gamma_replicate_Ib" #Mode FALSE TRUE #logical 1220 18010 for (i in 12:15){ print (colnames(x)[i]) tt=summary(x[,i]>1) print(tt) } #[1] "1000 * x[, 8]/x$genelength" #Mode FALSE TRUE #logical 3895 15335 #[1] "1000 * x[, 9]/x$genelength" #Mode FALSE TRUE #logical 3604 15626 #[1] "1000 * x[, 10]/x$genelength" #Mode FALSE TRUE #logical 2355 16875 #[1] "1000 * x[, 11]/x$genelength" #Mode FALSE TRUE #logical 3687 15543 summary(de <- decideTestsDGE(et, p=0.05)) # IL4-IFN #Down 113 #NotSig 18957 #Up 160 col_ordering = c(9,10,7,8)#KO_M1vsM2 poly =poly0[,col_ordering]; poly = poly[(rowSums(poly[,1:2])>=2)&(rowSums(poly[,3:4])>=2),] conditions = factor(c("IFN","IFN","IL4","IL4")); exp_study = DGEList(counts=poly, 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=merge(x=as.data.frame(tTags),y=gl,by="row.names",all.x=TRUE) x=x[,!(names(x) %in% "ID")] rownames(x)=x$Row.names; x=x[,!(names(x) %in% "Row.names")] x2=merge(x=x,y=c,by="row.names",all.x=TRUE) x=x2[,!(names(x2) %in% "ID")] x=cbind(x,1000*x[,8]/x$genelength,1000*x[,9]/x$genelength,1000*x[,10]/x$genelength,1000*x[,11]/x$genelength) x=cbind(x,0.5*(x[,12]+x[,13]),0.5*(x[,14]+x[,15])) write.table(x, file="edgeRgenesPolyExpr-KO_M1vsM2.csv") x_m2=x for (i in 8:11){ print (colnames(x)[i]) tt=summary(x[,i]>0) print(tt) } #[1] "KO_IL4_replicate_IIb" #Mode FALSE TRUE #logical 1415 20451 #[1] "KO_IL4_replicate_Ib" #Mode FALSE TRUE #logical 800 21066 #[1] "KO_IFN.gamma_replicate_IIb" #Mode FALSE TRUE #logical 1012 20854 #[1] "KO_IFN_replicate_Ib" #Mode FALSE TRUE #logical 299 21567 for (i in 12:15){ print (colnames(x)[i]) tt=summary(x[,i]>1) print(tt) } #[1] "1000 * x[, 8]/x$genelength" #Mode FALSE TRUE #logical 3464 18402 #[1] "1000 * x[, 9]/x$genelength" #Mode FALSE TRUE #logical 3370 18496 #[1] "1000 * x[, 10]/x$genelength" #Mode FALSE TRUE #logical 3353 18513 #[1] "1000 * x[, 11]/x$genelength" #Mode FALSE TRUE #logical 2615 19251 #write.table(x, file="edgeRgenesExpr-IL4.csv") summary(de <- decideTestsDGE(et, p=0.05)) # IL4-IFN #Down 377 #NotSig 20981 #Up 508 # merge rnaseq and polysomeseq require(psych) all_m1=merge(rx_m1,x_m1,by="Row.names") all_m2=merge(rx_m2,x_m2,by="Row.names") > + + + + [1] "KOIL4_1_trm_k20.bam" Mode FALSE TRUE NA's logical 514 20577 0 [1] "KOIL4_2_trm_k20.bam" Mode FALSE TRUE NA's logical 518 20573 0 [1] "KOIFN_1_trm_k20.bam" Mode FALSE TRUE NA's logical 757 20334 0 [1] "KOIFN_2_trm_k20.bam" Mode FALSE TRUE NA's logical 922 20169 0 > > > > > > > > > > > > > > > + + + + [1] "1000 * rx[, 8]/rx$genelength" Mode FALSE TRUE NA's logical 8092 12999 0 [1] "1000 * rx[, 9]/rx$genelength" Mode FALSE TRUE NA's logical 8086 13005 0 [1] "1000 * rx[, 10]/rx$genelength" Mode FALSE TRUE NA's logical 8706 12385 0 [1] "1000 * rx[, 11]/rx$genelength" Mode FALSE TRUE NA's logical 8690 12401 0 > > > > > > > > > > > > > > > > > > > > [1] "KO_0h_replicate_IIb" "KO_0h_replicate_Ib" [3] "KO_2h_replicate_IIb" "KO_2h_replicate_Ib" [5] "KO_6h_replicate_IIb" "KO_6h_replicate_Ib" [7] "KO_IFN.gamma_replicate_IIb" "KO_IFN_replicate_Ib" [9] "KO_IL4_replicate_IIb" "KO_IL4_replicate_Ib" [11] "WT_0h_replicate_IIb" "WT_0h_replicate_Ib" [13] "WT_2h_replicate_IIb" "WT_2h_replicate_Ib" [15] "WT_6h_replicate_IIb" "WT_6h_replicate_Ib" [17] "WT_IFN.gamma_replicate_IIb" "WT_IFN.gamma_replicate_Ib" [19] "WT_IL4_replicate_IIb" "WT_IL4_replicate_Ib" > > > > > > > > > > > > > > > > > > > > > > > > > > + + + + [1] "WT_IL4_replicate_IIb" Mode FALSE TRUE NA's logical 1379 17851 0 [1] "WT_IL4_replicate_Ib" Mode FALSE TRUE NA's logical 514 18716 0 [1] "WT_IFN.gamma_replicate_IIb" Mode FALSE TRUE NA's logical 616 18614 0 [1] "WT_IFN.gamma_replicate_Ib" Mode FALSE TRUE NA's logical 1220 18010 0 > > > > > > > > > > > > > > > + + + + [1] "1000 * x[, 8]/x$genelength" Mode FALSE TRUE NA's logical 3895 15335 0 [1] "1000 * x[, 9]/x$genelength" Mode FALSE TRUE NA's logical 3604 15626 0 [1] "1000 * x[, 10]/x$genelength" Mode FALSE TRUE NA's logical 2355 16875 0 [1] "1000 * x[, 11]/x$genelength" Mode FALSE TRUE NA's logical 3687 15543 0 > > > > > > > > > > > > > > > [,1] -1 113 0 18957 1 160 > > > > > > > > > > > > > > > > > > > > > > > > + + + + [1] "KO_IL4_replicate_IIb" Mode FALSE TRUE NA's logical 1415 20451 0 [1] "KO_IL4_replicate_Ib" Mode FALSE TRUE NA's logical 800 21066 0 [1] "KO_IFN.gamma_replicate_IIb" Mode FALSE TRUE NA's logical 1012 20854 0 [1] "KO_IFN_replicate_Ib" Mode FALSE TRUE NA's logical 299 21567 0 > > > > > > > > > > > > > > > + + + + [1] "1000 * x[, 8]/x$genelength" Mode FALSE TRUE NA's logical 3464 18402 0 [1] "1000 * x[, 9]/x$genelength" Mode FALSE TRUE NA's logical 3370 18496 0 [1] "1000 * x[, 10]/x$genelength" Mode FALSE TRUE NA's logical 3353 18513 0 [1] "1000 * x[, 11]/x$genelength" Mode FALSE TRUE NA's logical 2615 19251 0 > > > > > > > > > > > > > > > > [,1] -1 377 0 20981 1 508 > > > > > > > > Loading required package: psych > > str(all_m1) 'data.frame': 14671 obs. of 33 variables: $ Row.names :Class 'AsIs' chr [1:14671] "ENSMUSG00000000001" "ENSMUSG00000000028" "ENSMUSG00000000031" "ENSMUSG00000000037" ... $ logFC.x : num 1.02 -1.55 -1.66 -1.81 -2.58 ... $ logCPM.x : num 7.508 2.638 -3.576 -0.569 4.12 ... $ PValue.x : num 1.06e-14 8.00e-13 1.95e-01 3.05e-06 9.20e-41 ... $ FDR.x : num 7.44e-14 4.94e-12 2.93e-01 1.06e-05 2.03e-39 ... $ geneID.x : Factor w/ 45583 levels " 0610005C13Rik",..: 26497 5087 26921 39501 31916 4679 28785 39499 5892 8318 ... $ genelength.x : num 3253 1560 1265 3276 1780 ... $ wtIL4_1_trm_k20.bam : num 119.725 9.01 0.144 0.906 30.019 ... $ wtIL4_2_trm_k20.bam : num 121.185 9.4 0.023 1.067 29.312 ... $ wtIFN_1_trm_k20.bam : num 233.0333 3.1644 0.0301 0.211 4.7164 ... $ wtIFN_2_trm_k20.bam : num 254.1179 3.1077 0.0226 0.339 5.2209 ... $ 1000 * rx[, 8]/rx$genelength : num 36.804 5.777 0.114 0.277 16.868 ... $ 1000 * rx[, 9]/rx$genelength : num 37.2534 6.0267 0.0182 0.3259 16.4706 ... $ 1000 * rx[, 10]/rx$genelength: num 71.6364 2.0289 0.0238 0.0644 2.6502 ... $ 1000 * rx[, 11]/rx$genelength: num 78.118 1.9925 0.0179 0.1035 2.9336 ... $ 0.5 * (rx[, 12] + rx[, 13]) : num 37.0289 5.9018 0.0658 0.3013 16.6691 ... $ 0.5 * (rx[, 14] + rx[, 15]) : num 74.8772 2.0107 0.0209 0.084 2.7919 ... $ logFC.y : num -0.04481 -0.33237 0.91204 -0.00543 -1.56371 ... $ logCPM.y : num 6.52 2.45 4.85 3.29 3.41 ... $ PValue.y : num 0.9597 1 0.2025 1 0.0818 ... $ FDR.y : num 1 1 0.854 1 0.619 ... $ geneID.y : Factor w/ 45583 levels " 0610005C13Rik",..: 26497 5087 26921 39501 31916 4679 28785 39499 5892 8318 ... $ genelength.y : num 3253 1560 1265 3276 1780 ... $ WT_IL4_replicate_IIb : num 80.08 2.81 15.45 7.02 14.05 ... $ WT_IL4_replicate_Ib : num 100.87 3.81 20.94 7.61 10.47 ... $ WT_IFN.gamma_replicate_IIb : num 70.6 2.66 54.62 10.66 6.66 ... $ WT_IFN.gamma_replicate_Ib : num 105.1 2.66 14.63 3.99 1.33 ... $ 1000 * x[, 8]/x$genelength : num 24.62 1.8 12.22 2.14 7.89 ... $ 1000 * x[, 9]/x$genelength : num 31.01 2.44 16.55 2.32 5.88 ... $ 1000 * x[, 10]/x$genelength : num 21.7 1.71 43.19 3.25 3.74 ... $ 1000 * x[, 11]/x$genelength : num 32.31 1.706 11.573 1.218 0.748 ... $ 0.5 * (x[, 12] + x[, 13]) : num 27.81 2.12 14.39 2.23 6.89 ... $ 0.5 * (x[, 14] + x[, 15]) : num 27.01 1.71 27.38 2.24 2.25 ... > getwd() [1] "/mnt/max/b/genomics_facility/DKlab/mr/rnaseq-vs-polysomeseq" > png("deregulogram_m1_rna_vs_poly.png",width=1024,height=1024,pointsize = 22) plot(all_m1$logFC.x,all_m1$logFC.y) dev.off() > null device 1 > translatome_only=subset(all_m1,(abs(all_m1$logFC.x)>=1.8)&(abs(all_m1$logFC.y)<1.8)); transcriptiome_only=subset(all_m1,(abs(all_m1$logFC.y)>=1.8)&(abs(all_m1$logFC.x)<1.8)); png("deregulogram_m1_rna_vs_poly.png",width=1024,height=1024,pointsize = 22) plot(all_m1$logFC.x,all_m1$logFC.y) points(translatome_only$logFC.x,translatome_only$logFC.y,col="yellow") points(transcriptome_only$logFC.x,transcriptome_only$logFC.y,col="blue") dev.off() > Error in points(transcriptome_only$logFC.x, transcriptome_only$logFC.y, : object 'transcriptome_only' not found > > null device 1 > transcriptome_only=subset(all_m1,(abs(all_m1$logFC.y)>=1.8)&(abs(all_m1$logFC.x)<1.8)); png("deregulogram_m1_rna_vs_poly.png",width=1024,height=1024,pointsize = 22) plot(all_m1$logFC.x,all_m1$logFC.y) points(translatome_only$logFC.x,translatome_only$logFC.y,col="yellow") points(transcriptome_only$logFC.x,transcriptome_only$logFC.y,col="blue") dev.off() > null device 1 > translatome_only=subset(all_m1,(abs(all_m1$logFC.x)>=1.8)&(abs(all_m1$logFC.y)<1.8)); transcriptome_only=subset(all_m1,(abs(all_m1$logFC.y)>=1.8)&(abs(all_m1$logFC.x)<1.8)); homodirectional=subset(all_m1,((all_m1$logFC.x)>=1.8)&((all_m1$logFC.y)>=1.8)); #add down+down case antidirectional=subset(all_m1,((all_m1$logFC.x)>=1.8)&((all_m1$logFC.y)<1.8)); # add down+up case png("deregulogram_m1_rna_vs_poly.png",width=1024,height=1024,pointsize = 22) plot(all_m1$logFC.x,all_m1$logFC.y) points(translatome_only$logFC.x,translatome_only$logFC.y,col="yellow") points(transcriptome_only$logFC.x,transcriptome_only$logFC.y,col="blue") points(homodirectional$logFC.x,homodirectional$logFC.y,col="green") points(antidirectional$logFC.x,antidirectional$logFC.y,col="red") dev.off() > null device 1 > translatome_only=subset(all_m1,(abs(all_m1$logFC.y)>=1.8)&(abs(all_m1$logFC.x)<1.8)); transcriptome_only=subset(all_m1,(abs(all_m1$logFC.x)>=1.8)&(abs(all_m1$logFC.y)<1.8)); homodirectional=subset(all_m1,((all_m1$logFC.x)>=1.8)&((all_m1$logFC.y)>=1.8)); #add down+down case antidirectional=subset(all_m1,((all_m1$logFC.x)>=1.8)&((all_m1$logFC.y)<-1.8)); # add down+up case png("deregulogram_m1_rna_vs_poly.png",width=1024,height=1024,pointsize = 22) plot(all_m1$logFC.x,all_m1$logFC.y) points(translatome_only$logFC.x,translatome_only$logFC.y,col="yellow") points(transcriptome_only$logFC.x,transcriptome_only$logFC.y,col="blue") points(homodirectional$logFC.x,homodirectional$logFC.y,col="green") points(antidirectional$logFC.x,antidirectional$logFC.y,col="red") dev.off() #write.csv(all_m1, file="all_WT.csv",row.names=F) > > > Error in (all_m1$logFC.y) <- 1.8 : could not find function "(<-" > > > > > > > > > null device 1 > translatome_only=subset(all_m1,(abs(all_m1$logFC.y)>=1.8)&(abs(all_m1$logFC.x)<1.8)); transcriptome_only=subset(all_m1,(abs(all_m1$logFC.x)>=1.8)&(abs(all_m1$logFC.y)<1.8)); homodirectional=subset(all_m1,((all_m1$logFC.x)>=1.8)&((all_m1$logFC.y)>=1.8)); #add down+down case antidirectional=subset(all_m1,((all_m1$logFC.x)>=1.8)&((all_m1$logFC.y)<(-1.8))); # add down+up case png("deregulogram_m1_rna_vs_poly.png",width=1024,height=1024,pointsize = 22) plot(all_m1$logFC.x,all_m1$logFC.y) points(translatome_only$logFC.x,translatome_only$logFC.y,col="yellow") points(transcriptome_only$logFC.x,transcriptome_only$logFC.y,col="blue") points(homodirectional$logFC.x,homodirectional$logFC.y,col="green") points(antidirectional$logFC.x,antidirectional$logFC.y,col="red") dev.off() > null device 1 > getwd(() Error: unexpected ')' in "getwd(()" > getwd() [1] "/mnt/max/b/genomics_facility/DKlab/mr/rnaseq-vs-polysomeseq" > setwd("../tzanos/mm10-GRCm38.p6/rnaseq-vs-polysomeseq/") > translatome_only=subset(all_m1,(abs(all_m1$logFC.y)>=1.8)&(abs(all_m1$logFC.x)<1.8)); transcriptome_only=subset(all_m1,(abs(all_m1$logFC.x)>=1.8)&(abs(all_m1$logFC.y)<1.8)); homodirectional=subset(all_m1,((all_m1$logFC.x)>=1.8)&((all_m1$logFC.y)>=1.8)); #add down+down case antidirectional=subset(all_m1,((all_m1$logFC.x)>=1.8)&((all_m1$logFC.y)<(-1.8))); # add down+up case png("deregulogram_m1_rna_vs_poly.png",width=1024,height=1024,pointsize = 22) plot(all_m1$logFC.x,all_m1$logFC.y) points(translatome_only$logFC.x,translatome_only$logFC.y,col="yellow") points(transcriptome_only$logFC.x,transcriptome_only$logFC.y,col="blue") points(homodirectional$logFC.x,homodirectional$logFC.y,col="green") points(antidirectional$logFC.x,antidirectional$logFC.y,col="red") dev.off() > null device 1 > > head (tTags) Error: Two subscripts required > head (tTags,5) Error: Two subscripts required > str(tTags) Formal class 'TopTags' [package "edgeR"] with 1 slot ..@ .Data:List of 4 .. ..$ :'data.frame': 21866 obs. of 4 variables: .. .. ..$ logFC : num [1:21866] 5.81 6.23 5.53 6.33 6.47 ... .. .. ..$ logCPM: num [1:21866] 11.23 9.66 9.36 7.18 7.94 ... .. .. ..$ PValue: num [1:21866] 1.74e-40 2.75e-40 2.56e-34 3.14e-28 1.04e-26 ... .. .. ..$ FDR : num [1:21866] 3.00e-36 3.00e-36 1.87e-30 1.72e-24 4.56e-23 ... .. ..$ : chr "fdr" .. ..$ : chr [1:2] "IFN" "IL4" .. ..$ : chr "exact" > head(rownames(poly)) [1] "ENSMUSG00000051951" "ENSMUSG00000103377" "ENSMUSG00000104017" [4] "ENSMUSG00000103025" "ENSMUSG00000103201" "ENSMUSG00000103161" > poly[rownames(poly)=="ENSMUSG00000051951",] KO_IL4_replicate_IIb KO_IL4_replicate_Ib ENSMUSG00000051951 5 5 KO_IFN.gamma_replicate_IIb KO_IFN_replicate_Ib ENSMUSG00000051951 9 10 > head (r) [1] "ENSMUSG00000020641" "ENSMUSG00000022126" "ENSMUSG00000020638" [4] "ENSMUSG00000030156" "ENSMUSG00000041827" "ENSMUSG00000038179" > poly[rownames(poly)=="ENSMUSG00000020641",] KO_IL4_replicate_IIb KO_IL4_replicate_Ib ENSMUSG00000020641 47 81 KO_IFN.gamma_replicate_IIb KO_IFN_replicate_Ib ENSMUSG00000020641 4046 6115 > head(tTags) Error: Two subscripts required > tTags Comparison of groups: IL4-IFN logFC logCPM PValue FDR ENSMUSG00000020641 5.8099539420 11.226557 1.736083e-40 3.004988e-36 ENSMUSG00000022126 6.2315874981 9.663629 2.748548e-40 3.004988e-36 ENSMUSG00000020638 5.5310326460 9.361139 2.561904e-34 1.867286e-30 ENSMUSG00000030156 6.3335683618 7.180348 3.142576e-28 1.717889e-24 > 2^5.8 [1] 55.71524 > 4046/47 [1] 86.08511 > 6115/81 [1] 75.49383 > exp_study$samples$group [1] IFN IFN IL4 IL4 Levels: IFN IL4 > group=as.factor(exp_study$samples$group) > group [1] IFN IFN IL4 IL4 Levels: IFN IL4 > pair=1:2 > levs.group <- levels(group) if(is.numeric(pair)) else pair <- as.character(pair) > + > Error: unexpected 'else' in "else" > > > pair <- levs.group[pair] > pair [1] NA NA > levs.group [1] "IFN" "IL4" > pair=1:2 > pair <- levs.group[pair] > pair [1] "IFN" "IL4" > j <- group %in% pair > j [1] TRUE TRUE TRUE TRUE > j1 <- group==pair[1] > > j1 [1] TRUE TRUE FALSE FALSE > str(exp_study$counts) int [1:21866, 1:4] 5 10 3 2 1 0 2 2 11 2 ... - attr(*, "dimnames")=List of 2 ..$ : chr [1:21866] "ENSMUSG00000051951" "ENSMUSG00000103377" "ENSMUSG00000104017" "ENSMUSG00000103025" ... ..$ : chr [1:4] "KO_IL4_replicate_IIb" "KO_IL4_replicate_Ib" "KO_IFN.gamma_replicate_IIb" "KO_IFN_replicate_Ib" > getwd() [1] "/mnt/max/b/genomics_facility/DKlab/mr/tzanos/mm10-GRCm38.p6/rnaseq-vs-polysomeseq" > conditions = factor(c("IL4","IL4","IFN","IFN")); exp_study = DGEList(counts=poly, group=conditions);exp_study = calcNormFactors(exp_study);exp_study = estimateCommonDisp(exp_study);exp_study = estimateTagwiseDisp(exp_study);et = exactTest(exp_study,pair=c("IFN","IL4")) tTags = topTags(et,n=NULL,adjust.method = "fdr") > Comparison of groups: IL4-IFN logFC logCPM PValue FDR ENSMUSG00000020641 -5.8099539420 11.226557 1.736083e-40 3.004988e-36 ENSMUSG00000022126 -6.2315874981 9.663629 2.748548e-40 3.004988e-36 ENSMUSG00000020638 -5.5310326460 9.361139 2.561904e-34 1.867286e-30