reczko@max:/data/images/proton/run296$ source ../bam2fastq.sh 
mv IonXpress_001*fq R_001_GK3R1-T1.fq
mv IonXpress_002*fq R_002_GK3R2-T2.fq
mv IonXpress_003*fq R_003_GK3R3-T3.fq
mv IonXpress_004*fq R_004_GK3R4-U1.fq
mv IonXpress_006*fq R_006_GK3R6-U3.fq
mv IonXpress_007*fq R_007_GK3R7-G1.fq
mv IonXpress_008*fq R_008_GK3R8-G2.fq
mv IonXpress_009*fq R_009_GK3R9-G3.fq
mv IonXpress_010*fq R_010_GK3R5-U2.fq


reczko@max:/data/images/proton/run296$ ../rnaSeq-pip6fast-mm10.sh
reczko@max:/data/images/proton/run296$ ./pip.sh &> pip.log &
#unsampled:
reczko@max:/data/images/proton/run296/www$ source ../../make-stranded-mm10-tracks.sh
#normalized
reczko@max:/data/images/proton/run296/www$ awk -f ../../make-stranded-sampled-mm10-tracks1.awk reads.txt > make-stranded-sampled-mm10-tracks1.sh
reczko@max:/data/images/proton/run296/www$ source make-stranded-sampled-mm10-tracks1.sh
GK3R1-T1.bam ...


reczko@max:/data/images/proton/run296/www$ 
./cufflinks.sh &> cufflinks.log3
./cuffdiff.sh  &> cuffdiff.log &
./cuffdiff_G_vs_U.sh &> cuffdiff_G_vs_U.log &


library(cummeRbund)
cuff<-readCufflinks()
genes.MDS.rep<-MDSplot(genes(cuff),replicates=T)
genes.MDS.rep
png("T_vs_U_genes_MDSplot.png")
genes.MDS.rep
dev.off()
genes.MDS.rep<-MDSplot(isoforms(cuff),replicates=T)
genes.MDS.rep
png("T_vs_U_isoforms_MDSplot.png")
genes.MDS.rep
dev.off()
sig2 <- getSig(cuff, alpha=0.002, level='genes')
length(sig2)
[1] 1242
sigGenes2 <- getGenes(cuff, sig2)
h<-csHeatmap(sigGenes2,cluster='both')
png("T_vs_U_genes_alpha0.002.png",width=1024,height=20000,pointsize = 10)
h
dev.off()

known=(!is.na(sigGenes2@annotation$gene_short_name))
k=sig2[known]
sigGenes2 <- getGenes(cuff, k)
h<-csHeatmap(sigGenes2,cluster='both')

png("T_vs_U_known_genes_alpha0.002.png",width=400,height=2400,pointsize = 18)
h
dev.off()
h<-csHeatmap(sigGenes2,cluster='both',replicates=T)
png("T_vs_U_known_genes_alpha0.002_with_replicates.png",width=900,height=2400,pointsize = 18)
h
dev.off()

v<-csVolcano(genes(cuff),"Tcondition","U", showSignificant=TRUE, alpha=0.05,xlimits=c(-13,13))
v
png("T_vs_U_genes_alpha0.05_volcano.png")
v
dev.off()

#@
reczko@max:/data/images/proton/run296/www$ ./cuffdiff_G_vs_U.sh &> cuffdiff_G_vs_U.log &

setwd("/data/images/proton/run296/www/cuffdiff_G_vs_U")
cuff2<-readCufflinks()
genes.MDS.rep<-MDSplot(isoforms(cuff2),replicates=T)
png("G_vs_U_isoforms_MDSplot.png",width=1024,height=900)
genes.MDS.rep
dev.off()

sig2 <- getSig(cuff2, alpha=0.05, level='genes')
sigGenes2 <- getGenes(cuff2, sig2)
known=(!is.na(sigGenes2@annotation$gene_short_name))
k=sig2[known]
sigGenes2 <- getGenes(cuff, k)
h<-csHeatmap(sigGenes2,cluster='both',replicates=T)
png("G_vs_U_known_genes_alpha0.05_with_replicates.png")
h
dev.off()


v<-csVolcano(genes(cuff2),"Gcondition","U", showSignificant=TRUE, alpha=0.05,xlimits=c(-13,13))
v
png("G_vs_U_genes_alpha0.05_volcano.png")
v
dev.off()

#@ add biotypes
reczko@max:/data/images/proton/run296/www/merged_asm$ awk -f /data/results/tools/rnaseq/get-xloc-tcons-oID-biotype-from-merged1.awk merged.gtf |uniq > merged.gtf.bt
Biotypes read from /data/results/reference/mmu/mm10/Mus_musculus.GRCm38.84.bt
reczko@max:/data/images/proton/run296/www/merged_asm$ cd ../cuffdiff_T_vs_U/
reczko@max:/data/images/proton/run296/www/cuffdiff_T_vs_U$ awk -f /data/results/tools/rnaseq/get-biotypes-for-isoforms1.awk isoform_exp.diff > isoform_exp.diff.csv
reczko@max:/data/images/proton/run296/www/cuffdiff_T_vs_U$ cd ../cuffdiff_G_vs_U/
reczko@max:/data/images/proton/run296/www/cuffdiff_G_vs_U$ awk -f /data/results/tools/rnaseq/get-biotypes-for-isoforms1.awk isoform_exp.diff > isoform_exp.diff.csv


#@ metaseqR analysis
/home/moulos/Rbase/bin/R

library(metaseqR)
the.path="/data/images/proton/run296/www"
metaseqr(
    sample.list=file.path(the.path,"targets1.txt"),
    contrast=c("T_vs_U"),
    org="mm10",
    refdb="refseq",
    normalization="deseq",
    statistics="deseq",
    export.where=file.path(the.path,"metaseqr_T_vs_U"),
    qc.plots=c(
        "mds","biodetection","countsbio","saturation","readnoise",
        "correl","pairwise","boxplot","volcano","biodist","deheatmap"
    ),
    pcut=0.05,gene.filters=list(length=list(length=500),avg.reads=list(average.per.bp=100,quantile=0.25),expression=list(median=FALSE,mean=FALSE,quantile=0.25,known=NA,custom=NA),biotype=NULL),
    exon.filters=NULL,
    export.what=c("flags","annotation","p.value","adj.p.value","fold.change","stats",
        "counts"),
    export.scale=c("log2"),
    export.values="normalized",
    export.counts.table=TRUE,
    restrict.cores=0.3,
fig.format=c("png","pdf")
)

metaseqr(
    sample.list=file.path(the.path,"targets1.txt"),
    contrast=c("G_vs_U"),
    org="mm10",
    refdb="refseq",
    normalization="deseq",
    statistics="deseq",
    export.where=file.path(the.path,"metaseqr_G_vs_U"),
    qc.plots=c(
        "mds","biodetection","countsbio","saturation","readnoise",
        "correl","pairwise","boxplot","volcano","biodist","deheatmap"
    ),
    pcut=0.05,gene.filters=list(length=list(length=500),avg.reads=list(average.per.bp=100,quantile=0.25),expression=list(median=FALSE,mean=FALSE,quantile=0.25,known=NA,custom=NA),biotype=NULL),
    exon.filters=NULL,
    export.what=c("flags","annotation","p.value","adj.p.value","fold.change","stats",
        "counts"),
    export.scale=c("log2"),
    export.values="normalized",
    export.counts.table=TRUE,
    restrict.cores=0.3,
fig.format=c("png","pdf")
)