Merge pull request #146 from RBL-NCI/activeDev

slsevilla · web-flow · commit 1388a97eb5c0 · 2023-01-20T12:47:51.000-05:00
Active dev
diff --git a/config/cluster_config.yaml b/config/cluster_config.yaml
@@ -9,8 +9,8 @@ __default__:
     error: .%j.{wildcards}.err
 
 qc_barcode:
-    threads: 4
-    mem: 100g
+    threads: 8
+    mem: 75g
     time: 00-04:00:00
 
 demultiplex:
@@ -34,8 +34,8 @@ qc_screen_validator:
 star:
     time: 04-00:00:00
     gres: lscratch:800
-    threads: 32
-    mem: 85g
+    threads: 16
+    mem: 120g
 
 index_stats:
     threads: 8
diff --git a/config/snakemake_config.yaml b/config/snakemake_config.yaml
@@ -22,6 +22,7 @@ contrastManifest: "OUTPUT_DIR/manifests/contrasts.tsv"
 multiplexflag: "Y" #whether samples are multiplexed ["Y","N"]
 umiSeparator: "rbc:" #required for nondemultiplexed samples to determine delimiter for deduplication [":", "_", "rbc:"]
 mismatch: 1 #required for multiplexed samples, number of bp mismatches allowed in demultiplexing [1,2,3]
+barcode_qc_flag: "PROCESS" #barcodes will undergo QC to ensure uniformity within samples; if set to IGNORE this qc step will be skipped
 reference: "" #reference organism ["hg38","mm10"]
 filterlength: 20 #minimum read length to include in analysis [any int >20]
 phredQuality: 20 #minimum quality score for 3’ end trimming
diff --git a/run_snakemake.sh b/run_snakemake.sh
@@ -35,7 +35,7 @@ fi
 # functions
 #########################################################
 load_modules(){
-  if [[ $1 =~ "python" ]]; then module load python/3.8; fi
+  if [[ $1 =~ "python" ]]; then module load python/3.7; fi
   if [[ $1 =~ "snakemake" ]]; then module load snakemake/7.19.1; fi
   if [[ $1 =~ "graphviz" ]]; then module load graphviz/2.40; fi
 }
@@ -194,7 +194,7 @@ check_manifest_qc(){
           echo "The barcode manifest check FAILED. Check "${bc_file}" for more information."
           exit 1
         else
-          echo "-- barcode check completed successfully (sample $mp_id)"
+          echo "-- barcode manifest check completed successfully (sample $mp_id)"
         fi
     done
   fi
diff --git a/workflow/Snakefile b/workflow/Snakefile
@@ -41,6 +41,7 @@ contrast_manifest = config['contrastManifest']
 # user parameters
 multiplex_flag = config['multiplexflag'].capitalize()
 mismatch = config['mismatch']
+barcode_qc_flag = config['barcode_qc_flag']
 species_ref = config['reference']
 filter_length = config['filterlength']
 phredQuality = config['phredQuality']
@@ -219,6 +220,28 @@ def rename_cmd(wildcards):
         cmd = 'mv ' + source + ' ' + destination + '; ' + cmd
     return(cmd)
 
+def get_input_qc_troubleshoot(wildcards):
+    print("HERE")
+    d=dict()
+
+    l = "log_list"
+    t = "txt_bc"
+    p = "png_bc"
+
+    if (multiplex_flag == 'Y'):
+        print("m flag on")
+        d[l] = expand(join(out_dir,'log','STAR','{sp}.log'),sp=sp_list),
+        d[t] = expand(join(out_dir,'{mp}', '01_qc_post','{mp}_barcode_passed.txt'),mp=samp_dict.keys())
+        if (barcode_qc_flag == "PROCESS"):
+            print("b flag on")
+            d[p] = expand(join(out_dir,'{mp}', '01_qc_post','{mp}_barcode.png'),mp=samp_dict.keys())
+    else:
+        print("m flag off")
+        d[l] = expand(join(out_dir,'log','STAR','{sp}.log'),sp=sp_list),
+    print(d)
+    d="YO"
+    return (d)
+
 #read in contrast list for manorm
 def get_de_list():
     
@@ -381,7 +404,7 @@ rule all:
         # multiqc
         join(out_dir,'qc','multiqc_report.html'),
 
-        # qc_troubleshoot
+        # # qc_troubleshoot
         join(out_dir,'qc','qc_report_troubleshooting.html'),
 
         # dedup
@@ -408,7 +431,7 @@ rule all:
         # # qc barcode
         # expand(join(out_dir,'{mp}','01_qc_post','{mp}_barcode_counts.txt'),mp=mp_list),
         # expand(join(out_dir,'{mp}','01_qc_post','{mp}_barcode.png'),mp=mp_list),
-        # expand(join(out_dir,'{mp}','01_qc_post','{mp}_barcode.txt'),mp=mp_list),
+        # expand(join(out_dir,'{mp}','01_qc_post','{mp}_barcode_passed.txt'),mp=mp_list),
 
         # # demux OR nondemux
         # expand(join(out_dir,'{mp}', '02_preprocess','ultraplex_demux_5bc_no_match.fastq.gz'), mp=samp_dict.keys()),
@@ -462,9 +485,11 @@ else:
 if (multiplex_flag == 'Y'):
     rule qc_barcode:
         """
-        generate counts of barcodes and output to text file
-        will run python script that determines barcode expected and generates mismatches based on input
-        output barplot with top barcode counts
+        if the flag in the  config file is set to process, barcodes will be reviewed to ensure uniformtiy amongst samples.
+        - generate counts of barcodes and output to text file
+        - will run python script that determines barcode expected and generates mismatches based on input
+        - output barplot with top barcode counts
+        If it's set to IGNORE, then it will skip this processing and create the needed QC file
         """
         input:
             fq = get_fq_names,
@@ -473,30 +498,40 @@ if (multiplex_flag == 'Y'):
             memG = getmemG_80perc("qc_barcode"),
             R = join(source_dir,'workflow', 'scripts', '02_barcode_qc.R'),
             base = join(out_dir,'{mp}','01_qc_post'),
+            qc_dir=join(out_dir,'qc'),
             mm = mismatch,
             bc_len = barcode_length,
-            start_pos = 6 if barcode_length==6 else 4
+            start_pos = 6 if barcode_length==6 else 4,
+            b_qc_flag=barcode_qc_flag
         threads: getthreads("qc_barcode")
         envmodules:
             config['R'],
         output:
             counts = temp(join(out_dir,'{mp}','01_qc_post','{mp}_barcode_counts.txt')),
             png = temp(join(out_dir,'{mp}','01_qc_post','{mp}_barcode.png')),
-            txt = temp(join(out_dir,'{mp}','01_qc_post','{mp}_barcode.txt'))
+            txt = join(out_dir,'{mp}','01_qc_post','{mp}_barcode_passed.txt')
         shell:
             """
             set -exo pipefail
-            gunzip -c  {input.fq} \\
-                | awk 'NR%4==2 {{print substr($0, {params.start_pos}, {params.bc_len});}}' \\
-                | LC_ALL=C sort --buffer-size={params.memG} --parallel={threads} --temporary-directory=/lscratch/${{SLURM_JOB_ID}} -n \\
-                | uniq -c > {output.counts};
-            
-            Rscript {params.R} --sample_manifest {sample_manifest} \\
-                --multiplex_manifest {multiplex_manifest} \\
-                --barcode_input {output.counts} \\
-                --mismatch {params.mm} \\
-                --mpid {wildcards.mp} \\
-                --output_dir {params.base}
+
+            if [[ {params.b_qc_flag} == "PROCESS" ]]; then
+                gunzip -c  {input.fq} \\
+                    | awk 'NR%4==2 {{print substr($0, {params.start_pos}, {params.bc_len});}}' \\
+                    | LC_ALL=C sort --buffer-size={params.memG} --parallel={threads} --temporary-directory=/lscratch/${{SLURM_JOB_ID}} -n \\
+                    | uniq -c > {output.counts};
+                
+                Rscript {params.R} --sample_manifest {sample_manifest} \\
+                    --multiplex_manifest {multiplex_manifest} \\
+                    --barcode_input {output.counts} \\
+                    --mismatch {params.mm} \\
+                    --mpid {wildcards.mp} \\
+                    --output_dir {params.base} \\
+                    --qc_dir {params.base}
+            else
+                echo "Barcode QC checking was skipped for this run" > {output.txt}
+                touch {output.png}
+                touch {output.counts}
+            fi
             """
 
     rule demultiplex:
@@ -785,7 +820,7 @@ rule star:
         --sjdbScore {params.s_sj} \
         --winAnchorMultimapNmax {params.s_anchor}
 
-        # sort file\
+        # sort file
         samtools sort -m 80G -T $tmp_dir $tmp_dir/{params.out_prefix}Aligned.out.bam -o $tmp_dir/{params.out_prefix}Aligned.sortedByCoord.out.bam
 
         # move STAR files and final log file to output
@@ -861,66 +896,44 @@ rule multiqc:
             -o {params.out_dir}
         """
 
-# Pipeline splits to handle multiplexed QC vs non-multiplexed QC
-if (multiplex_flag == 'Y'):
-    rule qc_troubleshoot:
-        """
-        generates a PDF of barcode plots and alignment plots for qc troubleshooting
+rule qc_troubleshoot:
+    """
+    generates a PDF of barcode plots and alignment plots for qc troubleshooting
+    """
+    input:
+        log_list=expand(join(out_dir,'log','STAR','{sp}.log'),sp=sp_list),
+    params:
+        rname = "08_qc_troubleshoot",
+        R = join(source_dir,'workflow','scripts','03_qc_report.Rmd'),
+        qc = join(source_dir,'workflow','scripts','03_qc_unmapped_check.sh'),
+        log_dir = join(out_dir,'log','STAR'),
+        single_qc = single_qc_threshold,
+        proj_qc = project_qc_threshold,
+        m_flag=multiplex_flag,
+        b_flag=barcode_qc_flag,
+        txt_bc=expand(join(out_dir,'{mp}', '01_qc_post','{mp}_barcode_passed.txt'),mp=samp_dict.keys())
+    envmodules:
+        config['R']
+    output:
+        html = join(out_dir,'qc','qc_report_troubleshooting.html'),
+        qc = join(out_dir, 'log','STAR','qc_pass.txt')
+    shell:
         """
-        input:
-            log_list = expand(join(out_dir,'log','STAR','{sp}.log'),sp=sp_list),
-            txt_bc = expand(join(out_dir,'{mp}', '01_qc_post','{mp}_barcode.txt'),mp=samp_dict.keys()),
-            png_bc = expand(join(out_dir,'{mp}', '01_qc_post','{mp}_barcode.png'),mp=samp_dict.keys())
-        params:
-            rname = "08_qc_troubleshoot",
-            R = join(source_dir,'workflow','scripts','03_qc_report.Rmd'),
-            qc = join(source_dir,'workflow','scripts','03_qc_unmapped_check.sh'),
-            log_dir = join(out_dir,'log','STAR'),
-            single_qc = single_qc_threshold,
-            proj_qc = project_qc_threshold
-        envmodules:
-            config['R']
-        output:
-            html = join(out_dir,'qc','qc_report_troubleshooting.html'),
-            qc = join(out_dir, 'log','STAR','qc_pass.txt')
-        shell:
-            """
+        if [[ {params.m_flag} == "Y" ]] && [[ {params.b_flag} == "PROCESS" ]] ; then
             Rscript -e 'library(rmarkdown); \
             rmarkdown::render("{params.R}",
                 output_file = "{output.html}", \
-                params= list(log_list = "{input.log_list}", \
-                    b_txt = "{input.txt_bc}"))'
-
-            sh {params.qc} {params.log_dir} {params.single_qc} {params.proj_qc}
-            """
-else:
-    rule qc_troubleshoot:
-        """
-        generates a PDF of barcode plots and alignment plots for qc troubleshooting
-        """
-        input:
-            log_list = expand(join(out_dir,'log','STAR','{sp}.log'),sp=sp_list),
-        params:
-            rname = "08_qc_troubleshoot",
-            R = join(source_dir,'workflow','scripts','03_qc_report.Rmd'),
-            qc = join(source_dir,'workflow','scripts','03_qc_unmapped_check.sh'),
-            log_dir = join(out_dir,'log','STAR'),
-            single_qc = single_qc_threshold,
-            proj_qc = project_qc_threshold
-        envmodules:
-            config['R']
-        output:
-            html = join(out_dir,'qc','qc_report_troubleshooting.html'),
-            qc = join(out_dir, 'log','STAR','qc_pass.txt')
-        shell:
-            """
+                    params= list(log_list = "{input.log_list}", \
+                        b_txt = "{params.txt_bc}"))'
+        else
             Rscript -e 'library(rmarkdown); \
-            rmarkdown::render("{params.R}",
+                rmarkdown::render("{params.R}",
                 output_file = "{output.html}", \
-                params= list(log_list = "{input.log_list}"))'
-            
-            sh {params.qc} {params.log_dir} {params.single_qc} {params.proj_qc}
-            """
+                    params= list(log_list = "{input.log_list}"))'
+        fi
+
+        sh {params.qc} {params.log_dir} {params.single_qc} {params.proj_qc}
+        """
 
 rule dedup:
     """
diff --git a/workflow/scripts/02_barcode_qc.R b/workflow/scripts/02_barcode_qc.R
@@ -20,6 +20,7 @@ parser$add_argument("-bu","--barcode_input", dest="barcode_input", required=TRUE
 parser$add_argument("-o","--output_dir", dest="output_dir", required=TRUE, help="dir for barcode plots and summary")
 parser$add_argument("-m","--mismatch", dest="mismatch", required=TRUE, help="number of mismatches allowed")
 parser$add_argument("-i","--mpid", dest="mpid", required=TRUE, help="sampleid")
+parser$add_argument("-q","--qc_dir", dest="qc_dir", required=TRUE, help="dir for failed QC manifests")
 
 args <- parser$parse_args()
 sample_manifest = args$sample_manifest
@@ -32,12 +33,12 @@ mpid = args$mpid
 #test input
 testing="N"
 if(testing=="Y"){
-  sample_manifest  = "~/../../Volumes/CCBR_Pipeliner/iCLIP/test/sample_manifests/sample_mm10_one.tsv"
-  multiplex_manifest = "~/../../Volumes/CCBR_Pipeliner/iCLIP/test/multiplex_manifests/multiplex_mm10_one.tsv"
-  barcode_input = "~/../../Volumes/sevillas2-1/test.barcode_counts.txt"
-  output_dir = "~/../../Volumes/sevillas2-1"
+  sample_manifest  = "~/../../Volumes/RBL_NCI/Wolin/mov10_par_Y_r2_01062023/manifests/samples.tsv"
+  multiplex_manifest = "~/../../Volumes/RBL_NCI/Wolin/mov10_par_Y_r2_01062023/manifests/multiplex.tsv"
+  barcode_input = "~/../../Volumes/RBL_NCI/Wolin/mov10_par_Y_r2_01062023/sam_testing/barcode_counts.txt"
+  output_dir = "~/../../Volumes/RBL_NCI/Wolin/mov10_par_Y_r2_01062023/sam_testing"
   mismatch = 1
-  mpid="test_1"
+  mpid="CLIP"
 }
 
 #Read sample/multiplex manifests
@@ -90,9 +91,7 @@ for (rowid in rownames(barcode_combo)){
 
 #if there were errors, print message and exit
 if(error_message!=""){
-  cat(paste0("The number of differences between barocdes is less than, or equal to, the number of mismatches
-in the following barcodes:\n",error_message))
-  #exit()
+  cat(paste0("The number of differences between barocdes is less than, or equal to, the number of mismatches in the following barcodes:\n",error_message))
 } else{
   print("Barcode schema passes QC")
 }
@@ -154,7 +153,7 @@ diff_lis=setdiff(paste0("*",bc_exp),df_counts_merged$barcode_id)
 if(length(diff_lis)==0){
   print("Barcodes observed pass QC")
   #print text file
-  file_save = paste0(output_dir,'/',mpid,'_barcode.txt')
+  file_save = paste0(output_dir,'/',mpid,'_barcode_passed.txt')
   fileConn<-file(file_save)
   line1 = paste0("\n* SampleID ",mpid, " \n")
   line2 = paste0("\t + Number of mismatches allowed ", mismatch, " \n")
@@ -179,7 +178,7 @@ if(length(diff_lis)==0){
   ggsave(file_save,p)
   
 } else{
-  file_save = paste0(output_dir,'/',mpid,'_barcode_errors.txt')
+  file_save = paste0(qc_dir,'/',mpid,'_barcode_errors.txt')
   fileConn<-file(file_save)
   line1 = paste0("\n* SampleID ",mpid, " \n")
   line2 = paste0("\t + Number of mismatches allowed ", mismatch, " \n")
@@ -194,4 +193,8 @@ if(length(diff_lis)==0){
   
   writeLines(paste(line1,line2,line3,line4,line5,sep=""), fileConn)
   close(fileConn)
+
+  print(paste0("Barcodes failed QC check. Review file located here:",file_save,"."))
+  print(paste0("To check the reads go to the /project/01_preprocess/01_fastq/ dir and run `zcat {name_of_sample} | wc -l. Divide this number by 4 to get the total number of reads"))
+  print(paste0("If satisfied with this read count change the config/snakemake_config file barcode_qc_flag to IGNORE"))
 }
diff --git a/workflow/scripts/06_annotation.Rmd b/workflow/scripts/06_annotation.Rmd
@@ -2,9 +2,9 @@
 title: ""
 output: html_document
 params:
-  samplename: "Ro_Clip_u"
-  peak_in: "/Users/homanpj/OneDrive - National Institutes of Health/Loaner/Wolin/CLIP/mESC_clip_3_clip3_u/04_annotation/02_peaks/Ro_Clip_u_ALLreadPeaks_annotation_complete.txt"
-  output_table: "/Users/homanpj/OneDrive - National Institutes of Health/Loaner/Wolin/CLIP/mESC_clip_3_clip3_u/04_annotation/Ro_Clip_u_peakannotation_finalx.txt"
+  samplename: "YDKO_IgG_500"
+  peak_in: "~/../../Volumes/RBL_NCI/Wolin/mov10_par_Y_r2_01062023/04_annotation/02_peaks/YDKO_IgG_500_ALLreadPeaks_annotation_complete.txt"
+  output_table: "~/../../Volumes/RBL_NCI/Wolin/mov10_par_Y_r2_01062023/04_annotation/YDKO_IgG_500_annotation_ALLreadPeaks_final_table.txt"
   readdepth: 3
   PeakIdnt: "ALL" #MultiMap, Unique, all
 # params:
@@ -20,7 +20,7 @@ params:
 #   readdepth: 5
 #   PeakIdnt: "ALL" #MultiMap, Unique, all
 ---
- 
+
  <!-- # Load libraries, set theme -->
 ```{r envirment, include=F ,echo=F,warning=F,message=FALSE}
 rm(list=setdiff(ls(), "params"))
diff --git a/workflow/scripts/create_error_report.sh b/workflow/scripts/create_error_report.sh
