cd /data/images/proton2/run415/www-GKlab
source make-links.sh 
ln -s /data/images/proton2/run416/tophat2_001/sort_uniq.bam GK3R227_J-wt1.bam
ln -s /data/images/proton2/run416/tophat2_005/sort_uniq.bam GK3R231_J-Dare2.bam
ln -s /data/images/proton2/run416/tophat2_006/sort_uniq.bam GK3R232_J-Dare3.bam
ln -s /data/images/proton2/run416/tophat2_007/sort_uniq.bam GK3R233_J-Dare4.bam
ln -s /data/images/proton2/run416/tophat2_008/sort_uniq.bam GK3R234_J-Daredd1.bam
ln -s /data/images/proton2/run416/tophat2_011/sort_uniq.bam GK3R237_J-Daredd4.bam
ln -s /data/images/proton2/run416/tophat2_013/sort_uniq.bam GK3R239_J-DareTie2.bam
ln -s /data/images/proton2/run416/tophat2_014/sort_uniq.bam GK3R240_J-DareTie3.bam

awk -f /data/images/proton/make-metaseqr-targets1.awk ../samples.txt > targets.txt
awk -f /data/images/proton/make-metaseqr-targets1.awk ../../run416/samples.txt > targets.txt2
#merge targets.txt and targets.txt2

the.contrasts.1 <- c(
"wt_vs_Dare",
"wt_vs_Daredd",
"wt_vs_DareTie",
"Dare_vs_Daredd",
"Dare_vs_DareTie",
"Daredd_vs_DareTie")
the.path="/data/images/proton2/run415/www-GKlab"

metaseqr(
    sample.list=file.path(the.path,"targets.txt"),
    contrast=the.contrasts.1,
    annotation="download",
    org="mm10",
    count.type="utr",
    normalization="deseq",
    statistics="deseq",
    fig.format=c("png","pdf"),
    export.where=file.path(the.path,"metaseqr_quantseq_MSakkou_run415_416"),
    restrict.cores=0.5,
    qc.plots=c(
        "mds","biodetection","countsbio","saturation","readnoise","filtered",
        "correl","pairwise","boxplot","gcbias","lengthbias","meandiff",
        "meanvar","biodist","volcano","deheatmap"
    ),
    exon.filters=NULL,
    gene.filters=list(
       length=list(
            length=500
        ),
        avg.reads=list(
            average.per.bp=100,
            quantile=0.25
        ),
        expression=list(
            median=TRUE,
            mean=FALSE,
            quantile=NA,
            known=NA,
            custom=NA
        ),
        biotype=get.defaults("biotype.filter","mm10")
    ),
    pcut=0.05,
    export.what=c("annotation","p.value","adj.p.value","fold.change",
		"counts","flags"),
    export.scale=c("log2","rpgm"),
    export.values="normalized",
    export.counts.table=TRUE,
	report.top=0.05
)