Merge pull request #44 from Dartmouth-Data-Analytics-Core/dev

#43 Improvements to PCA rule

Author: owenwilkins

Snakefile

@@ -13,10 +13,6 @@ sample_list = list(samples_df['sample_id'])
 
 print(config)
 
-
-
-
-
 #####~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 # define rules
 #####~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -35,7 +31,8 @@ rule all:
         expand("rsem/{sample}.genes.results", sample=sample_list) if config["run_rsem"] == "yes" else [],
         expand("rsem/{sample}.isoforms.results", sample=sample_list) if config["run_rsem"] == "yes" else [],
         "featurecounts/featurecounts.readcounts.tsv",
-        "plots/PCA_Variance_Bar_Plot.png",
+        "plots/PCA_top_PC1_vs_PC2.png",
+        "plots/PCA_top_PCA_variance_bar.png",
         "featurecounts/featurecounts.readcounts.ann.tsv",
         "featurecounts/featurecounts.readcounts_tpm.tsv",
         "featurecounts/featurecounts.readcounts_tpm.ann.tsv",
@@ -45,7 +42,7 @@ rule all:
         "featurecounts/featurecounts.readcounts_fpkm.ann.tsv",
     conda:
         "env_config/multiqc.yaml",
-    resources: cpus="10", maxtime="2:00:00", mem_mb="60gb",
+    resources: cpus="10", maxtime="2:00:00", mem_mb=60000,
 
     params:
         layout=config["layout"],
@@ -96,7 +93,7 @@ rule trimming:
         nextseq_flag = config["cutadapt_nextseq_flag"]
     conda:
         "env_config/cutadapt.yaml",
-    resources: cpus="10", maxtime="2:00:00", mem_mb="60gb",
+    resources: cpus="10", maxtime="2:00:00", mem_mb=60000,
 
     shell: """
         if  [ "{params.layout}" == "paired" ] 
@@ -141,7 +138,7 @@ if config["aligner_name"]=="star":
       conda:
           "env_config/alignment.yaml",
 
-      resources: cpus="10", maxtime="8:00:00", mem_mb="120gb",
+      resources: cpus="10", maxtime="8:00:00", mem_mb=120000,
 
       shell: """
         align_folder="sample_ref/STAR_index"
@@ -183,7 +180,7 @@ if config["aligner_name"]=="star":
       conda:
           "env_config/alignment.yaml",
 
-      resources: cpus="5", maxtime="8:00:00", mem_mb="100gb",
+      resources: cpus="5", maxtime="8:00:00", mem_mb=100000,
 
       shell: """
         align_folder=`cat alignment/index_status.txt`
@@ -238,7 +235,7 @@ if config["aligner_name"]=="hisat":
       conda:
           "env_config/alignment.yaml",
 
-      resources: cpus="4", maxtime="8:00:00", mem_mb="40gb",
+      resources: cpus="4", maxtime="8:00:00", mem_mb=40000,
 
       shell: """
           {params.hisat2} \
@@ -272,7 +269,7 @@ rule alignment_metrics:
     conda:
         "env_config/samtools.yaml",
 
-    resources: cpus="2", maxtime="8:00:00", mem_mb="20gb",
+    resources: cpus="2", maxtime="8:00:00", mem_mb=20000,
 
     shell: """
             {params.samtools} flagstat alignment/{params.sample}.srt.bam > alignment/stats/{params.sample}.srt.bam.flagstat
@@ -290,7 +287,7 @@ rule picard_markdup:
     conda:
         "env_config/picard.yaml",
 
-    resources: cpus="2", maxtime="30:00", mem_mb="20gb",
+    resources: cpus="2", maxtime="30:00", mem_mb=20000,
 
     shell: """
             {params.picard} -Xmx2G -Xms2G  \
@@ -319,7 +316,7 @@ rule picard_collectmetrics:
     conda:
         "env_config/picard.yaml",
 
-    resources: cpus="2", maxtime="8:00:00", mem_mb="20gb",
+    resources: cpus="2", maxtime="8:00:00", mem_mb=20000,
 
     shell: """
         {params.picard} -Xmx2G -Xms2G \
@@ -344,7 +341,7 @@ rule rsem:
         rsem_paired_flag = '--paired-end' if config["layout"]=='paired' else '',
     conda:
         "env_config/rsem.yaml",
-    resources: cpus="10", maxtime="8:00:00", mem_mb="60gb",
+    resources: cpus="10", maxtime="8:00:00", mem_mb=60000,
 
     shell: """   
         {params.rsem_calc_exp_path} \
@@ -380,7 +377,7 @@ rule featurecounts:
     conda:
         "env_config/featurecounts.yaml",
 
-    resources: cpus="10", maxtime="8:00:00", mem_mb="100gb",
+    resources: cpus="10", maxtime="8:00:00", mem_mb=100000,
 
     shell: """
         {params.featurecounts} -T 32 {params.pair_flag} -s {params.strand}  -a {params.gtf} -o featurecounts/featurecounts.readcounts.raw.tsv {input}
@@ -400,35 +397,21 @@ rule featurecounts:
         fi
     """
 
-# The number of genes compared for PCA, chosen by largest variance
-num_genes_compared = 500
-
 rule pca_plots:
     input: "featurecounts/featurecounts.readcounts.tsv",
 
     output:
-        "plots/Heatmap_scaled_"+str(num_genes_compared)+"_features.png",
-        # there potentially could be more, but this plot must exist. Make sure -p flag has number at least 2 if specified
-        "plots/PCA_1_vs_2.png",
-        "plots/PCA_Variance_Bar_Plot.png",
-        "plots/Gene_Variance_Plot.png",
-
+        "plots/PCA_top_PC1_vs_PC2.png",
+        "plots/PCA_top_PCA_variance_bar.png",
     params:
-        num_genes = num_genes_compared,
-        pca_plot_script = config['pca_plot_script'],
-        
+        pca_plot_script = config['pca_plot_script'],   
     conda:
-        # uses a subset of the packages that featurecounts does
         "env_config/pcaplot.yaml",
-
-    resources: cpus="1", maxtime="1:00:00", mem_mb="2gb",
-
+    resources: cpus="1", maxtime="1:00:00", mem_mb=2000,
     shell: """
         python {params.pca_plot_script} \
         featurecounts/featurecounts.readcounts.tsv \
-        plots \
-        --genes_considered {params.num_genes} \
-        --color_file sample_ref/sample_colors_hex.tsv
+        plots
     """
 
 
@@ -447,7 +430,7 @@ rule check_refs:
         picard_rrna_list = config["picard_rrna_list"],
         run_rsem = config["run_rsem"],
         rsem_ref = config["rsem_ref_path"],
-    resources: cpus="1", maxtime="1:00:00", mem_mb="2gb",
+    resources: cpus="1", maxtime="1:00:00", mem_mb=2000,
     shell: """   
         
         echo "\nChecking for reference annotation GTF file..."
@@ -534,7 +517,7 @@ rule build_refs:
         rsem_prepare_path = config["rsem_prep_ref_path"],
     conda:
           "env_config/build_refs.yaml",
-    resources: cpus="12", maxtime="8:00:00", mem_mb="48gb",
+    resources: cpus="12", maxtime="8:00:00", mem_mb=48000,
     shell: """
             REF_NAME=`basename {params.ref_fa} .fa`
             mkdir -p ref/pipeline_refs

cluster_profile/config.yaml

@@ -1,6 +1,6 @@
 jobs: 10
 jobname: "{rule}.{jobid}"
-cluster: "sbatch -t {resources.maxtime} --mem={resources.mem_mb} -c {resources.cpus} --output=log_{rule}_%j.out --mail-type=FAIL"
+cluster: "sbatch -t {resources.maxtime} --mem={resources.mem_mb} -c {resources.cpus} --output=slurm_logs/log_{rule}_%j.out --mail-type=FAIL"
 default-resources: [cpus=1, time_min=60]
 max-jobs-per-second: 5
 max-status-checks-per-second: 20

env_config/multiqc.yaml

@@ -6,6 +6,10 @@ channels:
 dependencies:
   - nomkl
   - python=3.10
-  - multiqc==1.12
+  - multiqc=1.24.1
+  - setuptools
+  - pip
+  - pip:
+    - kaleido<1
 variables:
   PYTHONIOENCODING: utf-8

env_config/pcaplot.yaml

@@ -1,13 +1,13 @@
-name: pcaplot
+name: pca_plotnine_env
 channels:
-  - conda-forge
-  - bioconda
   - defaults
 dependencies:
-  - nomkl
-  - numpy==1.23.1
-  - pandas==1.4.3
-  - seaborn==0.11.2
-  - scikit-learn==1.1.1
-  - matplotlib<3.7  # Pin to older matplotlib
+  - python=3.11
+  - scikit-learn
+  - matplotlib
+  - numpy
+  - seaborn
+  - plotnine
+  - pandas
+  - fastcluster
 

job.script.sh

@@ -25,6 +25,9 @@
 source /optnfs/common/miniconda3/etc/profile.d/conda.sh
 conda activate /dartfs/rc/nosnapshots/G/GMBSR_refs/envs/snakemake
 
+#----- Make slurm logs folder
+mkdir -p slurm_logs
+
 #----- Call Snakemake
 snakemake -s Snakefile  \
 	--conda-frontend conda \