From 6a10c4750cc0df94bd59affb70daa7bcf0a9abb6 Mon Sep 17 00:00:00 2001 From: Kelsy Cotto Date: Wed, 8 Jan 2020 17:13:24 -0600 Subject: [PATCH 1/3] cleanup docs --- docs/workflow.md | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/docs/workflow.md b/docs/workflow.md index abb23c7..fc0126c 100644 --- a/docs/workflow.md +++ b/docs/workflow.md @@ -1,6 +1,6 @@ # 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. +This is an example workflow for running RegTools on a cohort of samples. This analysis requires that there be a VCF and RNA bam file for each sample. The workflow described below was used to run our own analysis on TCGA data. By the end of the analysis, the directory structure should look like the example below. The `*` in the example below refers to the tag/parameter used to run `regtools cis-splice-effects identify` with. @@ -40,7 +40,9 @@ By the end of the analysis, the directory structure should look like the example - junction_pvalues_*.tsv ``` -### Set tag and parameter shell arguments +## RegTools commands + +**Set tag and parameter shell arguments** ```bash tag= @@ -48,57 +50,56 @@ param= # (e.g. tag=default param=""; tag=E param="-E"; tag=i20e5 param="-i 20 -e 5") ``` -### Run `regtools cis-splice-effects identify` with desired options for selecting variant and window size +**Run `regtools cis-splice-effects identify` with desired options for selecting variant and window size** ```bash for i in samples/*/; do 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 /reference.fa reference.gtf; done ``` -### Make `variant.bed` for each sample +**Make `variant.bed` for each sample** ```bash for i in samples/*/; do bash variants.sh ${i}/output/cse_identify_filtered_$tag.tsv ${i}/output/variants_$tag.bed; done ``` -### Combine each sample's `variant.bed` file per tag to get all variants that were deemed significant to splicing across all samples for a given tag +**Combine each sample's `variant.bed` file per tag to get all variants that were deemed significant to splicing across all samples for a given tag** ```bash 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 each sample's variants.vcf). Then, compress it and index it +**Make vcf of all variants across all samples (from each sample's variants.vcf). Then, compress it and index it.** ```bash vcf-concat samples/*/variants.vcf.gz | vcf-sort > all_variants_sorted.vcf -###### optional ###### bgzip all_variants_sorted.vcf tabix all_variants_sorted.vcf.gz ``` -### Run `regtools cis-splice effects identify` on all samples with all variants (with `$tag` options as example) +**Run `regtools cis-splice effects identify` on all samples with all variants (with `$tag` options as example)** -```bash +``` 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 all_variants_sorted.vcf.gz ${i}/tumor_rna_alignments.bam reference.fa reference.gtf; done ``` -## Beginning of compare junctions analysis +## Statistical analysis -### Make directory to store comparison results +**Make directory to store comparison results** ```bash mkdir -p compare_junctions/hist ``` -### Run `compare_junctions_hist.R` on sample data +**Run `compare_junctions_hist.R` on sample data** ```bash Rscript --vanilla compare_junctions_hist.R ``` -### Run `filter_and_BH.R` to adjust p values and filter results +**Run `filter_and_BH.R` to adjust p values and filter results** ```bash Rscript --vanilla filter_and_BH.R From 36170b108d7d0b214091dd92b1e37cd61e998c3f Mon Sep 17 00:00:00 2001 From: Yang-Yang Feng Date: Fri, 10 Jan 2020 13:56:37 -0600 Subject: [PATCH 2/3] added min pvalue --- scripts/compare_junctions_hist_v2.R | 56 +++++++++++++++++++++++------ 1 file changed, 45 insertions(+), 11 deletions(-) diff --git a/scripts/compare_junctions_hist_v2.R b/scripts/compare_junctions_hist_v2.R index 4a2e27b..176bedd 100644 --- a/scripts/compare_junctions_hist_v2.R +++ b/scripts/compare_junctions_hist_v2.R @@ -22,8 +22,8 @@ debug = F # } # } -tag = 'E' -input_file = 'all_splicing_variants_E.bed' +tag = 'default' +input_file = 'all_splicing_variants_default.bed' # All splicing relevant variants (union of rows from variants.bed files; add column with comma-separated list of sample names) all_splicing_variants = unique(data.table::fread(input_file), sep = '\t', header = T, stringsAsFactors = FALSE) @@ -187,14 +187,14 @@ get_mean <- function(x){ return(x) } -x <- mapply(get_mean, regtools_data$norm_scores_non) -regtools_data$mean_norm_score_non <- x - get_sd <- function(x){ x <- sd(as.numeric(x)) return(x) } +x <- mapply(get_mean, regtools_data$norm_scores_non) +regtools_data$mean_norm_score_non <- x + x <- mapply(get_sd, regtools_data$norm_scores_non) regtools_data$sd_norm_score_non <- x @@ -222,12 +222,21 @@ regtools_data$mean_norm_score_variant <- x x <- mapply(get_sd, regtools_data$norm_scores_variant) regtools_data$sd_norm_score_variant <- x +get_min <- function(x){ + x <- min(as.numeric(x)) + return(x) +} + +x <- mapply(get_min, regtools_data$norm_scores_variant) +regtools_data$min_norm_score_variant <- x + print("test7") ################ calculate p-values ############################################ +# calculate using mean a <- function(x){ - variant_norm_score = as.numeric(unlist(strsplit(x[['norm_scores_variant']], ',', fixed=TRUE))) + variant_norm_score = mean(as.numeric(unlist(strsplit(x[['norm_scores_variant']], ',', fixed=TRUE)))) if(length(x[['norm_scores_non']]) <= 1){ return(0) } @@ -235,8 +244,32 @@ a <- function(x){ all_norm_scores = c(x$norm_scores_non, variant_norm_score) countable = rank(all_norm_scores) num_samples = str_count(x$norm_scores_variant, ',') + 1 - non_variant_norm_scores_ranked = head(countable, (-1 * num_samples)) - variant_norm_score_ranked = tail(countable, num_samples) + non_variant_norm_scores_ranked = head(countable, -1) + variant_norm_score_ranked = tail(countable, 1) + histinfo = hist(non_variant_norm_scores_ranked, + breaks = seq(0.5, max(non_variant_norm_scores_ranked)+1.5, by=1), plot=F) + mids = histinfo$mids + cd = cumsum(histinfo$density) + underestimate = max(which(mids <= variant_norm_score_ranked)) + pvalue = 1-cd[underestimate] + return(pvalue) +} + +# calculate using min +b <- function(x){ + variant_norm_score = min(as.numeric(unlist(strsplit(x[['norm_scores_variant']], ',', fixed=TRUE)))) + if(variant_norm_score == 0){ + return(1) + } + if(length(x[['norm_scores_non']]) <= 1){ + return(0) + } + + all_norm_scores = c(x$norm_scores_non, variant_norm_score) + countable = rank(all_norm_scores) + num_samples = str_count(x$norm_scores_variant, ',') + 1 + non_variant_norm_scores_ranked = head(countable, -1) + variant_norm_score_ranked = tail(countable, 1) histinfo = hist(non_variant_norm_scores_ranked, breaks = seq(0.5, max(non_variant_norm_scores_ranked)+1.5, by=1), plot=F) mids = histinfo$mids @@ -246,7 +279,8 @@ a <- function(x){ return(pvalue) } -regtools_data$p_value <- apply(regtools_data, 1, a) +regtools_data$p_value_mean <- apply(regtools_data, 1, a) +regtools_data$p_value_min <- apply(regtools_data, 1, b) print("Number of rows in data.table") print(length(regtools_data$samples)) @@ -258,12 +292,12 @@ regtools_data$norm_scores_non <- unlist(lapply(regtools_data$norm_scores_non,pas columns_to_keep = c('samples', 'variant_info.x', 'genes', 'sample', "chrom.x", "start.x", "end.x", 'strand.x', 'anchor.x', 'info', 'names', 'mean_norm_score_variant', 'sd_norm_score_variant', 'norm_scores_variant', 'total_score_variant', 'mean_norm_score_non', 'sd_norm_score_non', 'norm_scores_non', - 'total_score_non', 'p_value') + 'total_score_non', 'p_value_mean','p_value_min') regtools_data = subset(regtools_data, select=columns_to_keep) colnames(regtools_data) <- c('variant_samples', 'variant_info', 'genes', 'junction_samples', "chrom", "start", "end", 'strand', 'anchor', 'variant_junction_info', 'names', 'mean_norm_score_variant', 'sd_norm_score_variant', 'norm_scores_variant', 'total_score_variant', 'mean_norm_score_non', 'sd_norm_score_non', 'norm_scores_non', - 'total_score_non', 'pvalue') + 'total_score_non', 'p_value_mean','p_value_min') regtools_data$sd_norm_score_variant[is.na(regtools_data$sd_norm_score_variant)] = 0 regtools_data$mean_norm_score_non[is.na(regtools_data$mean_norm_score_non)] = 0 regtools_data$sd_norm_score_non[is.na(regtools_data$sd_norm_score_non)] = 0 From 889e48ff07825bd80cc65968a6dd1a21caa62343 Mon Sep 17 00:00:00 2001 From: Kelsy Cotto Date: Sat, 11 Jan 2020 15:24:44 -0600 Subject: [PATCH 3/3] update docs and add argparse back in for stats script --- docs/workflow.md | 6 ++--- scripts/compare_junctions_hist_v2.R | 34 ++++++++++++++--------------- 2 files changed, 20 insertions(+), 20 deletions(-) diff --git a/docs/workflow.md b/docs/workflow.md index fc0126c..25d40d4 100644 --- a/docs/workflow.md +++ b/docs/workflow.md @@ -81,7 +81,7 @@ tabix all_variants_sorted.vcf.gz **Run `regtools cis-splice effects identify` on all samples with all variants (with `$tag` options as example)** -``` +```bash 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 all_variants_sorted.vcf.gz ${i}/tumor_rna_alignments.bam reference.fa reference.gtf; done ``` @@ -93,10 +93,10 @@ for i in samples/*/; do bsub -oo $i/logs/regtools_compare_$tag.lsf regtools cis- mkdir -p compare_junctions/hist ``` -**Run `compare_junctions_hist.R` on sample data** +**Run `stats_wrapper.py` on sample data** ```bash -Rscript --vanilla compare_junctions_hist.R +python3 stats_wrapper.py ``` **Run `filter_and_BH.R` to adjust p values and filter results** diff --git a/scripts/compare_junctions_hist_v2.R b/scripts/compare_junctions_hist_v2.R index 176bedd..2645d51 100644 --- a/scripts/compare_junctions_hist_v2.R +++ b/scripts/compare_junctions_hist_v2.R @@ -9,21 +9,21 @@ library(tidyverse) debug = F -# system.time({ -# if (debug){ -# tag = paste("_", "default", sep="") -# } else { -# # get options tag -# args = commandArgs(trailingOnly = TRUE) -# tag = args[1] -# input_file = args[2] -# if ( substr(tag, 2, 3) == "--"){ -# stop("Please specify an option tag (e.g. \"default\", \"i20e5\")") -# } -# } - -tag = 'default' -input_file = 'all_splicing_variants_default.bed' +system.time({ +if (debug){ + tag = paste("_", "default", sep="") +} else { + # get options tag + args = commandArgs(trailingOnly = TRUE) + tag = args[1] + input_file = args[2] + if ( substr(tag, 2, 3) == "--"){ + stop("Please specify an option tag (e.g. \"default\", \"i20e5\")") + } +} + +# tag = 'E' +# input_file = '/Users/kcotto/Desktop/CHOL/all_splicing_variants_E.bed' # All splicing relevant variants (union of rows from variants.bed files; add column with comma-separated list of sample names) all_splicing_variants = unique(data.table::fread(input_file), sep = '\t', header = T, stringsAsFactors = FALSE) @@ -308,5 +308,5 @@ regtools_data = regtools_data %>% distinct() write.table(regtools_data, file=paste(input_file, "_out.tsv", sep=""), quote=FALSE, sep='\t', row.names = F) -# -# }) + +})