Merge pull request #36 from CCBR/v1.0.1-dev

V1.0.1 dev

aspen

@@ -198,13 +198,29 @@ sed -e "s/PIPELINE_HOME/${PIPELINE_HOME//\//\\/}/g" \
   -e "s/WORKDIR/${WORKDIR//\//\\/}/g" \
   -e "s/GENOME/${GENOME}/g" $f > $WORKDIR/$fbn
 
-for f in $MANIFEST ${PIPELINE_HOME}/resources/cluster.json ${PIPELINE_HOME}/resources/tools.yaml
+for f in ${PIPELINE_HOME}/resources/cluster.json ${PIPELINE_HOME}/resources/tools.yaml
 do
   echo "Copying essential file: $f"
   fbn=$(basename $f)
   cp $f $WORKDIR/$fbn
 done
 
+if [[ "$MANIFEST_SUPPLIED" == "true" ]];then
+  f=$MANIFEST
+  echo "Copying essential file: $f"
+  fbn=$(basename $f)
+  cp $f $WORKDIR/$fbn
+fi
+
+if [[ "$MANIFEST_SUPPLIED" == "false" ]];then
+f=$MANIFEST
+echo "Copying essential file: $f"
+fbn=$(basename $f)
+sed -e "s/PIPELINE_HOME/${PIPELINE_HOME//\//\\/}/g" \
+  -e "s/WORKDIR/${WORKDIR//\//\\/}/g" \
+  -e "s/GENOME/${GENOME}/g" $f > $WORKDIR/$fbn
+fi
+
 # copy essential folders
 for f in $ESSENTIAL_FOLDERS;do
   # rsync -az --progress ${PIPELINE_HOME}/$f $WORKDIR/
@@ -430,7 +446,7 @@ function runslurm() {
 function create_runinfo {
   modtime=$1
   if [ "$modtime" == "" ];then
-   modtime=$(stat ${WORKDIR}/runinfo.yaml|grep Modify|awk '{print $2,$3}'|awk -F"." '{print $1}'|sed "s/ //g"|sed "s/-//g"|sed "s/://g")
+   modtime=$(stat ${WORKDIR}/runinfo.yaml 2>/dev/null|grep Modify|awk '{print $2,$3}'|awk -F"." '{print $1}'|sed "s/ //g"|sed "s/-//g"|sed "s/://g")
   fi
   if [ -f ${WORKDIR}/runinfo.yaml ];then
     mv ${WORKDIR}/runinfo.yaml ${WORKDIR}/stats/runinfo.${modtime}.yaml
@@ -732,7 +748,10 @@ function main(){
   esac
   done
   WORKDIR=$(readlink -f "$WORKDIR")
+  MANIFEST_SUPPLIED="true"
+  # if manifest is empty ... aka not supplied at cli
   if [[ -z $MANIFEST ]];then
+    MANIFEST_SUPPLIED="false"
     MANIFEST=${PIPELINE_HOME}/config/samples.tsv
   fi
   echo "Working Dir       : $WORKDIR"

config/contrasts.tsv

@@ -0,0 +1 @@
+D4_Meso_iCre_Dox	iCre_D0

workflow/rules/diffatac.smk

@@ -139,7 +139,10 @@ if [ ! -d $outdir ];then mkdir $outdir;fi
 cd $outdir
 tail -n +2 {params.manifest} | cut -f1-2 | sort > {output.sampleinfo}
 
+cat {params.contrasts}
 while read g1 g2;do
+    echo "g1=$g1"
+    echo "g2=$g2"
     Rscript {params.scriptsdir}/DESeq2_runner.R \\
         --countsmatrix {input.counts} \\
         --genome {params.genome} \\

workflow/rules/init.smk

@@ -178,7 +178,11 @@ if contrastsfileexists:
             CONTRASTS = pd.DataFrame()
             print(contrastsfile + " is empty. No contrasts will be run.")
     except OSError: # contrast file is empty!
+        CONTRASTS = pd.DataFrame()
         print(contrastsfile + " is empty. No contrasts will be run.")
+else:
+    CONTRASTS = pd.DataFrame()
+    print(contrastsfile + " is absent. No contrasts will be run.")
 
 
 

workflow/scripts/DESeq2.Rmd

@@ -139,8 +139,12 @@ pander(sampleinfo,style="rmarkdown")
 dds1 <- DESeqDataSetFromMatrix(countData = as.matrix(countmatrix),
                                 colData = sampleinfo[,c("replicateName","sampleName")],
                                 design = ~ 1)
-dds1 <- DESeq(dds1)
-rld1 <- vst(dds1)
+dds1 <- DESeq(dds1) 
+if (nrow(dds1)<2000) {
+  rld1 <- varianceStabilizingTransformation(dds1, blind=TRUE)
+} else {
+  rld1 <- vst(dds1, blind=TRUE)
+}
 assayrld1 = as.data.frame(assay(rld1))
 assayrld1$row_variance = rowVars(as.matrix(assayrld1))
 assayrld1 = arrange(assayrld1,desc(row_variance))
@@ -178,7 +182,11 @@ dds2 <- DESeqDataSetFromMatrix(countData = as.matrix(countmatrix2),
                                 design = ~ sampleName)
 dds2$sampleName <- relevel(dds2$sampleName, ref=params$contrast_denominator)
 dds2 <- DESeq(dds2)
-rld2 <- vst(dds2)
+if (nrow(dds2)<2000) {
+  rld2 <- varianceStabilizingTransformation(dds2, blind=TRUE)
+} else {
+  rld2 <- vst(dds2, blind=TRUE)
+}
 pca2=prcomp(t(as.data.frame(assay(rld2))),scale.=T)
 m2.pc1 = round(pca2$sdev[1]^2/sum(pca2$sdev^2)*100,2)
 m2.pc2 = round(pca2$sdev[2]^2/sum(pca2$sdev^2)*100,2)

workflow/scripts/aggregate_results.Rmd

@@ -73,7 +73,11 @@ dds1 <- DESeqDataSetFromMatrix(countData = as.matrix(countmatrix),
                                 colData = sampleinfo[,c("replicateName","sampleName")],
                                 design = ~ 1)
 dds1 <- DESeq(dds1)
-rld1 <- vst(dds1)
+if (nrow(dds1)<2000) {
+  rld1 <- varianceStabilizingTransformation(dds1,blind=TRUE)
+} else {
+  rld1 <- vst(dds1, blind=TRUE)
+}
 assayrld1 = as.data.frame(assay(rld1))
 assayrld1$row_variance = rowVars(as.matrix(assayrld1))
 assayrld1 = arrange(assayrld1,desc(row_variance))
@@ -134,6 +138,7 @@ collist = c()
 for (i in 1:length(degfiles)){
   # i=1
   deg_file=paste(params$diffatacdir,degfiles[i],sep="/")
+  print(deg_file)
   contrast=contrasts[i]
   x <- read_deg_file(deg_file)
   all_sig_geneids <- c(all_sig_geneids,get_sig_geneids(x,params$FC_cutoff,params$FDR_cutoff))