Merge pull request #128 from griffithlab/improved_compare_junctions

update scripts

Author: kcotto

docs/workflow.md

@@ -0,0 +1,86 @@
+# RegTools example workflow
+
+This is an example workflow for running RegTools on a cohort of samples. This analysis requires that there be a vcf and RNA bam/cram file for each samples. The outline described below was used to run our own analysis on TCGA data.
+
+By the end of the analysis, the directory structure should look like so:
+
+- Project/ (SCLC/)
+  - all_splicing_variants*.bed
+  - paths.tsv
+  - make_vcfs.sh
+  - dir_names.tsv
+  - variants_all_sorted.vcf.gz
+  - variants_all_sorted.vcf.gz.tbi
+  - samples/
+    - Sample_1/
+      - tumor_rna_alignments.bam
+      - tumor_rna_alignments.bam.bai
+      - variants.per_gene.vep.vcf.gz
+      - variants.per_gene.vep.vcf.gz.tbi
+      - variants.ensembl
+      - logs/
+      - output/
+        - cse_identify_filtered_*
+        - cse_identify_filtered_compare_*
+        - variants*.bed
+  - compare_junctions/
+    - hist/
+      - junction_pvalues_*.tsv
+
+## Set tag and parameter shell arguments
+
+tag=<tag> param=<run option> (e.g. tag=default param=""; tag=E param="-E"; tag=i20e5 param="-i 20 -e 5")
+
+# run regtools cse identify with desired options for selecting variant and window size
+for i in samples/*/; do bsub -oo $i/logs/regtools_actual_$tag.lsf regtools cis-splice-effects identify $param -o ${i}/output/cse_identify_filtered_$tag.tsv -j ${i}/output/cse_identify_filtered_$tag.bed -v ${i}/output/cse_identify_filtered_$tag.vcf ${i}/variants.per_gene.vep.vcf.gz ${i}/tumor_rna_alignments.bam /gscmnt/gc2602/griffithlab/regtools/yafeng/GRCh37.fa /gscmnt/gc2602/griffithlab/regtools/GRCh37.87.exons.sorted.gtf; done
+
+# make variant.bed
+for i in samples/*/; do bsub -oo $i/logs/make_variant_bed_$tag.lsf bash /gscmnt/gc2602/griffithlab/regtools/yafeng/scripts/variants.sh ${i}/output/cse_identify_filtered_$tag.tsv ${i}/output/variants_$tag.bed; done 
+
+# make bed (really just tsv with columns: chrom start end samples) with all variants that were deemed significant to splicing across all samples
+echo -e 'chrom\tstart\tend\tsamples' > all_splicing_variants_$tag.bed
+for i in samples/*/; do j=${i##samples/}; uniq ${i}output/variants_$tag.bed | awk -v var=${j%%/} '{print $0 "\t" var}' >> all_splicing_variants_$tag.bed; done
+
+# make vcf of all variants across all samples (from variants.vcf)
+vcf-concat samples/*/variants.per_gene.vep.vcf.gz | vcf-sort > all_variants_sorted.vcf
+bgzip all_variants_sorted.vcf 
+tabix all_variants_sorted.vcf.gz 
+
+# run cis-splice effects identify on all samples with all variants (with $tag options as example)
+for i in samples/*/; do bsub -oo $i/logs/regtools_compare_$tag.lsf regtools cis-splice-effects identify $param -o ${i}/output/cse_identify_filtered_compare_$tag.tsv -j ${i}/output/cse_identify_filtered_compare_$tag.bed -v ${i}/output/cse_identify_filtered_compare_$tag.vcf variants_all_sorted.vcf.gz ${i}/tumor_rna_alignments.bam /gscmnt/gc2602/griffithlab/regtools/yafeng/GRCh37.fa /gscmnt/gc2602/griffithlab/regtools/GRCh37.87.exons.sorted.gtf; done
+
+<===== JUNCTION COMPARISON ANALYSIS =====>
+
+# make directories
+mkdir compare_junctions/
+mkdir compare_junctions/outlier
+mkdir compare_junctions/ratio
+
+# outlier analysis
+bsub -M 650000000 -R 'select[mem>65000] span[hosts=1] rusage[mem=65000]' /gscuser/zskidmor/R-3.3.0/bin/Rscript --vanilla /gscmnt/gc2602/griffithlab/regtools/yafeng/scripts/compare_junctions_outlier.R $tag
+
+# ratio analysis
+bsub -M 650000000 -R 'select[mem>65000] span[hosts=1] rusage[mem=65000]' /gscuser/zskidmor/R-3.3.0/bin/Rscript --vanilla /gscmnt/gc2602/griffithlab/regtools/yafeng/scripts/compare_junctions_ratio.R $tag
+
+# vep comparison (outlier)
+bsub /gscuser/zskidmor/R-3.3.0/bin/Rscript --vanilla /gscmnt/gc2602/griffithlab/regtools/yafeng/scripts/vep_compare.R outlier $tag
+
+# vep comparison (ratio)
+bsub /gscuser/zskidmor/R-3.3.0/bin/Rscript --vanilla /gscmnt/gc2602/griffithlab/regtools/yafeng/scripts/vep_compare.R ratio $tag
+
+NOTE: in the vep comparison, since regtools doesn't really have a concept of "variants" per se but rather "variant positions" (doesn't care about the actual alleles), we are really counting (positions of variants found splicing-significant by vep AND found significant by regtools) / (positions of variants found significant by found significant by regtools)
+
+To count:
+
+echo -e "anchor\tvep_significant_vars\ttotal_significant_vars\tpercent_vep_significant" > comparison_counts.tsv
+for i in default i20e5 i50e5 E I; do 
+	echo -n $i >> comparison_counts.tsv; 
+	echo -en "\t" >> comparison_counts.tsv; 
+	vep=$(cut -f 1,2,7 vep_comparison_$i.tsv | sort | uniq | grep "TRUE" | wc -l)
+	echo -n $vep >> comparison_counts.tsv; 
+	echo -en "\t" >> comparison_counts.tsv;
+	total=$(($(cut -f 1,2,7 vep_comparison_$i.tsv | sort | uniq | wc -l)-1)) 
+	echo -n $total >> comparison_counts.tsv;
+	echo -en "\t" >> comparison_counts.tsv;
+	echo $(bc -l <<< "$vep/$total") >> comparison_counts.tsv;
+done

scripts/compare_junctions_hist_v2.R

@@ -5,6 +5,7 @@
 # load libraries
 library(data.table)
 library(plyr)
+library(tidyverse)
 
 debug = F
 
@@ -20,12 +21,11 @@ if (debug){
     stop("Please specify an option tag (e.g. \"default\", \"i20e5\")")
   }
 }
- 
+
 # tag = 'I'
 # input_file = '~/Desktop/CHOL/all_splicing_variants_I.bed'
 
 # All splicing relevant variants (union of rows from variants.bed files; add column with comma-separated list of sample names)
-input_file = args[2]
 all_splicing_variants = unique(data.table::fread(input_file), sep = '\t', header = T, stringsAsFactors = FALSE)
 colnames(all_splicing_variants) <- c("chrom", "start", "end", "samples")
 
@@ -158,7 +158,6 @@ a <- function(x, y, z){
   return(x)
 }
 x <- mapply(a, regtools_data$norm_scores_non, length(all_samples), regtools_data$samples)
-x = split(x, rep(1:ncol(x), each = nrow(x)))
 regtools_data$norm_scores_non = x
 print("test7")
 
@@ -214,6 +213,6 @@ all_splicing_variants <- as.data.table(all_splicing_variants)
 regtools_data = regtools_data %>% distinct()
 
 
-write.table(regtools_data, file=paste(input_file, "_out.tsv", sep=""), quote=FALSE, sep='\t', row.names = F)
+write.table(regtools_data, file=paste(input_file, "_out_test.tsv", sep=""), quote=FALSE, sep='\t', row.names = F)
 
 })

scripts/stats_wrapper.py

@@ -17,9 +17,11 @@
 tag = args.tag
 cwd = os.getcwd()
 
-lines_per_file = 50000
+lines_per_file = 1000
 smallfile = None
+num_small_file = 0
 with open(f'all_splicing_variants_{tag}.bed', 'r') as bigfile:
+    num_small_file +=1
     for lineno, line in enumerate(bigfile):
         if lineno % lines_per_file == 0:
             if smallfile:
@@ -28,11 +30,13 @@
             smallfile = open(small_filename, "w")
         smallfile.write(line)
     if smallfile:
+        num_small_file += 1
         smallfile.close()
 #get chunks
 files = glob.glob('small_file_*')
+files.sort()
 for file in files:
-    subprocess.run(f'Rscript --vanilla /home/ec2-user/workspace/regtools/scripts/compare_junctions_hist_v2.R {tag} ~{cwd}/{file}', shell=True, check=True)
+    subprocess.run(f'Rscript --vanilla ~/Git/regtools/scripts/compare_junctions_hist_v2.R {tag} {file}', shell=True, check=False)
 output_files = glob.glob("*_out.tsv")
 output_files.sort()  # glob lacks reliable ordering, so impose your own if output order matters
 with open(f'junction_pvalues_{tag}.tsv', 'wb') as outfile:
@@ -43,5 +47,6 @@
             # Block copy rest of file from input to output without parsing
             shutil.copyfileobj(infile, outfile)
             print(fname + " has been imported.")
-os.remove('small_file*')
+for file in files:
+    os.remove(file)