Graphing changes

pirl-unc · Oct 21, 2024 · 90c6830 · 90c6830
1 parent 7ee5620
commit 90c6830
Show file tree

Hide file tree

Showing 24 changed files with 228 additions and 131 deletions.
diff --git a/arriba_helper.sh b/arriba_helper.sh
@@ -5,17 +5,17 @@
 #SBATCH --time 12:00:00
 #SBATCH --mem=128G
 
-DATADIR=/datastore/scratch/users/vantwisk/sim/shortreads_training${6}k
-MAPPING_DIR=/datastore/scratch/users/vantwisk/sim/shortreads_training${6}k_mapping
-ARRIBA_DIR=/datastore/scratch/users/vantwisk/sim/shortreads_training${6}k_arriba
+MAPPING_DIR=${ALIGNMENT_STORAGE_DIR}/shortreads_${6}k_star
+ARRIBA_DIR=${ARRIBA_STORAGE_DIR}/shortreads_${6}k_arriba
 
 [ ! -d ${ARRIBA_DIR} ] && mkdir ${ARRIBA_DIR}
 
-#for j in $(seq 1 10); do
-#singularity exec --pid --bind /datastore arriba_latest.sif \
 /home/vantwisk/arriba_v2.2.1/arriba \
     -x ${MAPPING_DIR}/fusions-${1}-${4}-${5}-Aligned.sortedByCoord.out.bam \
     -o ${ARRIBA_DIR}/fusions-${1}-${4}-${5}.tsv -O ${ARRIBA_DIR}/fusions-${1}-${4}-${5}.discarded.tsv \
-    -a Homo_sapiens.GRCh38.dna.primary_assembly.fa -g Homo_sapiens.GRCh38.105.chr.gtf \
-    -b /home/vantwisk/arriba_v2.2.1/database/blacklist_hg38_GRCh38_v2.2.1.tsv -k /home/vantwisk/arriba_v2.2.1/database/known_fusions_hg38_GRCh38_v2.2.1.tsv -t /home/vantwisk/arriba_v2.2.1/database/known_fusions_hg38_GRCh38_v2.2.1.tsv -p /home/vantwisk/arriba_v2.2.1/database/protein_domains_hg38_GRCh38_v2.2.1.gff3
-#done
+    -a ${DNA_REFERENCE} \
+    -g ${GTF_REFERENCE} \
+    -b /home/vantwisk/arriba_v2.2.1/database/blacklist_hg38_GRCh38_v2.2.1.tsv \
+    -k /home/vantwisk/arriba_v2.2.1/database/known_fusions_hg38_GRCh38_v2.2.1.tsv \
+    -t /home/vantwisk/arriba_v2.2.1/database/known_fusions_hg38_GRCh38_v2.2.1.tsv \
+    -p /home/vantwisk/arriba_v2.2.1/database/protein_domains_hg38_GRCh38_v2.2.1.gff3
diff --git a/art_helper.sh b/art_helper.sh
@@ -6,7 +6,7 @@
 #SBATCH -t 02:00:00
 
 
-OUTDIR=shortreads
+OUTDIR=${SIM_STORAGE_DIR}/shortreads_${6}k
 [ ! -d ${OUTDIR} ] && mkdir ${OUTDIR}
 
 echo $1
@@ -15,5 +15,5 @@ echo $3
 echo $4
 
 #for j in $(seq 1 10); do
-/home/vantwisk/art_bin_MountRainier/art_illumina --rndSeed $RANDOM -ss HS25 -i ${3} -o ${OUTDIR}/fusions-${1}-${2}-${4}- -l ${2} -f ${4} -p -m 500 -s 10
+/home/vantwisk/art_bin_MountRainier/art_illumina --rndSeed $RANDOM -ss HS25 -i ${FUSION_TRANSCRIPTOME} -o ${OUTDIR}/fusions-${1}-${4}-${5}- -l ${2} -f ${1} -p -m 500 -s 10
 #done
diff --git a/badread_helper.sh b/badread_helper.sh
@@ -5,7 +5,7 @@
 #SBATCH -n 1
 #SBATCH --time 24:00:00
 
-OUTDIR=/datastore/scratch/users/vantwisk/sim/longreads_training${8}k
+OUTDIR=${SIM_STORAGE_DIR}/longreads_${8}k
 
 [ ! -d ${OUTDIR} ] && mkdir ${OUTDIR}
 echo $1
@@ -17,10 +17,12 @@ echo $6
 echo $7
 echo $8
 
-#for j in $(seq 1 10); do
-rustyread --threads 32 simulate --reference ${7} \
+rustyread --threads 32 simulate --reference ${FUSION_TRANSCRIPTOME} \
     --quantity ${1}x \
     --qscore_model /home/vantwisk/Badread/badread/qscore_models/${6} --glitches 0,0,0 --junk_reads 0 --random_reads 0 \
     --error_model /home/vantwisk/Badread/badread/error_models/${6} --identity ${5} \
-    --chimera 0 --seed $RANDOM | gzip > ${OUTDIR}/fusions-${1}-${5}-${6}-${4}.fq.gz #pfun/fuse-${1}-${4}.fq.gz  #fuse-transcript-${1}-${4}.fq.gz
-#done
+    --chimera 0 --seed $RANDOM > ${OUTDIR}/fusions-${1}-${5}-${6}-${4}.fq #pfun/fuse-${1}-${4}.fq.gz  #fuse-transcript-${1}-${4}.fq.gz
+
+awk '{ if (NR%4==1) gsub(".*","@transcript/"NR,$1); print }' ${OUTDIR}/fusions-${1}-${5}-${6}-${4}.fq > ${OUTDIR}/fusions-${1}-${5}-${6}-${4}_2.fq
+mv ${OUTDIR}/fusions-${1}-${5}-${6}-${4}_2.fq ${OUTDIR}/fusions-${1}-${5}-${6}-${4}.fq
+#awk '{ if (NR%4==1) gsub(".*","@transcript/"NR,$1); print }' ${OUTDIR}/fusions-${1}-${5}-${6}-${4}.fq > ${OUTDIR}/fusions-${1}-${5}-${6}-${4}.fq
diff --git a/combined_graphs.R b/combined_graphs.R
@@ -1,30 +1,29 @@
 ## Args 1 <- Simulated Fusions
 
-
-valid <- read.table('longread-fusion-transcript-pipeline/test_fusions1.txt', stringsAsFactors=F, header=T)
-colnames(valid) <- c('V1', 'V2')
-
-valid <- read.table('longread-fusion-transcript-pipeline/training1k_fusions.txt', stringsAsFactors=F, header=T)
-colnames(valid) <- c('V1', 'V2')
-
-coverage <- c('3', '5', '10', '30', '50', '100')
-identities <- c('87,97,5', '75,90,8', '95,100,4')
-tech <- c('pacbio2016') #, 'nanopore2020')
-#tech <- c(2, 5, 10) #c(10, 30, 50, 100, 300, 500)
-replicates <- 10
+arg <- commandArgs(trailingOnly=TRUE)
+arg[1] <- '/home/vantwisk/vantwisk/fusions/seq_run3/ref/fusim_ref_100.txt'
+arg[2] <- '/home/vantwisk/vantwisk/fusions/seq_run3/results/genion/longreads_1k_genion'
+arg[3] <- '/home/vantwisk/vantwisk/fusions/seq_run3/results/jaffal/longreads_1k_jaffal'
+arg[4] <- '/home/vantwisk/vantwisk/fusions/seq_run3/results/longgf/longreads_1k_longgf'
+
+valid <- read.table(arg[1], sep="\t", stringsAsFactors=F, header=T)
+
+coverage <- c('3', '10')
+identities <- c('87,97,5', '95,100,4')
+tech <- c('pacbio2016', 'nanopore2020')
+replicates <- 1
 #identities <- c("", "-minsup5", "_minsup10")
 
 fuseg <- lapply(coverage, function(i) {
   res1 <- lapply(identities, function(x) {
     res2 <- lapply(tech, function(j) {
       res3 <- lapply(1:replicates, function(r) {
-        filename <- paste0('/home/vantwisk/vantwisk/fusions/longread-fusion-transcript-pipeline/longreads_training1k_genion/fusions-', i,'-', x, '-', j, '-', r,'-genion-minsup-2.tsv')
+        filename <- paste0(arg[2], '/fusions-', i,'-', x, '-', j, '-', r,'-genion-minsup-1.tsv')
         message(filename)
         if (file.size(filename) == 0L){
           data.frame(tpos = 0, fpos=0, fneg=0, recall = 0, precision=0, coverage = i, quality=x, tech=j, replicate = r)
         } else {
         tab <- read.csv(filename, header=F, sep="\t", stringsAsFactors = F)
-        #tab <- tab[tab[,6] == "PASS:GF",]
         ss <- tab[,2]
         ss <- strsplit(ss, '::')
         taber <- do.call(rbind, ss)
@@ -40,10 +39,7 @@ fuseg <- lapply(coverage, function(i) {
         valid2['pass1'] <- ifelse(valid2[,1] %in% ll, 1, 0)
         valid2['pass2'] <- ifelse(valid2[,2] %in% ll, 1, 0)
         valid2['total'] <- ifelse(valid2['pass1'] + valid2['pass2'] > 0, 1, 0)
-        #dup <- duplicated(rbind(valid, tab))
-        #tab_total <- nrow(tab)
-        #total <- length(dup)
-        tpos <- sum(valid2["total"]) + 8
+        tpos <- sum(valid2["total"])
         fneg <- nrow(valid) - tpos
         fpos <- nrow(tab) - tpos
         fpos <- ifelse(fpos < 0, 0, fpos)
@@ -141,20 +137,13 @@ ggplot(fuseg, aes(x=coverage, y=tpos)) +
   # geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=.2,
   #               position=position_dodge(.9))
   ylab('mean true positives') + ggtitle('Nanopore2020 Mean found of control transcriptome with LongGF at coverage levels and read quality')
-valid <- read.table('training1k_fusions.txt', stringsAsFactors=F, header=T)
-colnames(valid) <- c('V1', 'V2')
 
-coverage <- c('5', '10', '30', '50', '100')
-identities <- c( '75,90,8', '87,97,5', '95,100,4')
-tech <- c('pacbio2016') # 'nanopore2020')
-replicates <- 10
-#identities <- c(0, 5, 10)
+if (FALSE){
 
 fusej <- lapply(coverage, function(i) {
   res1 <- lapply(identities, function(x) {
     res2 <- lapply(tech, function(j) {
-        filename <- paste0('/home/vantwisk/vantwisk/fusions/longread-fusion-transcript-pipeline/longreads1k_training_jaffal/fusions-', i,'-', x, '-', j,'-', 1,'-jaffal_out/jaffa_results.csv')
-        message(filename)
+        filename <- paste0(arg[3], '/fusions-', i,'-', x, '-', j,'-', 1,'-jaffal_out/jaffa_results.csv')
     tab <- read.csv(filename, sep=",", stringsAsFactors = F)
     tab <- tab[tab$spanning.reads > 2,]
     if(nrow(tab) == 0) {
@@ -175,10 +164,7 @@ fusej <- lapply(coverage, function(i) {
     valid2['pass1'] <- ifelse(valid2[,1] %in% ll, 1, 0)
     valid2['pass2'] <- ifelse(valid2[,2] %in% ll, 1, 0)
     valid2['total'] <- ifelse(valid2['pass1'] + valid2['pass2'] > 0, 1, 0)
-#    dup <- duplicated(rbind(valid, tab))
-#    tab_total <- nrow(tab)
-#    total <- length(dup)
-    tpos <- sum(valid2["total"]) + 6
+    tpos <- sum(valid2["total"])
     fneg <- nrow(valid) - tpos
     fpos <- nrow(tab) - tpos
     fpos <- ifelse(fpos < 0, 0, fpos)
@@ -266,7 +252,7 @@ fusec <- lapply(coverage, function(i) {
     res2 <- lapply(tech, function(j) {
       res3 <- lapply(1:replicates, function(r) {
         #res4 <- lapply(one, function(n) {
-        tab <- readLines(con=paste0('/home/vantwisk/vantwisk/fusions/longread-fusion-transcript-pipeline/longreads_training1k_longgf_control/fusions-', i,'-', x, '-', j,'-', r, '-', 100, '-', 50, '-', 100,'.log'))
+        tab <- readLines(con=paste0(arg[4], '/fusions-', i,'-', x, '-', j,'-', r, '-', 100, '-', 50, '-', 100,'.log'))
         sw <- startsWith(tab, 'GF')
         ss <- strsplit(tab[sw], ' ')
         ss <- vapply(ss, function(x) x[1], character(1))
@@ -321,7 +307,7 @@ fusel <- lapply(coverage, function(i) {
     res2 <- lapply(tech, function(j) {
       res3 <- lapply(1:replicates, function(r) {
         #res4 <- lapply(one, function(n) {
-  tab <- readLines(con=paste0('/home/vantwisk/vantwisk/fusions/longread-fusion-transcript-pipeline/longreads_training1k_longgf/fusions-', i,'-', x, '-', j, '-', r,'-', 100 ,'-', 50,'-', 100,'.log'))
+  tab <- readLines(con=paste0(arg[5], '/fusions-', i,'-', x, '-', j, '-', r,'-', 100 ,'-', 50,'-', 100,'.log'))
   sw <- startsWith(tab, 'GF')
   ss <- strsplit(tab[sw], ' ')
   ss <- vapply(ss, function(x) x[1], character(1))
@@ -955,3 +941,5 @@ ggplot(df, aes(x = coverage, y = precision, color = tool, group=tool)) +
   geom_errorbar(aes(ymin=precision-precision_sd, ymax=precision+precision_sd), width=2,
                 position=position_dodge(0.05)) +
   ggtitle("Longread and Short read tools coverage and precision")
+
+}
diff --git a/environment.sh b/environment.sh
@@ -1,12 +1,12 @@
 
 ## SYSTEM OPTIONS
-export TF_BASH=bash
+export TF_BASH=sbatch
 #export TF_BASH=sbatch
 
 export THREADS=32
 
 ## SCRATCH STORAGE DIRECTORY
-export STORAGE_DIR=/home/vantwisk/vantwisk/fusions/seq_run
+export STORAGE_DIR=/home/vantwisk/vantwisk/fusions/seq_run3
 [ ! -d ${STORAGE_DIR} ] && mkdir ${STORAGE_DIR}
 
 ## BASE STORAGE DIRECTORIES
@@ -21,6 +21,9 @@ export RESULTS_STORAGE_DIR=${STORAGE_DIR}/results
 export ALIGNMENT_STORAGE_DIR=${STORAGE_DIR}/alignments
 [ ! -d ${ALIGNMENT_STORAGE_DIR} ] && mkdir ${ALIGNMENT_STORAGE_DIR}
 
+
+export STAR_INDEX_STORAGE_DIR=${REF_STORAGE_DIR}/hg38_star_index
+
 ## RESULTS STORAGE DIRECTORIES
 export JAFFAL_STORAGE_DIR=${RESULTS_STORAGE_DIR}/jaffal
 [ ! -d ${JAFFAL_STORAGE_DIR} ] && mkdir ${JAFFAL_STORAGE_DIR}
@@ -43,9 +46,17 @@ export GRAPHS_STORAGE_DIR=${RESULTS_STORAGE_DIR}/graphs
 [ ! -d ${GRAPHS_STORAGE_DIR} ] && mkdir ${GRAPHS_STORAGE_DIR}
 
 ## RESOURCES
-export DNA_REFERENCE=${REF_STORAGE_DIR}/Homo_sapiens.GRCh38.dna.primary_assembly.fa
-export CDNA_REFERENCE=${REF_STORAGE_DIR}/Homo_sapiens.GRCh38.cdna.all.fa
-export GTF_REFERENCE=${REF_STORAGE_DIR}/Homo_sapiens.GRCh38.105.gtf
+export DNA_REFERENCE_GEN=${REF_STORAGE_DIR}/GRCh38.primary_assembly.genome.fa
+export CDNA_REFERENCE_GEN=${REF_STORAGE_DIR}/gencode.v38.transcripts.fa
+export GTF_REFERENCE_GEN=${REF_STORAGE_DIR}/gencode.v38.annotation.gtf
+
+export DNA_REFERENCE=${REF_STORAGE_DIR}/GRCh38.primary_assembly.genome.fa
+export CDNA_REFERENCE=${REF_STORAGE_DIR}/gencode.v38.transcripts.fa
+export GTF_REFERENCE=${REF_STORAGE_DIR}/gencode.v38.annotation.gtf
+
+export DNA_REFERENCE_ENS=${REF_STORAGE_DIR}/Homo_sapiens.GRCh38.dna.primary_assembly.fa
+export CDNA_REFERENCE_ENS=${REF_STORAGE_DIR}/Homo_sapiens.GRCh38.cdna.all.fa
+export GTF_REFERENCE_ENS=${REF_STORAGE_DIR}/Homo_sapiens.GRCh38.105.gtf
 
 export PBMM2_MMI=${REF_STORAGE_DIR}/hg38_gencode.mmi
 
@@ -55,13 +66,14 @@ export GENOMIC_SUPER_DUPS=${REF_STORAGE_DIR}/genomicSuperDups.txt
 
 ## Annotion Limit Settings
 export TRANSCRIPT_LIMIT=1000
-export TRANSCRIPT_LIMITED_FILE=${REF_STORAGE_DIR}/Homo_sapiens.GRCh38.cdna.limited_${TRANSCRIPT_LIMIT}.fa
+export TRANSCRIPT_LIMITED_FILE=${REF_STORAGE_DIR}/Homo_sapiens_transcriptome_limited_${TRANSCRIPT_LIMIT}.fa
 
 ## FUSION SIMULATION SETTINGS
 export NFUSIONS=100
 export FUSIM_FASTA_FILE=${REF_STORAGE_DIR}/fusim_${NFUSIONS}.fasta
-export FUSIM_TXT_FILE=${REF_STORAGE_DIR}/fusim_${NFUSIONS}.fxt
-export FUSION_TRANSCRIPTOME=${REF_STORAGE_DIR}/Homo_sapiens.GRCh38.cdna.limited_${TRANSCRIPT_LIMIT}_fusions_${NFUSIONS}.fa
+export FUSIM_TXT_FILE=${REF_STORAGE_DIR}/fusim_${NFUSIONS}.txt
+export FUSIM_REF_FILE=${REF_STORAGE_DIR}/fusim_ref_${NFUSIONS}.txt
+export FUSION_TRANSCRIPTOME=${REF_STORAGE_DIR}/Homo_sapiens.GRCh38.cdna.limited_${TRANSCRIPT_LIMIT}_fusions_${NFUSIONS}_2.fa
 
 ## LONGREAD AND SHORTREAD READ SIMULATION SETTINGS
 export N_TRANSCRIPTS=('1') #('1' '2' '4' '8' '16')

diff --git a/fusionseeker_helper.sh b/fusionseeker_helper.sh
@@ -1,20 +1,27 @@
 #!/bin/bash
-#SBATCH --job-name longgf
+#SBATCH --job-name fusionseeker
 #SBATCH --partition allnodes
 #SBATCH --time UNLIMITED
 #SBATCH --cpus-per-task 1
 #SBATCH --mem-per-cpu 10g
 
-DATADIR_MINIMAP=/datastore/scratch/users/vantwisk/sim/longreads_training${12}k_minimap2
-FUSIONSEEKER_DIR=/datastore/scratch/users/vantwisk/sim/longreads_training${12}k_fusionseeker
+DATADIR_MINIMAP=${ALIGNMENT_STORAGE_DIR}/longreads_${9}k_minimap2
+FUSIONSEEKER_DIR=${FUSIONSEEKER_STORAGE_DIR}/longreads_${9}k_fusionseeker
 
 [ ! -d ${FUSIONSEEKER_DIR} ] && mkdir ${FUSIONSEEKER_DIR}
 
-#singularity exec --pid --bind /datastore longgf_0.1.2--h05f6578_1.sif \
-fusionseeker \
-  --tread 16 \
+/home/vantwisk/FusionSeeker/fusionseeker \
+  --thread 32 \
   --bam ${DATADIR_MINIMAP}/fusions-${1}-${5}-${6}-${4}-sorted.bam \
-  --gtf Homo_sapiens.GRCh38.105.gtf \
-  --ref ../hg38.fa \
-  -o ${FUSIONSEEKER_DIR}/fusions-${1}-${5}-${6}-${4}-fusionseeker \
+  --datatype nanopore \
+  --human38 \
+  --outpath ${FUSIONSEEKER_DIR}/fusions-${1}-${5}-${6}-${4}-fusionseeker \
   -s 2
+
+#fusionseeker \
+#  --thread 16 \
+#  --bam ${DATADIR_MINIMAP}/fusions-${1}-${5}-${6}-${4}-sorted.bam \
+#  --gtf ${GTF_REFERENCE} \
+#  --ref ${DNA_REFERENCE} \
+#  -o ${FUSIONSEEKER_DIR}/fusions-${1}-${5}-${6}-${4}-fusionseeker \
+#  -s 2
diff --git a/generate_annotation_resources.sh b/generate_annotation_resources.sh
@@ -1,16 +1,23 @@
 
-if [ ! -f ${DNA_REFERENCE} ]; then
-  wget -O ${DNA_REFERENCE}.gz http://ftp.ensembl.org/pub/release-105/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
-  gunzip -d ${DNA_REFERENCE}.gz
+if [ ! -f ${DNA_REFERENCE_GEN} ]; then
+  wget -O ${DNA_REFERENCE_GEN}.gz https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_38/GRCh38.primary_assembly.genome.fa.gz
+  gunzip -d ${DNA_REFERENCE_GEN}.gz
+fi
+
+if [ ! -f ${DNA_REFERENCE_ENS} ]; then
+  wget -O ${DNA_REFERENCE_ENS}.gz http://ftp.ensembl.org/pub/release-105/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
+  gunzip -d ${DNA_REFERENCE_ENS}.gz
 fi
 
 if [ ! -f ${CDNA_REFERENCE} ]; then
-  wget -O ${CDNA_REFERENCE}.gz http://ftp.ensembl.org/pub/release-105/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz
+  #wget -O ${CDNA_REFERENCE}.gz http://ftp.ensembl.org/pub/release-105/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz
+  wget -O ${CDNA_REFERENCE}.gz https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_38/gencode.v38.transcripts.fa.gz
   gunzip -d ${CDNA_REFERENCE}.gz
 fi
 
 if [ ! -f ${GTF_REFERENCE} ]; then
-  wget -O ${GTF_REFERENCE}.gz http://ftp.ensembl.org/pub/release-105/gtf/homo_sapiens/Homo_sapiens.GRCh38.105.gtf.gz
+  #wget -O ${GTF_REFERENCE}.gz http://ftp.ensembl.org/pub/release-105/gtf/homo_sapiens/Homo_sapiens.GRCh38.105.gtf.gz
+  wget -O ${GTF_REFERENCE}.gz https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_38/gencode.v38.annotation.gtf.gz
   gunzip -d ${GTF_REFERENCE}.gz
 fi
 
@@ -20,5 +27,5 @@ if [ ! -f ${GENOMIC_SUPER_DUPS} ]; then
 fi
 
 if [ ! -f ${TRANSCRIPT_LIMITED_FILE} ]; then
-  Rscript generate_breakpoints.R ${GTF_REFERENCE} ${CDNA_REFERENCE} ${TRANSCRIPT_LIMITED_FILE} ${TRANSCRIPT_LIMIT}
+  Rscript generate_breakpoints.R ${GTF_REFERENCE_ENS} ${CDNA_REFERENCE_ENS} ${TRANSCRIPT_LIMITED_FILE} ${TRANSCRIPT_LIMIT}
 fi
diff --git a/generate_arriba.sh b/generate_arriba.sh
@@ -1,9 +1,15 @@
+N_TRANSCRIPTS=('1') #('1' '2' '4' '8' '16')
+REPLICATES=1 #10
+COVERAGE=(3 10) #(3 5 10 30 50 100)
+QUALITY=('87,97,5' '95,100,4')  #('75,90,8' '87,97,5'  '95,100,4')
+TECH=('pacbio2016' 'nanopore2020')
+READ_LENGTHS=(100 150)
 
 for i in $(seq 1 ${REPLICATES}); do
   for q in ${!COVERAGE[@]}; do
     for j in ${!READ_LENGTHS[@]}; do
       for n in ${!N_TRANSCRIPTS[@]}; do
-        eval ${TF_BASH} arriba_helper.sh ${COVERAGE[$q]} 1 1 ${READ_LENGTHS[$j]} ${i} ${N_TRANSCRIPTS[$n]}
+        sbatch arriba_helper.sh ${COVERAGE[$q]} 1 1 ${READ_LENGTHS[$j]} ${i} ${N_TRANSCRIPTS[$n]}
       done
     done
   done

diff --git a/generate_breakpoints.R b/generate_breakpoints.R
@@ -10,6 +10,11 @@ library(Biostrings)
 
 arg <- commandArgs(trailingOnly=TRUE)
 
+arg[1] <- '/home/vantwisk/vantwisk/fusions/seq_run2/ref/gencode.v38.annotation.gtf'
+arg[2] <- '/home/vantwisk/vantwisk/fusions/seq_run2/ref/gencode.v38.transcripts.fa'
+arg[3] <- '/home/vantwisk/vantwisk/fusions/seq_run2/ref/Homo_sapiens_transcriptome_limited_1000.fa'
+arg[4] <- 1000
+
 message(arg[1])
 message(arg[2])
 message(arg[3])
@@ -44,6 +49,7 @@ vals1 <- ele[ele$full_tx %in% tx1,]
 wh1 <- which(tx_act %in% ele$full_tx)
 fa1 <- fa[wh1]
 fa1 <- fa1[lengths(fa1) > 100]
+
 fa1 <- sample(fa1, arg[4], replace=F)
 
 writeXStringSet(fa1, fasta_out)

diff --git a/generate_fusim.sh b/generate_fusim.sh
@@ -0,0 +1,23 @@
+
+if [ ! -f ${REF_STORAGE_DIR}/refFlat.txt ]; then
+  wget -O ${REF_STORAGE_DIR}/refFlat.txt.gz http://hgdownload.cse.ucsc.edu/goldenPath/hg38/database/refFlat.txt
+  gunzip ${REF_STORAGE_DIR}/refFlat.txt.gz
+fi
+
+if [ ! -f ${REF_STORAGE_DIR}/hg38_chroma.fa ]; then
+  wget -O ${REF_STORAGE_DIR}/hg38.chromFa.tar.gz ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.chromFa.tar.gz
+  tar -xzf ${REF_STORAGE_DIR}/hg38.chromFa.tar.gz
+  cat ${REF_STORAGE_DIR}/chr*.fa > hg38_chroma.fa
+  samtools faidx ${REF_STORAGE_DIR}/hg38_chroma.fa
+fi
+
+java -jar /home/vantwisk/fusim-0.2.2/fusim.jar \
+  --gene-model=${REF_STORAGE_DIR}/refFlat.txt \
+  --fusions=${NFUSIONS} \
+  --reference=${REF_STORAGE_DIR}/hg38_chroma.fa \
+  --fasta-output=${FUSIM_FASTA_FILE} \
+  --text-output=${FUSIM_TXT_FILE}
+
+cat ${TRANSCIPT_LIMITED_FILE} ${FUSIM_FASTA_FILE} > ${FUSION_TRANSCRIPTOME}
+
+Rscript create_fusim_ref.R ${FUSIM_TXT_FILE} ${FUSIM_REF_FILE}
diff --git a/generate_fusionseeker.sh b/generate_fusionseeker.sh
@@ -1,10 +1,16 @@
+N_TRANSCRIPTS=('1') #('1' '2' '4' '8' '16')
+REPLICATES=1 #10
+COVERAGE=(3 10) #(3 5 10 30 50 100)
+QUALITY=('87,97,5' '95,100,4')  #('75,90,8' '87,97,5'  '95,100,4')
+TECH=('pacbio2016' 'nanopore2020')
+READ_LENGTHS=(100 150)
 
 for i in $(seq 1 ${REPLICATES}); do
   for q in ${!COVERAGE[@]}; do
     for j in ${!QUALITY[@]}; do
       for k in ${!TECH[@]}; do
         for n in ${!N_TRANSCRIPTS[@]}; do
-           eval ${TF_BASH} fusionseeker_helper.sh ${COVERAGE[$q]} 1 1 ${i} ${QUALITY[$j]} ${TECH[$k]} ax sam ${MIN_OVERLAP_LEN} ${BIN_SIZE} ${MIN_MAP_LENGTH} ${N_TRANSCRIPTS[$n]}
+           eval ${TF_BASH} fusionseeker_helper.sh ${COVERAGE[$q]} 1 1 ${i} ${QUALITY[$j]} ${TECH[$k]} ax sam ${N_TRANSCRIPTS[$n]}
         done
       done
     done