/data/results/tools/rnaseq/subread/subread-1.5.2-source/bin/featureCounts -O --fraction -M -s 1 -J -f -g gene_name -t exon -T 4 -a /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.gtf -o wt01.cnt wt01.bam wt02.bam WT_0h_replicate_I.bam WT_0h_replicate_II_.bam awk -f getRPM1.awk wt01.cnt > wt01.cnt.rpm.csv awk -f getRPMperGene1.awk wt01.cnt > wt01.cnt.rpmPerGene.csv Sample Total_mapped_reads wt01.bam 3.69769e+07 wt02.bam 3.72366e+07 WT_0h_replicate_I.bam 535162 WT_0h_replicate_II_.bam 382403 /data/results/tools/rnaseq/subread/subread-1.5.2-source/bin/featureCounts -O --fraction -M -s 1 -J -f -g gene_name -t exon -T 4 -a /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.gtf -o wt21.cnt wt21.bam wt22.bam WT_2h_replicate_I.bam WT_2h_replicate_II.bam awk -f getRPM1.awk wt21.cnt > wt21.cnt.rpm.csv awk -f getRPMperGene1.awk wt21.cnt > wt21.cnt.rpmPerGene.csv Sample Total_mapped_reads wt21.bam 3.46494e+07 wt22.bam 4.00449e+07 WT_2h_replicate_I.bam 501228 WT_2h_replicate_II.bam 279147 /data/results/tools/rnaseq/subread/subread-1.5.2-source/bin/featureCounts -O --fraction -M -s 1 -J -f -g gene_name -t exon -T 4 -a /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.gtf -o wt61.cnt wt61.bam wt62.bam WT_6h_replicate_I.bam WT_6h_replicate_II.bam awk -f getRPM1.awk wt61.cnt > wt61.cnt.rpm.csv awk -f getRPMperGene1.awk wt61.cnt > wt61.cnt.rpmPerGene.csv Sample Total_mapped_reads wt61.bam 3.95097e+07 wt62.bam 3.8493e+07 WT_6h_replicate_I.bam 387418 WT_6h_replicate_II.bam 1.05731e+06 /data/results/tools/rnaseq/subread/subread-1.5.2-source/bin/featureCounts -O --fraction -M -s 1 -J -f -g gene_name -t exon -T 4 -a /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.gtf -o wtIFN1.cnt wtIFN1.bam wtIFN2.bam WT_IFN-gamma_replicate_I.bam WT_IFN-gamma_replicate_II.bam awk -f getRPM1.awk wtIFN1.cnt > wtIFN1.cnt.rpm.csv awk -f getRPMperGene1.awk wtIFN1.cnt > wtIFN1.cnt.rpmPerGene.csv Sample Total_mapped_reads wtIFN1.bam 3.9288e+07 wtIFN2.bam 5.25123e+07 WT_IFN-gamma_replicate_I.bam 1.16697e+06 WT_IFN-gamma_replicate_II.bam 544481 /data/results/tools/rnaseq/subread/subread-1.5.2-source/bin/featureCounts -O --fraction -M -s 1 -J -f -g gene_name -t exon -T 4 -a /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.gtf -o wt01IL4.cnt wt01IL4.bam wt02IL4.bam WT_IL4_replicate_I.bam WT_IL4_replicate_II.bam awk -f getRPM1.awk wt01IL4.cnt > wt01IL4.cnt.rpm.csv awk -f getRPMperGene1.awk wt01IL4.cnt > wt01IL4.cnt.rpmPerGene.csv Sample Total_mapped_reads wt01IL4.bam 4.57558e+07 wt02IL4.bam 3.52183e+07 WT_IL4_replicate_I.bam 1.44613e+06 WT_IL4_replicate_II.bam 829652 /data/results/tools/rnaseq/subread/subread-1.5.2-source/bin/featureCounts -O --fraction -M -s 1 -J -f -g gene_name -t exon -T 4 -a /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.gtf -o Ko01.cnt Ko01.bam Ko02.bam KO_0h_replicate_I.bam KO_0h_replicate_II.bam awk -f getRPM1.awk Ko01.cnt > Ko01.cnt.rpm.csv awk -f getRPMperGene1.awk Ko01.cnt > Ko01.cnt.rpmPerGene.csv Sample Total_mapped_reads Ko01.bam 4.01576e+07 Ko02.bam 5.13778e+07 KO_0h_replicate_I.bam 825450 KO_0h_replicate_II.bam 1.01729e+06 /data/results/tools/rnaseq/subread/subread-1.5.2-source/bin/featureCounts -O --fraction -M -s 1 -J -f -g gene_name -t exon -T 4 -a /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.gtf -o Ko21.cnt Ko21.bam Ko22.bam KO_2h_replicate_I.bam KO_2h_replicate_II.bam awk -f getRPM1.awk Ko21.cnt > Ko21.cnt.rpm.csv awk -f getRPMperGene1.awk Ko21.cnt > Ko21.cnt.rpmPerGene.csv Sample Total_mapped_reads Ko21.bam 4.89587e+07 Ko22.bam 4.84856e+07 KO_2h_replicate_I.bam 625229 KO_2h_replicate_II.bam 522928 /data/results/tools/rnaseq/subread/subread-1.5.2-source/bin/featureCounts -O --fraction -M -s 1 -J -f -g gene_name -t exon -T 4 -a /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.gtf -o Ko61.cnt Ko61.bam Ko62.bam KO_6h_replicate_I.bam KO_6h_replicate_II.bam awk -f getRPM1.awk Ko61.cnt > Ko61.cnt.rpm.csv awk -f getRPMperGene1.awk Ko61.cnt > Ko61.cnt.rpmPerGene.csv Sample Total_mapped_reads Ko61.bam 4.64522e+07 Ko62.bam 5.48543e+07 KO_6h_replicate_I.bam 1.13254e+06 KO_6h_replicate_II.bam 752620 /data/results/tools/rnaseq/subread/subread-1.5.2-source/bin/featureCounts -O --fraction -M -s 1 -J -f -g gene_name -t exon -T 4 -a /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.gtf -o KOIFN1.cnt KoIFN1.bam KoIFN2.bam KO_IFN-gamma_replicate_I.bam KO_IFN-gamma_replicate_II.bam awk -f getRPM1.awk KOIFN1.cnt > KOIFN1.cnt.rpm.csv awk -f getRPMperGene1.awk KOIFN1.cnt > KOIFN1.cnt.rpmPerGene.csv Sample Total_mapped_reads KoIFN1.bam 5.11802e+07 KoIFN2.bam 4.43768e+07 KO_IFN-gamma_replicate_I.bam 1.26831e+06 KO_IFN-gamma_replicate_II.bam 908339 /data/results/tools/rnaseq/subread/subread-1.5.2-source/bin/featureCounts -O --fraction -M -s 1 -J -f -g gene_name -t exon -T 4 -a /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.gtf -o Ko1IL4.cnt Ko1IL4.bam Ko2IL4.bam KO_IL4_replicate_I.bam KO_IL4_replicate_II.bam awk -f getRPM1.awk Ko1IL4.cnt > Ko1IL4.cnt.rpm.csv awk -f getRPMperGene1.awk Ko1IL4.cnt > Ko1IL4.cnt.rpmPerGene.csv Sample Total_mapped_reads Ko1IL4.bam 5.02855e+07 Ko2IL4.bam 5.17752e+07 KO_IL4_replicate_I.bam 809723 KO_IL4_replicate_II.bam 629921 #@ zcat mm9_merged_exons_as_genes.txt.gz | head chromosome start end gene_id gc_content strand gene_name biotype chr3 107910198 107912234 ENSMUSG00000000001_MEX_1 0 - Gnai3 protein_coding chr3 107912321 107912530 ENSMUSG00000000001_MEX_2 0 - Gnai3 protein_coding chromosome start end transcript_id gene_id strand gene_name biotype chr1 3073823 3075077 ENSMUSG00000102693_MET_1 ENSMUSG00000102693 + 4933401J01Rik TEC chr1 3101626 3102880 ENSMUSG00000064842_MET_1 ENSMUSG00000064842 + Gm26206 snRNA reczko@max:/data/images/proton/DKlab/orsalia/perExon$ zcat /data/results/tools/rnaseq/metaseqr/mm10/gene_data_mm10.txt.gz | head chromosome start end gene_id gc_content strand gene_name biotype chr1 3073253 3074322 ENSMUSG00000102693 34.21 + 4933401J01Rik TEC chr1 3102016 3102125 ENSMUSG00000064842 36.36 + Gm26206 snRNA chr1 3205901 3671498 ENSMUSG00000051951 38.51 - Xkr4 protein_coding chr1 3252757 3253236 ENSMUSG00000102851 39.79 + Gm18956 processed_pseudogene chr1 3365731 3368549 ENSMUSG00000103377 40.79 - Gm37180 TEC chr1 3375556 3377788 ENSMUSG00000104017 36.99 - Gm37363 TEC reczko@max:/data/images/proton/DKlab/orsalia/perExon$ zcat mm9_merged_exons_as_genes.txt.gz |awk -f exon2geneFMT.awk > mm9_merged_exons_as_genes2.txt /data/results/tools/rnaseq/subread/subread-1.5.2-source/bin/featureCounts -O --fraction -M -s 1 -J -f -g gene_name -t exon -T 4 -a /data/results/reference/hsa/gencode/gencode.v26lift37.annotation_exons.gtf -o /data/images/proton2/run341/rpkm/exonGencode.cnt PHR23r-FUBP2.bam PHR22r-TCF4.bam PHR21r-bcatMono.bam PHR20r-bcatPoly.bam PHR24r-IgG.bam PHR19r-Input.bam /data/images/proton/run15/www/con_b.bam /data/images/proton/run15/www/con_a.bam cd /data/images/proton2/run341/rpkm/all /data/images/proton2/run341/rpkm/getRPKM_with_control.awk awk -f ../getRPKM_with_control.awk ../exonGencode.cnt > getRPKM_with_control.log #edit header /data/images/proton2/run341/rpkm/all/process-RPKM2.r library(hexbin) t=read.table("/data/images/proton/DKlab/orsalia/perExon/wt01.cnt.rpm.csv",header=T) png("wt0h-PolysomeVSRNAseq.png",width=1024,height=1024,pointsize = 22) plot (log2(0.5*(t$wt01.bam_rpm+t$wt02.bam_rpm)),log2(0.5*(t$WT_0h_replicate_I.bam_rpm+t$WT_0h_replicate_II_.bam_rpm)),pch=".",xlab="reads_per_exon_per_million_RNAseq",ylab="reads_per_exon_per_million_PolysomesSeq") dev.off() x=log2(0.5*(t$wt01.bam_rpm+t$wt02.bam_rpm));y=log2(0.5*(t$WT_0h_replicate_I.bam_rpm+t$WT_0h_replicate_II_.bam_rpm)); x[mapply(is.infinite, x)] <- NA;y[mapply(is.infinite, y)] <- NA; bin<-hexbin(x,y, xbins=50);bin@count=log2(bin@count); png("wt0h-PolysomeVSRNAseq-hist.png",width=1024,height=1024,pointsize = 22) plot(bin, main="wt0h",xlab="reads_per_exon_per_million_RNAseq",ylab="reads_per_exon_per_million_PolysomesSeq") dev.off() t=read.table("/data/images/proton/DKlab/orsalia/perExon/wt21.cnt.rpm.csv",header=T) png("wt2h-PolysomeVSRNAseq.png",width=1024,height=1024,pointsize = 22) plot (log2(0.5*(t$wt21.bam_rpm+t$wt22.bam_rpm)),log2(0.5*(t$WT_2h_replicate_I.bam_rpm+t$WT_2h_replicate_II.bam_rpm)),pch=".",xlab="reads_per_exon_per_million_RNAseq",ylab="reads_per_exon_per_million_PolysomesSeq") dev.off() x=log2(0.5*(t$wt21.bam_rpm+t$wt22.bam_rpm));y=log2(0.5*(t$WT_2h_replicate_I.bam_rpm+t$WT_2h_replicate_II.bam_rpm)); x[mapply(is.infinite, x)] <- NA;y[mapply(is.infinite, y)] <- NA; bin<-hexbin(x,y, xbins=50);bin@count=log2(bin@count); png("wt2h-PolysomeVSRNAseq-hist.png",width=1024,height=1024,pointsize = 22) plot(bin, main="wt2h",xlab="reads_per_exon_per_million_RNAseq",ylab="reads_per_exon_per_million_PolysomesSeq") dev.off() t=read.table("/data/images/proton/DKlab/orsalia/perExon/wt61.cnt.rpm.csv",header=T) png("wt6h-PolysomeVSRNAseq.png",width=1024,height=1024,pointsize = 22) plot (log2(0.5*(t$wt61.bam_rpm+t$wt62.bam_rpm)),log2(0.5*(t$WT_6h_replicate_I.bam_rpm+t$WT_6h_replicate_II.bam_rpm)),pch=".",xlab="reads_per_exon_per_million_RNAseq",ylab="reads_per_exon_per_million_PolysomesSeq") dev.off() x=log2(0.5*(t$wt61.bam_rpm+t$wt62.bam_rpm));y=log2(0.5*(t$WT_6h_replicate_I.bam_rpm+t$WT_6h_replicate_II.bam_rpm)); x[mapply(is.infinite, x)] <- NA;y[mapply(is.infinite, y)] <- NA; bin<-hexbin(x,y, xbins=50);bin@count=log2(bin@count); png("wt6h-PolysomeVSRNAseq-hist.png",width=1024,height=1024,pointsize = 22) plot(bin, main="wt6h",xlab="reads_per_exon_per_million_RNAseq",ylab="reads_per_exon_per_million_PolysomesSeq") dev.off() t=read.table("/data/images/proton/DKlab/orsalia/perExon/wtIFN1.cnt.rpm.csv",header=T) png("wtIFN-PolysomeVSRNAseq.png",width=1024,height=1024,pointsize = 22) plot (log2(0.5*(t$wtIFN1.bam_rpm+t$wtIFN2.bam_rpm)),log2(0.5*(t$WT_IFN.gamma_replicate_I.bam_rpm+t$WT_IFN.gamma_replicate_II.bam_rpm)),pch=".",xlab="reads_per_exon_per_million_RNAseq",ylab="reads_per_exon_per_million_PolysomesSeq") dev.off() x=log2(0.5*(t$wtIFN1.bam_rpm+t$wtIFN2.bam_rpm));y=log2(0.5*(t$WT_IFN.gamma_replicate_I.bam_rpm+t$WT_IFN.gamma_replicate_II.bam_rpm)); x[mapply(is.infinite, x)] <- NA;y[mapply(is.infinite, y)] <- NA; bin<-hexbin(x,y, xbins=50);bin@count=log2(bin@count); png("wtIFN-PolysomeVSRNAseq-hist.png",width=1024,height=1024,pointsize = 22) plot(bin, main="wtIFN",xlab="reads_per_exon_per_million_RNAseq",ylab="reads_per_exon_per_million_PolysomesSeq") dev.off() t=read.table("/data/images/proton/DKlab/orsalia/perExon/wt01IL4.cnt.rpm.csv",header=T) png("wtIL4-PolysomeVSRNAseq.png",width=1024,height=1024,pointsize = 22) plot (log2(0.5*(t$wt01IL4.bam_rpm+t$wt02IL4.bam_rpm)),log2(0.5*(t$WT_IL4_replicate_I.bam_rpm+t$WT_IL4_replicate_II.bam_rpm)),pch=".",xlab="reads_per_exon_per_million_RNAseq",ylab="reads_per_exon_per_million_PolysomesSeq") dev.off() x=log2(0.5*(t$wt01IL4.bam_rpm+t$wt02IL4.bam_rpm));y=log2(0.5*(t$WT_IL4_replicate_I.bam_rpm+t$WT_IL4_replicate_II.bam_rpm)); x[mapply(is.infinite, x)] <- NA;y[mapply(is.infinite, y)] <- NA; bin<-hexbin(x,y, xbins=50);bin@count=log2(bin@count); png("wtIL4-PolysomeVSRNAseq-hist.png",width=1024,height=1024,pointsize = 22) plot(bin, main="wtIL4",xlab="reads_per_exon_per_million_RNAseq",ylab="reads_per_exon_per_million_PolysomesSeq") dev.off() t=read.table("/data/images/proton/DKlab/orsalia/perExon/Ko01.cnt.rpm.csv",header=T) png("Ko0h-PolysomeVSRNAseq.png",width=1024,height=1024,pointsize = 22) plot (log2(0.5*(t$Ko01.bam_rpm+t$Ko02.bam_rpm)),log2(0.5*(t$KO_0h_replicate_I.bam_rpm+t$KO_0h_replicate_II.bam_rpm)),pch=".",xlab="reads_per_exon_per_million_RNAseq",ylab="reads_per_exon_per_million_PolysomesSeq") dev.off() x=log2(0.5*(t$Ko01.bam_rpm+t$Ko02.bam_rpm));y=log2(0.5*(t$KO_0h_replicate_I.bam_rpm+t$KO_0h_replicate_II.bam_rpm)); x[mapply(is.infinite, x)] <- NA;y[mapply(is.infinite, y)] <- NA; bin<-hexbin(x,y, xbins=50);bin@count=log2(bin@count); png("Ko0h-PolysomeVSRNAseq-hist.png",width=1024,height=1024,pointsize = 22) plot(bin, main="Ko0h",xlab="reads_per_exon_per_million_RNAseq",ylab="reads_per_exon_per_million_PolysomesSeq") dev.off() t=read.table("/data/images/proton/DKlab/orsalia/perExon/Ko21.cnt.rpm.csv",header=T) png("Ko2h-PolysomeVSRNAseq.png",width=1024,height=1024,pointsize = 22) plot (log2(0.5*(t$Ko21.bam_rpm+t$Ko22.bam_rpm)),log2(0.5*(t$KO_2h_replicate_I.bam_rpm+t$KO_2h_replicate_II.bam_rpm)),pch=".",xlab="reads_per_exon_per_million_RNAseq",ylab="reads_per_exon_per_million_PolysomesSeq") dev.off() x=log2(0.5*(t$Ko21.bam_rpm+t$Ko22.bam_rpm));y=log2(0.5*(t$KO_2h_replicate_I.bam_rpm+t$KO_2h_replicate_II.bam_rpm)); x[mapply(is.infinite, x)] <- NA;y[mapply(is.infinite, y)] <- NA; bin<-hexbin(x,y, xbins=50);bin@count=log2(bin@count); png("Ko2h-PolysomeVSRNAseq-hist.png",width=1024,height=1024,pointsize = 22) plot(bin, main="Ko2h",xlab="reads_per_exon_per_million_RNAseq",ylab="reads_per_exon_per_million_PolysomesSeq") dev.off() t=read.table("/data/images/proton/DKlab/orsalia/perExon/Ko61.cnt.rpm.csv",header=T) png("Ko6h-PolysomeVSRNAseq.png",width=1024,height=1024,pointsize = 22) plot (log2(0.5*(t$Ko61.bam_rpm+t$Ko62.bam_rpm)),log2(0.5*(t$KO_6h_replicate_I.bam_rpm+t$KO_6h_replicate_II.bam_rpm)),pch=".",xlab="reads_per_exon_per_million_RNAseq",ylab="reads_per_exon_per_million_PolysomesSeq") dev.off() x=log2(0.5*(t$Ko61.bam_rpm+t$Ko62.bam_rpm));y=log2(0.5*(t$KO_6h_replicate_I.bam_rpm+t$KO_6h_replicate_II.bam_rpm)); x[mapply(is.infinite, x)] <- NA;y[mapply(is.infinite, y)] <- NA; bin<-hexbin(x,y, xbins=50);bin@count=log2(bin@count); png("Ko6h-PolysomeVSRNAseq-hist.png",width=1024,height=1024,pointsize = 22) plot(bin, main="Ko6h",xlab="reads_per_exon_per_million_RNAseq",ylab="reads_per_exon_per_million_PolysomesSeq") dev.off() t=read.table("/data/images/proton/DKlab/orsalia/perExon/KOIFN1.cnt.rpm.csv",header=T) png("KOIFN-PolysomeVSRNAseq.png",width=1024,height=1024,pointsize = 22) plot (log2(0.5*(t$KoIFN1.bam_rpm+t$KoIFN2.bam_rpm)),log2(0.5*(t$KO_IFN.gamma_replicate_I.bam_rpm+t$KO_IFN.gamma_replicate_II.bam_rpm)),pch=".",xlab="reads_per_exon_per_million_RNAseq",ylab="reads_per_exon_per_million_PolysomesSeq") dev.off() x=log2(0.5*(t$KoIFN1.bam_rpm+t$KoIFN2.bam_rpm));y=log2(0.5*(t$KO_IFN.gamma_replicate_I.bam_rpm+t$KO_IFN.gamma_replicate_II.bam_rpm)); x[mapply(is.infinite, x)] <- NA;y[mapply(is.infinite, y)] <- NA; bin<-hexbin(x,y, xbins=50);bin@count=log2(bin@count); png("KOIFN-PolysomeVSRNAseq-hist.png",width=1024,height=1024,pointsize = 22) plot(bin, main="KOIFN",xlab="reads_per_exon_per_million_RNAseq",ylab="reads_per_exon_per_million_PolysomesSeq") dev.off() t=read.table("/data/images/proton/DKlab/orsalia/perExon/Ko1IL4.cnt.rpm.csv",header=T) png("KoIL4-PolysomeVSRNAseq.png",width=1024,height=1024,pointsize = 22) plot (log2(0.5*(t$Ko1IL4.bam_rpm+t$Ko2IL4.bam_rpm)),log2(0.5*(t$KO_IL4_replicate_I.bam_rpm+t$KO_IL4_replicate_II.bam_rpm)),pch=".",xlab="reads_per_exon_per_million_RNAseq",ylab="reads_per_exon_per_million_PolysomesSeq") dev.off() x=log2(0.5*(t$Ko1IL4.bam_rpm+t$Ko2IL4.bam_rpm));y=log2(0.5*(t$KO_IL4_replicate_I.bam_rpm+t$KO_IL4_replicate_II.bam_rpm)); x[mapply(is.infinite, x)] <- NA;y[mapply(is.infinite, y)] <- NA; bin<-hexbin(x,y, xbins=50);bin@count=log2(bin@count); png("KoIL4-PolysomeVSRNAseq-hist.png",width=1024,height=1024,pointsize = 22) plot(bin, main="KO_IL4",xlab="reads_per_exon_per_million_RNAseq",ylab="reads_per_exon_per_million_PolysomesSeq") dev.off() s=t[(t$RPKM_PHR23r.FUBP2>=1)&(t$PHR23r.FUBP2>=10),] s=s[!(duplicated(s)),] dim(s) ss=s[order(s$Geneid,-(s$fc_vs_PHR24r.IgG)),]#sort by id and max reads s=ss[!(duplicated(ss$Geneid)),] selectDE=order(s$fc_vs_PHR24r.IgG, decreasing=T)#[1:100] write.table(s[selectDE,],"RPKM_PHR23r-FUBP2_vs_PHR24r-IgG.csv",row.names=FALSE,col.names=T,quote=FALSE,sep="\t") see: /data/images/proton/DKlab/orsalia/perExon/doubleDeltaExpression1.r tWT0h=read.table("/data/images/proton/DKlab/orsalia/perExon/wt01.cnt.rpmPerGene.csv",header=T) xWT0hRNA=log2(0.5*(tWT0h$wt01.bam_rpm+tWT0h$wt02.bam_rpm));yWT0hPoly=log2(0.5*(tWT0h$WT_0h_replicate_I.bam_rpm+tWT0h$WT_0h_replicate_II_.bam_rpm)); xWT0hRNA[mapply(is.infinite, xWT0hRNA)] <- NA;yWT0hPoly[mapply(is.infinite, yWT0hPoly)] <- NA; tKO0h=read.table("/data/images/proton/DKlab/orsalia/perExon/Ko01.cnt.rpmPerGene.csv",header=T) xKO0hRNA=log2(0.5*(tKO0h$Ko01.bam_rpm+tKO0h$Ko02.bam_rpm));yKO0hPoly=log2(0.5*(tKO0h$KO_0h_replicate_I.bam_rpm+tKO0h$KO_0h_replicate_II.bam_rpm)); xKO0hRNA[mapply(is.infinite, xKO0hRNA)] <- NA;yKO0hPoly[mapply(is.infinite, yKO0hPoly)] <- NA; dd1=xWT0hRNA-yWT0hPoly; dd2=xKO0hRNA-yKO0hPoly; png("0h-wtVSko-PolysomeVSRNAseq_perGene.png",width=2048,height=2048,pointsize = 22) plot (dd1,dd2,pch=".",xlab="log2_reads_per_gene_per_million_RNA_minus_Polysomes_wt0h",ylab="delta_log2_reads_per_gene_per_million_RNA_minus_Polysomes_ko0h",xlim=c(-20,20),ylim=c(-20,20)) points(mean(na.omit(cbind(dd1))),mean(na.omit(cbind(dd2))),col="green") abline(lm(dd1~dd2), col="red") # regression line (y~x) abline(0,1, col="blue") # regression line (y~x) grid(nx=NULL, ny=NULL) set=as.logical(abs(dd1-dd2)>=5) text(dd1[set],dd2[set], labels = tWT0h$Geneid[set], pos = 3,cex=0.5) dev.off()