diff --git a/.gitignore b/.gitignore
index e0d25df6cc8d76bda2280298fe657320fa3acc60..9d78a4386152bf03172641ac43f375e606d26068 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,3 +6,7 @@
# Python egg metadata, regenerated from source files by setuptools.
/*.egg-info
+
+# OS generated file
+
+.DS_Store
diff --git a/README.md b/README.md
index 867ca096ed249a853215b902be50fff38e71a438..3d3cddb6f261df52f9257475703b71de720c7236 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,28 @@
+__Changes since jSplice 1.1.0 (request to original git sent as issue)__
+This lists just a couple of quick fixes to bedtools coverage...
+
+* addition of `-g` for genomeFile in BEDtools when `-sorted` is used (for BEDtools > 2.24)
+* `--covgenome` argument - path to genome file for BEDtools genomeFile. To generate this please check the following biostars post (https://www.biostars.org/p/70795/) or
+
+```
+# fasta index from genome (generates file genome.fa.fai)
+
+samtools faidx genome.fa
+
+# take first two columns from index above
+# option 1: using `awk`
+
+awk -v OFS='\t' {'print $1,$2'} genome.fa.fai > genome.txt
+
+# or option 2: using `cut`. Note, .fai file should already be tab-delimited, but if needed can use `tr` command with cut
+
+cut -f1,2 genome.fa.fai > genome.txt
+```
+
+The resulting file ` <genome.txt>` from either `awk` or `cut` is used as input to the `--covgenome` argument.
+
+_____
jSplice 1.1.0
-------------
@@ -48,7 +72,7 @@ Required parameters
-------------------
`-o, --outdir`: Output directory. Any file already present in this directory will be overwritten. Ideally, the directory should identify your experiment. After execution, the directory will contain a jSplice object, a BED file of the trimmed genome, and all files generated by coverageBed (if any). Note that the directory will be used by runJSplice.py later on.
-`-d, --design`: Path to the experimental design file. The text file should contain one line per experiment and four columns: experiment, condition, junction file, and BAM file. If `–b` is set, then a BED file is expected in the fourth column instead of the BAM file. If `–j` is set, then the fourth column can be skipped.
+`-d, --design`: Path to the experimental design file. The text file should contain one line per experiment and four columns: experiment, condition, junction file, and BAM file. If `–b` is set, then a BED file is expected in the fourth column instead of the BAM file. If `–j` is set, then the fourth column can be skipped.
Important: file paths in the experimental design file can be either absolute or relative but cannot contain `~`.
This parameter is only required for the first run.
@@ -57,7 +81,7 @@ Optional parameters - controlling the thresholds
------------------------------------------------
`-c, --count`: Read count threshold. Any junction or exon that does not reach the threshold in at least one condition for all experiment is discarded. Default is 20.
-`-r, --relfc`: Relative fold-change threshold. The ratio difference of two element of an ASM is computed as the difference between the fold-change of condition A to B. The pair of elements that maximizes the fold-change difference (or relative fold-change) across all experiments is selected as representative of the ASM. The average log2 ratio difference of the representative pair is used to sort ASMs. If its value does not meet the threshold then the ASM is discarded. Note that a fold-change ratio is expected, not its log2 transform.
+`-r, --relfc`: Relative fold-change threshold. The ratio difference of two element of an ASM is computed as the difference between the fold-change of condition A to B. The pair of elements that maximizes the fold-change difference (or relative fold-change) across all experiments is selected as representative of the ASM. The average log2 ratio difference of the representative pair is used to sort ASMs. If its value does not meet the threshold then the ASM is discarded. Note that a fold-change ratio is expected, not its log2 transform.
(2 by default)
`-k, --rpkm`: RPKM threshold for gene expression. This parameter allows a simple filtering of non-expressed genes. To be retained, a gene has to be above the threshold in all experiments and all conditions. (1 by default)
diff --git a/jsplice/lib/jsplicelib.py b/jsplice/lib/jsplicelib.py
index f5962bd5b08359c67ac3b4479371b67a746baacf..d85aea66d657df9c2fb9c8eb68cb8cedf50a5bd9 100644
--- a/jsplice/lib/jsplicelib.py
+++ b/jsplice/lib/jsplicelib.py
@@ -38,19 +38,19 @@ def FDR(x):
return l
#Run coverageBed. Assumes that BAM files are sorted
-def runCovBed(q, bedFile, strand_arg):
+def runCovBed(q, bedFile, strand_arg, bed_genome_arg):
#get bedtools version
cmd = "coverageBed -h 2>&1 >/dev/null | grep Version | awk '{ print $2 }'"
p=subprocess.Popen(cmd,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out, _ = p.communicate()
bedtools_version = float('.'.join(out[1:].split('.')[:2]))
-
+
while True:
try:
exp = q.get(False)
except Queue.Empty:
break
-
+
if bedtools_version<2.24:
if strand_arg is None:
cmd='coverageBed -split -abam "'+exp.bamFile+'" -b "'+bedFile+'" > "'+exp.bedCountFile+'"'
@@ -58,10 +58,10 @@ def runCovBed(q, bedFile, strand_arg):
cmd='coverageBed '+strand_arg+' -split -abam "'+exp.bamFile+'" -b "'+bedFile+'" > "'+exp.bedCountFile+'"'
else: #starting on version 2.24, the a and b file are interverted...
if strand_arg is None:
- cmd='coverageBed -sorted -split -a "'+bedFile+'" -b "'+exp.bamFile+'" > "'+exp.bedCountFile+'"'
+ cmd='coverageBed -sorted -split -g '+bed_genome_arg+' -a "'+bedFile+'" -b "'+exp.bamFile+'" > "'+exp.bedCountFile+'"'
else:
- cmd='coverageBed '+strand_arg+' -sorted -split -a "'+bedFile+'" -b "'+exp.bamFile+'" > "'+exp.bedCountFile+'"'
-
+ cmd='coverageBed '+strand_arg+' -sorted -split -g '+bed_genome_arg+' -a "'+bedFile+'" -b "'+exp.bamFile+'" > "'+exp.bedCountFile+'"'
+
print cmd
p=subprocess.Popen(cmd,shell=True) #Not the safest solution but cannot find any other... :(
p.wait()
@@ -714,9 +714,9 @@ def correctAnnotations(junctions, genome):
for chrom in jxnChromList:
if chrom not in annChromList:
unmatchedChrom.append(chrom)
-
+
print 'Unmatched chrom. in annotation: '+str(unmatchedChrom)
-
+
if len(unmatchedChrom)>0:
print 'WARNING: Chromosome names in annotation file are different from the junctions\'.'
@@ -724,7 +724,7 @@ def correctAnnotations(junctions, genome):
print 'Junctions chrom. scheme: '+jxnScheme
annScheme = getChromosomeScheme(annChromList)
print 'Annotation chrom. scheme: '+annScheme
-
+
if jxnScheme is not None and annScheme is not None and jxnScheme!=annScheme:
print 'jSplice will attempt to correct the annotations by adding or removing "chr" in front of the chromosome name.'
@@ -910,7 +910,7 @@ def runPreparationStep(args):
flag = ExperimentalDesign.NOBAM
if args.jxnonly:
flag = ExperimentalDesign.JXNONLY
-
+
if not args.design:
print 'ERROR: No experimental design provided!'
sys.exit(0)
@@ -978,13 +978,16 @@ def runPreparationStep(args):
print 'New genome is saved in "'+os.path.join(args.outdir,'jSplice_genome.bed')+'"'
#Running coverageBed (if needed)
+
+ BED_genomeFile = args.covgenome
+
if not args.jxnonly and not args.nobam:
strand = None
if args.samestrand:
strand = '-s'
elif args.diffstrand:
strand = '-S'
-
+
if args.nbcores==0:
args.nbcores = len(expDesign)
print 'Running coverageBed on '+str(args.nbcores)+' threads.. (might take a long time)'
@@ -997,7 +1000,7 @@ def runPreparationStep(args):
tList = list()
for _ in range(args.nbcores):
- t = threading.Thread(target=runCovBed, args=(q,genomeFile,strand))
+ t = threading.Thread(target=runCovBed, args=(q,genomeFile,strand,BED_genomeFile))
t.daemon=True #Thread dies if the main process dies
tList.append(t)
t.start()
@@ -1102,11 +1105,3 @@ def runAnalysisStep(genome, expDesign, junctionCounts, exonCounts, totCount, jxn
printText(os.path.join(args.outdir,'jSplice_results.txt'),srtAvgRDiffs, rdiffs, asmFDRs)
print str(len(avgRDiffs))+' ASMs selected. (See file "'+os.path.join(args.outdir,'jSplice_results.html')+'" for details.)'
-
-
-
-
-
-
-
-
diff --git a/jsplice/run.py b/jsplice/run.py
index b23f2a1ff2c489cf990f9c09e0f3edacfe49b5ca..592b27c3ac1c3e75f0344634bd966fe8f750ef81 100644
--- a/jsplice/run.py
+++ b/jsplice/run.py
@@ -30,6 +30,7 @@ def main():
#User parameters
argparser = argparse.ArgumentParser(description='jSplice: a fast method to detect differential alternative splicing events.')
argparser.add_argument('-o','--outdir',help='Output directory', required=True)
+ argparser.add_argument('--covgenome',help='Path to genome file for BEDtools - Tmp fix.')
argparser.add_argument('-a','--annotation',help='GTF genome annotation file.')
argparser.add_argument('-d','--design', help='Experimental design filename. (Required for the first execution.)')
argparser.add_argument('-t','--type',help='Exon keyword in GTF annotation file (e.g. "exon" in Ensembl GTF files, default)',default='exon')
@@ -72,6 +73,13 @@ def main():
else:
os.makedirs(args.outdir)
+ # check if genomeFile for BEDtools exists if user provides one
+ # (which would be a requirement for use with -sorted with bedtools > 2.24. Would need to check prior versions...)
+
+ if args.covgenome != None:
+ if not os.path.isfile(args.covgenome):
+ sys.exit("---ERROR: Path to BEDtools genome file is not a valid directory, or file---")
+
#Check if first run
firstrun = not os.path.exists(os.path.join(args.outdir,'jSplice.dat'))
if firstrun and not args.design:
@@ -113,6 +121,5 @@ def main():
runAnalysisStep(genome, expDesign, jxnCounts, exnCounts, totCounts, jxn2gene, args)
-
if __name__ == "__main__":
- main()
\ No newline at end of file
+ main()