/data/images/proton2/rnaSeq-pip-dme1.sh

cd wwwESlab
/data/images/proton/make-bam-links-tophat2.sh
source make-links.sh
awk -f /data/images/proton2/make-stranded-sampled-dme-tracks.awk  reads.txt  > make-stranded-sampled-dme-tracks.sh

@@ change/add chr to chromosomes in bam file:
http://seqanswers.com/forums/archive/index.php/t-22504.html
for i in *.bam
do
echo $i
samtools view -h $i | awk 'BEGIN{FS=OFS="\t"} (/^@/ && !/@SQ/){print $0} $2~/^SN:[1-9]|^SN:X|^SN:Y|^SN:MT|^SN:M/{print $0} $3~/^[1-9]|X|Y|MT|M/{$3="chr"$3; print $0} ' | sed 's/SN:/SN:chr/g' | sed 's/chrMT/chrM/g' | samtools view -bS - > $i.bam2
done


library(metaseqR)

the.path <- "/data/images/proton2/run430/wwwESlab"

the.contrasts.1 <- c("Control_vs_Mutant")

metaseqr(
    sample.list=file.path(the.path,"targets.txt"),
    contrast=the.contrasts.1,
    annotation="download",
    org="dm3",
    refdb="ensembl",
    count.type="utr",
    normalization="edger",
    statistics="edger",
#    normalization="deseq",
#    statistics="deseq",
    fig.format=c("png","pdf"),
    export.where=file.path(the.path,"metaseqr_ESlab_run430"),
    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
)


Normalizing with: edger
Applying gene filter length...
  Threshold below which ignored: 500
Applying gene filter avg.reads...
  Threshold below which ignored: -Inf
Applying gene filter expression...
  Threshold below which ignored: NA
Applying gene filter biotype...
  Biotypes ignored: rRNA, IG_V_pseudogene, TR_V_pseudogeneError in temp.matrix[the.dead, ] : no 'dimnames' attribute for array
In addition: Warning messages:
1: In if (whats.strand %in% c("yes", "no")) { :
  the condition has length > 1 and only the first element will be used
2: "yes" and "no" for read strandedness have been deprecated. Please use "forward", "forward" or "no". Replacing "yes" with "forward"... 
3: In DGEList(counts = gene.counts, group = classes) :
  library size of zero detected
4: In max(cumDim[cumDim <= lstats]) :
  no non-missing arguments to max; returning -Inf
5: In max(apply(avg.mat, 2, quantile, gene.filters$avg.reads$quantile)) :
  no non-missing arguments to max; returning -Inf