diff --git a/LICENSE.md b/LICENSE.md
index 33e2553..7e92416 100644
--- a/LICENSE.md
+++ b/LICENSE.md
@@ -1,4 +1,4 @@
- Copyright (c) Minh Duc Cao, Monash Uni & UQ, All rights reserved.
+ Copyright (c) the author(s), All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
diff --git a/Makefile b/Makefile
index b8e0687..3462121 100644
--- a/Makefile
+++ b/Makefile
@@ -52,7 +52,8 @@ JAR_FILE=japsa.jar
LIB_DIR=libs
#List of external libraries
-EXT_LIBS=colloquial.jar commons-math3-3.0.jar jhdf5.jar jhdfobj.jar htsjdk-1.126.jar guava-18.0.jar jcommon-1.0.23.jar jfreechart-1.0.19.jar JRIEngine.jar JRI.jar
+##EXT_LIBS=colloquial.jar commons-math3-3.0.jar jhdf5.jar jhdfobj.jar htsjdk-2.10.1.jar guava-18.0.jar jcommon-1.0.23.jar jfreechart-1.0.19.jar JRIEngine.jar JRI.jar
+EXT_LIBS=colloquial.jar commons-math3-3.0.jar jhdf5.jar jhdfobj.jar htsjdk-2.10.1.jar guava-18.0.jar jcommon-1.0.23.jar jfreechart-1.0.19.jar JRIEngine.jar JRI.jar gs-core-1.3.jar gs-ui-1.3.jar gs-algo-1.3.jar javax.json-api-1.0.jar javax.json-1.0.4.jar slf4j-api-1.7.25.jar slf4j-simple-1.7.25.jar jfxrt.jar
###########################################################################
##What this scripts does:
@@ -97,7 +98,11 @@ classes: $(CLASS_FILES)
$(CLASS_DIR):
mkdir -p $(CLASS_DIR)
-$(JAR_FILE):classes
+images:
+ cp -r $(SRC_DIR)/japsa/bio/misc/dnaPlatform/gui/images $(CLASS_DIR)/japsa/bio/misc/dnaPlatform/gui/
+ cp -r $(SRC_DIR)/japsa/seq/nanopore/icons $(CLASS_DIR)/japsa/seq/nanopore/
+
+$(JAR_FILE):classes images
jar cf $(JAR_FILE) -C $(CLASS_DIR) .
jar: $(JAR_FILE)
@@ -132,19 +137,29 @@ else
O_MXMEM=
endif
+
+ifdef SERVER
+O_SERVER=--server=${SERVER}
+else
+O_SERVER=
+endif
+
+docs: jar
+ java -cp $(JAR_FILE):$(LIB_DIR)/guava-18.0.jar:$(LIB_DIR)/slf4j-api-1.7.25.jar:$(LIB_DIR)/slf4j-simple-1.7.25.jar japsa.util.deploy.GenDocs && cd docs && make html
+
RELEASE=JapsaRelease
pre-install: jar
- @@echo "java -cp $(JAR_FILE):$(LIB_DIR)/guava-18.0.jar japsa.util.deploy.Deploy --mode install --libs $(subst $(SPACE),:, $(EXT_LIBS)) $(O_INS_DIR) $(O_JLP) $(O_MXMEM) --compiler \"`$(JCC) -version 2>&1`\"" > install.sh && \
+ @@echo "java -cp $(JAR_FILE):$(LIB_DIR)/guava-18.0.jar:$(LIB_DIR)/slf4j-api-1.7.25.jar:$(LIB_DIR)/slf4j-simple-1.7.25.jar japsa.util.deploy.Deploy --mode install --libs $(subst $(SPACE),:, $(EXT_LIBS)) $(O_INS_DIR) $(O_JLP) $(O_MXMEM) $(O_SERVER) --compiler \"`$(JCC) -version 2>&1`\"" > install.sh && \
chmod u+x install.sh && \
- echo "java -cp $(JAR_FILE);$(LIB_DIR)\guava-18.0.jar japsa.util.deploy.Deploy --mode install --libs $(subst $(SPACE),:, $(EXT_LIBS)) $(O_INS_DIR) $(O_JLP) $(O_MXMEM) --compiler \"`$(JCC) -version 2>&1`\"" > install.bat && \
+ echo "java -cp $(JAR_FILE);$(LIB_DIR)\guava-18.0.jar;$(LIB_DIR)\slf4j-api-1.7.25.jar;$(LIB_DIR)\slf4j-simple-1.7.25.jar japsa.util.deploy.Deploy --mode install --libs $(subst $(SPACE),:, $(EXT_LIBS)) $(O_INS_DIR) $(O_JLP) $(O_MXMEM) $(O_SERVER) --compiler \"`$(JCC) -version 2>&1`\"" > install.bat && \
echo "Installation scripts created"
install: pre-install
./install.sh
uninstall:
- @@java -cp $(JAR_FILE):$(LIB_DIR)/guava-18.0.jar japsa.util.deploy.Deploy --mode uninstall --libs $(subst $(SPACE),:, $(EXT_LIBS)) ${O_INS_DIR} && echo "Japsa uninstalled!"
+ @@java -cp $(JAR_FILE):$(LIB_DIR)/guava-18.0.jar:$(LIB_DIR)/slf4j-api-1.7.25.jar:$(LIB_DIR)/slf4j-simple-1.7.25.jar japsa.util.deploy.Deploy --mode uninstall --libs $(subst $(SPACE),:, $(EXT_LIBS)) ${O_INS_DIR} && echo "Japsa uninstalled!"
release: pre-install
@@mkdir -p $(RELEASE)/libs/ && \
diff --git a/README.md b/README.md
index 8313ddc..ad7c72b 100644
--- a/README.md
+++ b/README.md
@@ -19,18 +19,37 @@ Details of installation (including for Windows) and usage of Japsa can be found
in its documentation hosted on [ReadTheDocs](http://japsa.readthedocs.org/en/latest/index.html)
### Authors and Contributors
-Japsa was mainly developed and curated by Minh Duc Cao (@mdcao). The following
-people have contributed to the development of Japsa, including ideas,
+Japsa was maintained by Minh Duc Cao (@mdcao). The following
+people (in alphatical order) have contributed to the development of Japsa, including ideas,
algorithms, implementation, documentation and feedback:
-* Lachlan Coin
-* Son Hoang Nguyen
-* Mikael Boden
-* Lloyd Allison
-* Trevor I Dix
+* Bhuvan Sankar
+* David Powell
+* Devika Ganesamoorthy
* Hoang Anh Nguyen
* Julia Bernal
-* David Powell
+* [Lachlan Coin](http://www.imb.uq.edu.au/lachlan-coin)
+* [Lloyd Allison](http://www.allisons.org/ll/)
+* Michael Hall
+* Mikael Boden
+* [Son Hoang Nguyen](https://github.com/hsnguyen)
+* Trevor I Dix
+
+
+### Other projects based on Japsa
+
+
+* [eXpert Model](https://github.com/mdcao/xm): The expert model compression model
+* [XMas](https://github.com/mdcao/XMas): Phylogenetic distance method using information theory
+* [capsim](https://github.com/mdcao/capsim): Simulation of capture sequencing
+* [npScarf](https://github.com/mdcao/npScarf): Scaffold and Complete assemblies in real-time fashion
+* [npAnalysis](https://github.com/mdcao/npAnalysis): Realtime identification of bacterial sample
+* [npReader](https://github.com/mdcao/npReader): Real-time extraction and analysis Oxford Nanopore sequencing data
+* [npBarcode](https://github.com/hsnguyen/npBarcode): Demultiplex barcoded Oxford Nanopore sequencing
+* [PhageXpress](https://github.com/mdcao/phagexpress)
+
+and more to come.
+
### License
Japsa is released under the accompanying BSD-like license.
diff --git a/bin/.gitignore b/bin/.gitignore
new file mode 100644
index 0000000..efaee1f
--- /dev/null
+++ b/bin/.gitignore
@@ -0,0 +1,2 @@
+/japsa/
+/japsadev/
diff --git a/docs/source/install.rst b/docs/source/install.rst
index d07a0ee..b488392 100644
--- a/docs/source/install.rst
+++ b/docs/source/install.rst
@@ -4,7 +4,7 @@ Installation
There are two methods to install Japsa in your computer. The first method
-(using pre-compiled package in JDK 1.6) is straight-forward and can be used for
+(using pre-compiled package in JDK 1.8) is straight-forward and can be used for
any operating systems, including Windows. The second method (compile from source
code) requires some extra tools (make and JDK) but may yield better runtime
performance as the package will be compiled with the same version of the Java
@@ -43,7 +43,9 @@ agree with its suggestion, just type Enter. The questions are:
HDFViewer (https://www.hdfgroup.org/products/java/release/download.html)
installed, and enter the path to file *libjhdf5.so* (on Linux/Unix/Mac) or
to *jhdf5.dll* (Windows). This is only required if you intend to use npReader(
- jsa.np.f5reader).
+ jsa.np.f5reader). Note that we tested with hdfj version 2.10.1, which you can
+ download from https://support.hdfgroup.org/ftp/HDF5/releases/HDF-JAVA/hdf-java-2.10.1/bin/
+
* *Path to JRI library:* Enter path to JRI library. This is required only for
running real-time species typing (jsa.np.speciesTyping).
diff --git a/docs/source/tools/index.rst b/docs/source/tools/index.rst
index 9158a85..b48509b 100644
--- a/docs/source/tools/index.rst
+++ b/docs/source/tools/index.rst
@@ -16,13 +16,20 @@ This chapter presents the list of tools provided by Japsa.
jsa.seq.annovcf.rst
jsa.seq.gff2fasta.rst
jsa.seq.emalign.rst
- jsa.np.f5reader.rst
+ jsa.hts.countReads.rst
+ jsa.hts.errorAnalysis.rst
+ jsa.hts.n50.rst
+ jsa.np.npreader.rst
jsa.np.filter.rst
jsa.np.rtSpeciesTyping.rst
jsa.np.rtMLST.rst
jsa.np.rtStrainTyping.rst
jsa.np.rtResistGenes.rst
+ jsa.np.npscarf.rst
+ jsa.np.barcode.rst
jsa.util.streamServer.rst
jsa.util.streamClient.rst
jsa.phylo.xmas.rst
jsa.phylo.normalise.rst
+ jsa.sim.capsim.rst
+ jsa.xm.compress.rst
diff --git a/docs/source/tools/jsa.hts.countReads.rst b/docs/source/tools/jsa.hts.countReads.rst
new file mode 100644
index 0000000..dc42dea
--- /dev/null
+++ b/docs/source/tools/jsa.hts.countReads.rst
@@ -0,0 +1,49 @@
+------------------------------------------------
+*jsa.hts.countReads*: Count reads from bam files
+------------------------------------------------
+
+
+~~~~~~~~
+Synopsis
+~~~~~~~~
+
+*jsa.hts.countReads*: Count the number of reads in some regions from a sorted, indexed bam file
+
+~~~~~
+Usage
+~~~~~
+::
+
+ jsa.hts.countReads [options]
+
+~~~~~~~
+Options
+~~~~~~~
+ --bamFile=s Name of the bam file
+ (REQUIRED)
+ --bedFile=s Name of the regions file in bed format
+ (REQUIRED)
+ --output=s Name of output file, - for from standard out.
+ (default='-')
+ --flanking=i Size of the flanking regions, effectively expand the region by flanking
+ (default='0')
+ --qual=i Minimum quality
+ (default='0')
+ --filterBits=i Filter reads based on flag. Common values:
+ 0 no filter
+ 256 exclude secondary alignment
+ 1024 exclude PCR/optical duplicates
+ 2048 exclude supplementary alignments
+ (default='0')
+ --contained Count reads contained in the region
+ (default='false')
+ --overlap Count number of read overlap with the region
+ (default='false')
+ --span Count reads span the region
+ (default='false')
+ --help Display this usage and exit
+ (default='false')
+
+
+
+
diff --git a/docs/source/tools/jsa.hts.errorAnalysis.rst b/docs/source/tools/jsa.hts.errorAnalysis.rst
new file mode 100644
index 0000000..275d0b7
--- /dev/null
+++ b/docs/source/tools/jsa.hts.errorAnalysis.rst
@@ -0,0 +1,39 @@
+----------------------------------------------------------
+*jsa.hts.errorAnalysis*: Error analysis of sequencing data
+----------------------------------------------------------
+
+*jsa.hts.errorAnalysis* assesses the error profile of sequencing data by getting the numbers
+of errors (mismatches, indels etc) from a bam file. Obviously, it does not distinguish
+sequencing errors from mutations, and hence consider mutations as errors. It is best to use
+with the bam file from aligning sequencing reads to a reliable assembly of the sample.
+
+~~~~~~~~
+Synopsis
+~~~~~~~~
+
+*jsa.hts.errorAnalysis*: Error analysis of sequencing data
+
+~~~~~
+Usage
+~~~~~
+::
+
+ jsa.hts.errorAnalysis [options]
+
+~~~~~~~
+Options
+~~~~~~~
+ --bamFile=s Name of bam file
+ (REQUIRED)
+ --reference=s Name of reference genome
+ (REQUIRED)
+ --pattern=s Pattern of read name, used for filtering
+ (default='null')
+ --qual=i Minimum quality required
+ (default='0')
+ --help Display this usage and exit
+ (default='false')
+
+
+
+
diff --git a/docs/source/tools/jsa.hts.n50.rst b/docs/source/tools/jsa.hts.n50.rst
new file mode 100644
index 0000000..bac2936
--- /dev/null
+++ b/docs/source/tools/jsa.hts.n50.rst
@@ -0,0 +1,28 @@
+-----------------------------------------
+*jsa.hts.n50*: Compute N50 of an assembly
+-----------------------------------------
+
+~~~~~~~~
+Synopsis
+~~~~~~~~
+
+*jsa.hts.n50*: Compute N50 of an assembly
+
+~~~~~
+Usage
+~~~~~
+::
+
+ jsa.hts.n50 [options]
+
+~~~~~~~
+Options
+~~~~~~~
+ --input=s Name of the file
+ (REQUIRED)
+ --help Display this usage and exit
+ (default='false')
+
+
+
+
diff --git a/docs/source/tools/jsa.np.barcode.rst b/docs/source/tools/jsa.np.barcode.rst
new file mode 100644
index 0000000..52d1a31
--- /dev/null
+++ b/docs/source/tools/jsa.np.barcode.rst
@@ -0,0 +1,128 @@
+---------------------------------------------------------------------------
+*barcode*: real-time de-multiplexing Nanopore reads from barcode sequencing
+---------------------------------------------------------------------------
+
+*barcode* (jsa.np.barcode) is a program that demultiplex the nanopore reads from
+Nanopore barcode sequencing. Downstream analysis can be invoked concurrently by an input script.
+
+*barcode* is included in the `Japsa package `_.
+
+~~~~~~~~
+Synopsis
+~~~~~~~~
+
+*jsa.np.barcode*: Clustering nanopore sequences based on barcode
+
+~~~~~
+Usage
+~~~~~
+::
+
+ jsa.np.barcode [options]
+
+~~~~~~~
+Options
+~~~~~~~
+ --bcFile=s Barcode file
+ (REQUIRED)
+ --seqFile=s Nanopore sequences file
+ (REQUIRED)
+ --scriptRun=s Invoke command script to run npScarf
+ (default='null')
+ --threshold=d Minimum identity(%) for barcode alignment
+ (default='70.0')
+ --distance=d Minimum identity(%) distance between the best alignment to others
+ (default='4.0')
+ --twoends Whether a read must contain barcode sequence from both ends or just one end (default)
+ (default='false')
+ --print Print out demultiplexed reads to corresponding FASTA file or not.
+ (default='false')
+ --help Display this usage and exit
+ (default='false')
+
+
+
+
+~~~~~~~~~~~~~~
+Usage examples
+~~~~~~~~~~~~~~
+
+A summary of *barcode* usage can be obtained by invoking the --help option::
+
+ jsa.np.barcode --help
+
+=====
+Input
+=====
+ *barcode* takes 2 files as required input::
+
+ jsa.np.barcode -seq -bc
+
+<*nanopore reads*> is either the long reads in FASTA/FASTQ file (after MinION sequencing is
+finished) or standard input ( specified by "-", for real-time analysis).
+
+<*barcode.fasta*> is the FASTA file of barcode sequences (given by ONT) with name correspond to the assigned sample id.
+
+Missing any file would break down the whole pipeline.
+
+In addition, one can provide <*analysis_script*> which is the script call for further action on the de-multiplexed reads. It always take one argument and be
+executable by invoking::
+
+ ./analysis_script
+
+in which <*id*> is the identifier of a sample as given in the . The script should read the standard input
+of long-read streams to do further analysis.
+
+*barcode* allows user to set the minimum criteria of a hit with barcode reference to be considered valid. The default value
+is 70% for minimum identity. At the same time, 4% distance between the best hit and the second best is necessary for differentiation.
+Decreasing the thresholds will lead to more reads being clustered but with higher risk of false positive while more stringent parameters
+will generate less but more confident of demultiplexed reads.
+
+User can also have control on the matching condition for barcode detection, either one-end match or both-end match. For the first case (default), only the
+a legal maximal hit from one end of a read is enough to label it while in the later case, we take into account a pair from both 5' and 3'terminus.
+Thus the input for each use case should be different. The one-end option can take the simple FASTA file of Nanopore barcodes while the two-end need pairs of
+barcode to be specified (e.g. with _F and _R suffix). One of a typical use case for two-end matching is when we want to detect the super-barcode which includes
+also tail- and primer-sequences in pre-defined orientation.
+
+======
+Output
+======
+*barcode* output depends on the <*analysis script*> because the de-multiplexed reads are streamed directly to its dedicated process.
+If ones only interest in de-multiplexing alone, then the script should be as simple as to write stream to file. For example:
+
+.. code-block:: bash
+ :linenos:
+
+ #!/bin/bash
+ while read line
+ do
+ echo "$line"
+ done >> ${1}_script.fasta
+
+This is equivalent to enable the *-p* option::
+
+ jsa.np.barcode -seq -bc -script -p
+
+that would print out de-multiplexed FASTA sequences \_clustered.fasta
+
+============================================
+Real-time scaffolding for barcode sequencing
+============================================
+
+One use-case for barcode sequencing is to run *npscarf* on the resulted de-multiplexed reads. This could be done by calling a script
+that can take an output folder of long reads from a sample to scaffold its corresponding short-reads (e.g. SPAdes) assembly.
+E.g.
+
+.. code-block:: bash
+ :linenos:
+
+ #!/bin/bash
+ dirname=`find /coin/barcode/ -maxdepth 1 -type d -name "*${1}*" -print -quit`
+
+ bwa index ${dirname}/contigs.fasta
+
+ bwa mem -t 16 -k11 -W20 -r10 -A1 -B1 -O1 -E1 -L0 -a -Y -K 10000 ${dirname}/contigs.fasta - 2> /dev/null | \
+ jsa.np.npscarf -realtime -read 100 -time 1 -b - -seq ${dirname}/contigs.fasta -spadesDir ${dirname} -prefix ${1} > ${1}.log 2>&1
+
+In this scenario, we assume the output SPAdes folders locate in one directory and the folder names contain the ID of the corresponding samples.
+
diff --git a/docs/source/tools/jsa.np.filter.rst b/docs/source/tools/jsa.np.filter.rst
index 7dbf187..8e1fd5a 100644
--- a/docs/source/tools/jsa.np.filter.rst
+++ b/docs/source/tools/jsa.np.filter.rst
@@ -3,7 +3,7 @@
---------------------------------------
*jsa.np.filter* filters sequencing data based on sequence read type, length and
-quality. Examples of its usage can be found on jsa.np.f5reader_.
+quality. Examples of its usage can be found on jsa.np.npreader_.
~~~~~~~~
Synopsis
@@ -33,6 +33,8 @@ Options
(default='0.0')
--qualMax=d Maximum average quality
(default='1000.0')
+ --group=s Group need to be extracted, leave blank for selecting all groups
+ (default='')
--excl2D Exclude 2D reads
(default='false')
--exclTemp Exclude template reads
@@ -49,9 +51,9 @@ Options
See also
~~~~~~~~
-jsa.np.f5reader_, jsa.util.streamServer_, jsa.util.streamClient_
+jsa.np.npreader_, jsa.util.streamServer_, jsa.util.streamClient_
-.. _jsa.np.f5reader: jsa.np.f5reader.html
+.. _jsa.np.npreader: jsa.np.npreader.html
.. _jsa.util.streamServer: jsa.util.streamServer.html
.. _jsa.util.streamClient: jsa.util.streamClient.html
diff --git a/docs/source/tools/jsa.np.f5reader.rst b/docs/source/tools/jsa.np.npreader.rst
similarity index 80%
rename from docs/source/tools/jsa.np.f5reader.rst
rename to docs/source/tools/jsa.np.npreader.rst
index 7c7af5b..911ffc3 100644
--- a/docs/source/tools/jsa.np.f5reader.rst
+++ b/docs/source/tools/jsa.np.npreader.rst
@@ -2,7 +2,7 @@
*npReader*: real-time conversion and analysis of Nanopore sequencing data
-------------------------------------------------------------------------
-*npReader* (jsa.np.f5reader) is a program that extracts Oxford Nanopore
+*npReader* (jsa.np.npreader) is a program that extracts Oxford Nanopore
sequencing data from FAST5 files, performs an initial analysis of the date and
streams them to real-time analysis pipelines. These pipelines can run on the
same computer or on computing clouds/high performance clusters.
@@ -34,18 +34,22 @@ Linux distribution software repository, such as::
The library is typically installed to */usr/lib/jni*. Enter this path when
prompted for "Path to HDF library" during installation of Japsa.
+HDF-View (https://www.hdfgroup.org/products/java/release/download.html) also
+contains the neccessary library. Please install HDF-2.10.1 instead of the
+latest version.
+
~~~~~~~~
Synopsis
~~~~~~~~
-*jsa.np.f5reader*: Extract and stream Oxford Nanopore sequencing data in real-time
+*jsa.np.npreader*: Extract and stream Oxford Nanopore sequencing data in real-time. Demultiplexe included.
~~~~~
Usage
~~~~~
::
- jsa.np.f5reader [options]
+ jsa.np.npreader [options]
~~~~~~~
Options
@@ -65,13 +69,17 @@ Options
--format=s Format of sequence reads (fastq or fasta)
(default='fastq')
--minLength=i Minimum read length
- (default='0')
+ (default='1')
--number Add a unique number to read name
(default='false')
--stats Generate a report of read statistics
(default='false')
- --time Extract the sequencing time of each read -- only work with Metrichor > 1.12
+ --time Extract the sequencing time of each read -- experimental
+ (default='false')
+ --exhaustive Whether to traverse the input directory exhaustively (albacore) or lazily (metrichor)
(default='false')
+ --barcode=s The file containing all barcode sequences for demultiplexing.
+ (default='null')
--help Display this usage and exit
(default='false')
@@ -97,16 +105,16 @@ Usage examples
A summary of npReader usage can be obtained by invoking the --help option::
- jsa.np.f5reader --help
+ jsa.np.npreader --help
The simplest way to run *npReader* in GUI mode is by typing::
- jsa.np.f5reader -GUI -realtime
+ jsa.np.npreader -GUI -realtime
and specify various options in the GUI. All of these options can be specified
from the command line::
- jsa.np.f5reader -GUI -realtime -folder c:\Downloads\ -fail -output myrun.fastq --minLength 200 --stats
+ jsa.np.npreader -GUI -realtime -folder c:\Downloads\ -fail -output myrun.fastq --minLength 200 --stats
npReader can run natively on a Windows laptop that runs the Metrichor agent. It
can stream sequence data to multiple analysis pipelines on the same computer
@@ -139,7 +147,7 @@ analysis pipelines, such as::
Once these pipelines are ready, npReader can start streaming data off the
MinION and the Metrichor agent to these pipelines::
- jsa.np.f5reader -realtime -folder c:\Downloads\ -fail -output myrun.fastq \
+ jsa.np.npreader -realtime -folder c:\Downloads\ -fail -output myrun.fastq \
--minLength 200 --streams server1IP:3456,server2IP:3457
One can run *npReader* on a computing cloud if the download folder (containing
@@ -147,22 +155,30 @@ base-called data) can be mounted to the cloud. In such case, npReader can
direct stream data to the pipelines without the need of
*jsa.util.streamServer*::
- jsa.np.f5reader -realtime -folder c:\Downloads\ -fail -output - | \
+ jsa.np.npreader -realtime -folder c:\Downloads\ -fail -output - | \
bwa mem -t 8 -k11 -W20 -r10 -A1 -B1 -O1 -E1 -L0 -Y -K 10000 index - | \
jsa.np.speciesTyping -bam - --index speciesIndex -output output.dat
+*npReader* now supports barcode sequencing demultiplex. For this analysis, it
+requires a FASTA file of barcode tag sequences and will classify output sequences
+based on alignment. User can specify the threshold for alignment confidence from
+the GUI. Demultiplexing results are illustrated as prefix Barcode::|
+added to each output sequence name.
+
+ jsa.np.npreader -GUI -barcode barcode.fasta
+
Japsa also provides *jsa.np.filter*, a tool to bin sequence data in groups of
the user's liking. Like any other streamline tools, jsa.np.filter can run
behind *jsa.util.streamServer* on a remote machine, or can get data directly
from npReader via pipe::
- jsa.np.f5reader -realtime -folder c:\Downloads\ -fail -output - | \
+ jsa.np.npreader -realtime -folder c:\Downloads\ -fail -output - | \
jsa.np.filter -input - -lenMin 2000 --qualMin 10 -output goodreads.fq
One can also use *tee* to group data into different bins *in real-time* with
*jsa.np.filter*::
- jsa.np.f5reader -realtime -folder c:\Downloads\ -fail -output - | \
+ jsa.np.npreader -realtime -folder c:\Downloads\ -fail -output - | \
tee >(jsa.np.filter -input - -lenMax 2000 -output 0k2k.fq) \
>(jsa.np.filter -lenMin 2000 -lenMax 4000 -input - -output 2k4k.fq) \
>(jsa.np.filter -lenMin 4000 -lenMax 6000 -input - -output 4k6k.fq) \
diff --git a/docs/source/tools/jsa.np.npscarf.rst b/docs/source/tools/jsa.np.npscarf.rst
new file mode 100644
index 0000000..9f0f4f7
--- /dev/null
+++ b/docs/source/tools/jsa.np.npscarf.rst
@@ -0,0 +1,187 @@
+----------------------------------------------------------------------------------
+*npScarf*: real-time scaffolder using SPAdes contigs and Nanopore sequencing reads
+----------------------------------------------------------------------------------
+
+*npScarf* (jsa.np.npscarf) is a program that connect contigs from a draft genomes
+to generate sequences that are closer to finish. These pipelines can run on a single laptop
+for microbial datasets. In real-time mode, it can be integrated with simple structural
+analyses such as gene ordering, plasmid forming.
+
+*npScarf* is included in the `Japsa package `_.
+
+~~~~~~~~
+Synopsis
+~~~~~~~~
+
+*jsa.np.npscarf*: Experimental Scaffold and finish assemblies using Oxford Nanopore sequencing reads
+
+~~~~~
+Usage
+~~~~~
+::
+
+ jsa.np.npscarf [options]
+
+~~~~~~~
+Options
+~~~~~~~
+ --seqFile=s Name of the assembly file (sorted by length)
+ (REQUIRED)
+ --input=s Name of the input file, - for stdin
+ (REQUIRED)
+ --format=s format of the input fastq/fasta or sam/bam
+ (default='fastq/fasta')
+ --bwaExe=s Path to bwa
+ (default='bwa')
+ --bwaThread=i Theads used by bwa
+ (default='4')
+ --long Whether report all sequences, including short/repeat contigs (default) or only long/unique/completed sequences.
+ (default='false')
+ --spadesDir=s Name of the output folder by SPAdes: assembly graph and paths will be used for better gap-filling.
+ (default='null')
+ --prefix=s Prefix for the output files
+ (default='out')
+ --genes=s Realtime annotation: name of annotated genes in GFF 3.0 format
+ (default='null')
+ --resistGene=s Realtime annotation: name of antibiotic resistance gene fasta file
+ (default='null')
+ --insertSeq=s Realtime annotation: name of IS fasta file
+ (default='null')
+ --oriRep=s Realtime annotation: name of fasta file containing possible origin of replication
+ (default='null')
+ --minContig=i Minimum contigs length that are used in scaffolding.
+ (default='300')
+ --maxRepeat=i Maximum length of repeat in considering species.
+ (default='7500')
+ --cov=d Expected average coverage of Illumina, <=0 to estimate
+ (default='0.0')
+ --qual=i Minimum quality
+ (default='1')
+ --support=i Minimum supporting long read needed for a link between markers
+ (default='1')
+ --realtime Process in real-time mode. Default is batch mode (false)
+ (default='false')
+ --read=i Minimum number of reads between analyses
+ (default='50')
+ --time=i Minimum number of seconds between analyses
+ (default='10')
+ --verbose Turn on debugging mode
+ (default='false')
+ --help Display this usage and exit
+ (default='false')
+
+
+~~~~~~~~
+See also
+~~~~~~~~
+
+jsa.np.npreader_, jsa.util.streamServer_, jsa.util.streamClient_
+
+.. _jsa.np.npreader: jsa.np.npreader.html
+.. _jsa.util.streamServer: jsa.util.streamServer.html
+.. _jsa.util.streamClient: jsa.util.streamClient.html
+
+
+
+~~~~~~~~~~~~~~
+Usage examples
+~~~~~~~~~~~~~~
+
+A summary of *npScarf* usage can be obtained by invoking the --help option::
+
+ jsa.np.npscarf --help
+
+Input
+=====
+ *npScarf* takes two files as required input::
+
+ jsa.np.npscarf -seq -input
+
+<*draft*> input is the FASTA file containing the pre-assemblies. Normally this
+is the output from running SPAdes on Illumina MiSeq paired end reads.
+
+<*nanopore*> is either the long reads in FASTA/FASTQ file or SAM/BAM formated alignments
+between them to <*draft*> file. We use BWA-MEM as the recommended aligner
+with the fixed parameter set as follow::
+
+ bwa mem -k11 -W20 -r10 -A1 -B1 -O1 -E1 -L0 -a -Y >
+
+The input file format is specified by option --format. The default is FASTA/FASTQ in which
+the path to BWA version 0.7.11 or newer is required. Remember to always *INDEXING* the
+reference before running BWA::
+
+ bwa index
+
+Missing this step would break down the whole pipeline.
+
+Output
+=======
+*npScarf* output is specified by *-prefix* option. The default prefix is \'out\'.
+Normally the tool generate two files: *prefix*.fin.fasta and *prefix*.fin.japsa which
+indicate the result scaffolders in FASTA and JAPSA format.
+
+In realtime mode, if any annotation analysis is enabled, a file named
+*prefix*.anno.japsa is generated instead. This file contains features detected after
+scaffolding.
+
+Real-time scaffolding
+=====================
+To run *npScarf* in streaming mode::
+
+ jsa.np.npscarf -realtime [options]
+
+In this mode, the <*bam*> file will be processed block by block. The size of block
+(number of BAM/SAM records) can be manipulated through option *-read* and *-time*.
+
+The idea of streaming mode is when the input <*nanopore*> file is retrieved in stream.
+npReader is the module that provides such data from fast5 files returned from the real-time
+base-calling cloud service Metrichor. Ones can run::
+
+ jsa.np.npreader -realtime -folder c:\Downloads\ -fail -output - | \
+ jsa.np.npscarf --realtime -bwaExe= -bwaThread=10 -input - -seq > log.out 2>&1
+
+For the same purpose, you can also invoke BWA-MEM explicitly as in the old version of *npScarf*,
+In this case, option --format=SAM must be presented as follow:
+
+ jsa.np.npreader -realtime -folder c:\Downloads\ -fail -output - | \
+ bwa mem -t 10 -k11 -W20 -r10 -A1 -B1 -O1 -E1 -L0 -a -Y -K 3000 - 2> /dev/null | \
+ jsa.np.npscarf --realtime -input - -format=SAM -seq > log.out 2>&1
+
+or if you have the whole set of Nanopore long reads already and want to emulate the
+streaming mode::
+
+ jsa.np.timeEmulate -s 100 -i -output - | \
+ jsa.np.npscarf --realtime -bwaExe= -bwaThread=10 -input - -seq > log.out 2>&1
+
+Note that jsa.np.timeEmulate based on the field *timestamp* located in the read name line to
+decide the order of streaming data. So if your input <*nanopore*> already contains the field,
+you have to sort it::
+
+ jsa.seq.sort -i -o -sortKey=timestamp
+
+or if your file does not have the *timestamp* data yet, you can manually make ones. For example::
+
+ cat | \
+ awk 'BEGIN{time=0.0}NR%4==1{printf "%s timestamp=%.2f\n", $0, time; time++}NR%4!=1{print}' \
+ > <*nanopore-with-time*>
+
+Real-time annotation
+====================
+
+The tool includes usecase for streaming annotation. Ones can provides database of antibiotic
+resistance genes and/or Origin of Replication in FASTA format for the analysis of gene ordering
+and/or plasmid identifying respectively::
+
+ jsa.np.timeEmulate -s 100 -i -output - | \
+ jsa.np.npscarf --realtime -bwaExe= -input - -seq -resistGene -oriRep > log.out 2>&1
+
+Assembly graph
+==============
+
+*npScarf* can read the assembly graph info from SPAdes to make the results more precise.
+The results might be slightly deviate from the old version in term of number of final contigs::
+
+ jsa.np.npscarf --spadesFolder=
+
+where SPAdes_output_directory indicates the result folder of SPAdes, containing files such as contigs.fasta,
+contigs.paths and assembly_graph.fastg.
diff --git a/docs/source/tools/jsa.np.rtMLST.rst b/docs/source/tools/jsa.np.rtMLST.rst
index d2e041a..cabc124 100644
--- a/docs/source/tools/jsa.np.rtMLST.rst
+++ b/docs/source/tools/jsa.np.rtMLST.rst
@@ -42,9 +42,9 @@ Options
See also
~~~~~~~~
-jsa.np.f5reader_, jsa.np.rtSpeciesTyping_, jsa.np.rtStrainTyping_, jsa.np.rtResistGenes_, jsa.util.streamServer_, jsa.util.streamClient_
+jsa.np.npreader_, jsa.np.rtSpeciesTyping_, jsa.np.rtStrainTyping_, jsa.np.rtResistGenes_, jsa.util.streamServer_, jsa.util.streamClient_
-.. _jsa.np.f5reader: jsa.np.f5reader.html
+.. _jsa.np.npreader: jsa.np.npreader.html
.. _jsa.np.rtSpeciesTyping: jsa.np.rtSpeciesTyping.html
.. _jsa.np.rtStrainTyping: jsa.np.rtStrainTyping.html
.. _jsa.np.rtResistGenes: jsa.np.rtResistGenes.html
diff --git a/docs/source/tools/jsa.np.rtResistGenes.rst b/docs/source/tools/jsa.np.rtResistGenes.rst
index 425b561..12ea1c3 100644
--- a/docs/source/tools/jsa.np.rtResistGenes.rst
+++ b/docs/source/tools/jsa.np.rtResistGenes.rst
@@ -51,9 +51,9 @@ Options
See also
~~~~~~~~
-jsa.np.f5reader_, jsa.np.rtSpeciesTyping_, jsa.np.rtStrainTyping_, jsa.util.streamServer_, jsa.util.streamClient_
+jsa.np.npreader_, jsa.np.rtSpeciesTyping_, jsa.np.rtStrainTyping_, jsa.util.streamServer_, jsa.util.streamClient_
-.. _jsa.np.f5reader: jsa.np.f5reader.html
+.. _jsa.np.npreader: jsa.np.npreader.html
.. _jsa.np.rtSpeciesTyping: jsa.np.rtSpeciesTyping.html
.. _jsa.np.rtStrainTyping: jsa.np.rtStrainTyping.html
.. _jsa.util.streamServer: jsa.util.streamServer.html
@@ -65,6 +65,7 @@ jsa.np.f5reader_, jsa.np.rtSpeciesTyping_, jsa.np.rtStrainTyping_, jsa.util.stre
Setting up
~~~~~~~~~~
-Refer to real-time analyais page at https://github.com/mdcao/npAnalysis/
+Refer to the documentation at https://github.com/mdcao/npAnalysis/ for more
+details.
diff --git a/docs/source/tools/jsa.np.rtSpeciesTyping.rst b/docs/source/tools/jsa.np.rtSpeciesTyping.rst
index 27fc192..b8b1725 100644
--- a/docs/source/tools/jsa.np.rtSpeciesTyping.rst
+++ b/docs/source/tools/jsa.np.rtSpeciesTyping.rst
@@ -6,15 +6,11 @@
using Oxford Nanopore sequencing in real-time. It reads data in SAM/BAM format
of the alignments of sequence reads to a collection of species genomes.
-We provide a genome collection of nearly 1500 bacterial species http://genomicsresearch.org/public/researcher/npAnalysis/SpeciesTyping.tar.gz
-(or https://swift.rc.nectar.org.au:8888/v1/AUTH_15574c7fb24c44b3b34069185efba190/npAnalysis/SpeciesTyping.tar.gz).
-Obtain them by::
-
- wget https://swift.rc.nectar.org.au:8888/v1/AUTH_15574c7fb24c44b3b34069185efba190/npAnalysis/SpeciesTyping.tar.gz.tar.gz
- tar zxvf SpeciesTyping.tar.gz
-
-which will generate three folders for the three species.
-
+We provide a genome collection of nearly 1500 bacterial species
+on http://data.genomicsresearch.org/Projects/npAnalysis/.
+Refer to the documentation at https://github.com/mdcao/npAnalysis/ for more
+details.
+
~~~~~~~~
Synopsis
~~~~~~~~
@@ -31,7 +27,7 @@ Usage
~~~~~~~
Options
~~~~~~~
- --output=s Output file
+ --output=s Output file, - for standard output
(default='output.dat')
--bamFile=s The bam file
(REQUIRED)
@@ -45,6 +41,10 @@ Options
(default='50')
--time=i Minimum number of seconds between analyses
(default='30')
+ --web Whether to use Web visualization.
+ (default='false')
+ --log Whether to write mapping details to species2reads.map.
+ (default='false')
--help Display this usage and exit
(default='false')
@@ -53,9 +53,9 @@ Options
See also
~~~~~~~~
-jsa.np.f5reader_, jsa.np.rtStrainTyping_, jsa.np.rtResistGenes_, jsa.util.streamServer_, jsa.util.streamClient_
+jsa.np.npreader_, jsa.np.rtStrainTyping_, jsa.np.rtResistGenes_, jsa.util.streamServer_, jsa.util.streamClient_
-.. _jsa.np.f5reader: jsa.np.f5reader.html
+.. _jsa.np.npreader: jsa.np.npreader.html
.. _jsa.np.rtStrainTyping: jsa.np.rtStrainTyping.html
.. _jsa.np.rtResistGenes: jsa.np.rtResistGenes.html
.. _jsa.util.streamServer: jsa.util.streamServer.html
@@ -87,6 +87,6 @@ to listen on port 3456
and streams data to this pipeline using npReader:
::
- jsa.np.f5reader -GUI -realtime -folder -fail -output data.fastq -stream serverAddress:3456
+ jsa.np.npreader -GUI -realtime -folder -fail -output data.fastq -stream serverAddress:3456
diff --git a/docs/source/tools/jsa.np.rtStrainTyping.rst b/docs/source/tools/jsa.np.rtStrainTyping.rst
index acc678f..b83b0c6 100644
--- a/docs/source/tools/jsa.np.rtStrainTyping.rst
+++ b/docs/source/tools/jsa.np.rtStrainTyping.rst
@@ -9,14 +9,10 @@ gene presence, it makes an inference of the strain, together with the confidence
interval of 95%.
We provide the gene databases for three bacterial species K. pneumoniae,
-E. coli and S. aureus on http://genomicsresearch.org/public/researcher/npAnalysis/StrainTyping.tar.gz
-(or https://swift.rc.nectar.org.au:8888/v1/AUTH_15574c7fb24c44b3b34069185efba190/npAnalysis/StrainTyping.tar.gz).
-Obtain them by::
+E. coli and S. aureus on http://data.genomicsresearch.org/Projects/npAnalysis/.
+Refer to the documentation at https://github.com/mdcao/npAnalysis/ for more
+details.
- wget https://swift.rc.nectar.org.au:8888/v1/AUTH_15574c7fb24c44b3b34069185efba190/npAnalysis/StrainTyping.tar.gz
- tar zxvf StrainTyping.tar.gz
-
-which will generate three folders for the three species.
~~~~~~~~
Synopsis
@@ -56,9 +52,9 @@ Options
See also
~~~~~~~~
-jsa.np.f5reader_, jsa.np.rtSpeciesTyping_, jsa.np.rtResistGenes_, jsa.util.streamServer_, jsa.util.streamClient_
+jsa.np.npreader_, jsa.np.rtSpeciesTyping_, jsa.np.rtResistGenes_, jsa.util.streamServer_, jsa.util.streamClient_
-.. _jsa.np.f5reader: jsa.np.f5reader.html
+.. _jsa.np.npreader: jsa.np.npreader.html
.. _jsa.np.rtSpeciesTyping: jsa.np.rtSpeciesTyping.html
.. _jsa.np.rtResistGenes: jsa.np.rtResistGenes.html
.. _jsa.util.streamServer: jsa.util.streamServer.html
@@ -88,6 +84,6 @@ to listen on port 3457
and streams data to this pipeline using npReader:
::
- jsa.np.f5reader -GUI -realtime -folder -fail -output data.fastq -stream serverAddress:3457
+ jsa.np.npreader -GUI -realtime -folder -fail -output data.fastq -stream serverAddress:3457
diff --git a/docs/source/tools/jsa.seq.emalign.rst b/docs/source/tools/jsa.seq.emalign.rst
index 8982a7b..80ad1d3 100644
--- a/docs/source/tools/jsa.seq.emalign.rst
+++ b/docs/source/tools/jsa.seq.emalign.rst
@@ -19,6 +19,8 @@ Usage
~~~~~~~
Options
~~~~~~~
+ --iteration=i Number of iteration
+ (default='5')
--help Display this usage and exit
(default='false')
diff --git a/docs/source/tools/jsa.seq.gff2fasta.rst b/docs/source/tools/jsa.seq.gff2fasta.rst
index 47a5eb7..7fbbefa 100644
--- a/docs/source/tools/jsa.seq.gff2fasta.rst
+++ b/docs/source/tools/jsa.seq.gff2fasta.rst
@@ -1,8 +1,8 @@
-------------------------------------------
- *jsa.seq.gff2fasta*: Extract gene sequences
+*jsa.seq.gff2fasta*: Extract gene sequences
-------------------------------------------
- *jsa.seq.gff2fasta* extract the functional sequences (genes, CDS, etc) from
+*jsa.seq.gff2fasta* extract the functional sequences (genes, CDS, etc) from
a gff file and a sequence file.
~~~~~~~~
diff --git a/docs/source/tools/jsa.sim.capsim.rst b/docs/source/tools/jsa.sim.capsim.rst
new file mode 100644
index 0000000..dace12e
--- /dev/null
+++ b/docs/source/tools/jsa.sim.capsim.rst
@@ -0,0 +1,57 @@
+----------------------------------------------------------------------------
+*capsim*: Simulating the Dynamics of Targeted Capture Sequencing with CapSim
+----------------------------------------------------------------------------
+
+*capsim* (jsa.sim.capsim) is a tool to simulate target capture sequencing. Its
+simulates the dynamics of capture process
+
+~~~~~~~~
+Synopsis
+~~~~~~~~
+
+*jsa.sim.capsim*: Simulate capture sequencing
+
+~~~~~
+Usage
+~~~~~
+::
+
+ jsa.sim.capsim [options]
+
+~~~~~~~
+Options
+~~~~~~~
+ --reference=s Name of genome to be
+ (REQUIRED)
+ --probe=s File containing probes mapped to the reference in bam format
+ (default='null')
+ --logFile=s Log file
+ (default='-')
+ --ID=s A unique ID for the data set
+ (default='')
+ --miseq=s Name of read file if miseq is simulated
+ (default='null')
+ --pacbio=s Name of read file if pacbio is simulated
+ (default='null')
+ --fmedian=i Median of fragment size at shearing
+ (default='2000')
+ --fshape=d Shape parameter of the fragment size distribution
+ (default='6.0')
+ --num=i Number of fragments
+ (default='1000000')
+ --pblen=i PacBio: Average (polymerase) read length
+ (default='30000')
+ --illen=i Illumina: read length
+ (default='300')
+ --seed=i Random seed, 0 for a random seed
+ (default='0')
+ --help Display this usage and exit
+ (default='false')
+
+
+
+
+~~~~~~~~~~~~~
+Usage samples
+~~~~~~~~~~~~~
+
diff --git a/docs/source/tools/jsa.util.streamClient.rst b/docs/source/tools/jsa.util.streamClient.rst
index 835aa16..19adc2c 100644
--- a/docs/source/tools/jsa.util.streamClient.rst
+++ b/docs/source/tools/jsa.util.streamClient.rst
@@ -33,11 +33,11 @@ Options
See also
~~~~~~~~
-jsa.util.streamServer_, jsa.np.filter_, jsa.np.f5reader_
+jsa.util.streamServer_, jsa.np.filter_, jsa.np.npreader_
.. _jsa.util.streamServer: jsa.util.streamServer.html
.. _jsa.np.filter: jsa.np.filter.html
-.. _jsa.np.f5reader: jsa.np.f5reader.html
+.. _jsa.np.npreader: jsa.np.npreader.html
diff --git a/docs/source/tools/jsa.util.streamServer.rst b/docs/source/tools/jsa.util.streamServer.rst
index e0a2e3e..75a9a04 100644
--- a/docs/source/tools/jsa.util.streamServer.rst
+++ b/docs/source/tools/jsa.util.streamServer.rst
@@ -34,11 +34,11 @@ Options
See also
~~~~~~~~
-jsa.util.streamClient_, jsa.np.filter_, jsa.np.f5reader_
+jsa.util.streamClient_, jsa.np.filter_, jsa.np.npreader_
.. _jsa.util.streamClient: jsa.util.streamClient.html
.. _jsa.np.filter: jsa.np.filter.html
-.. _jsa.np.f5reader: jsa.np.f5reader.html
+.. _jsa.np.npreader: jsa.np.npreader.html
diff --git a/docs/source/tools/jsa.xm.compress.rst b/docs/source/tools/jsa.xm.compress.rst
new file mode 100644
index 0000000..0ba2a67
--- /dev/null
+++ b/docs/source/tools/jsa.xm.compress.rst
@@ -0,0 +1,75 @@
+---------------------------------------------------------
+*Expert Model*: tool for compression of genomic sequences
+---------------------------------------------------------
+
+*jsa.xm.compress* in the implementation of the expert model (XM) algorithm for
+compression of genomics sequences. The source code is included in the
+`Japsa package `_.
+Please see check the installation_ page for instructions.
+
+.. _installation: ../install.html
+
+~~~~~~~~
+Synopsis
+~~~~~~~~
+
+*jsa.xm.compress*: Compression of DNA/protein sequences
+
+~~~~~
+Usage
+~~~~~
+::
+
+ jsa.xm.compress [options] file1 file2 ...
+
+~~~~~~~
+Options
+~~~~~~~
+ --hashSize=i Hash size
+ (default='11')
+ --context=i Length of the context
+ (default='15')
+ --limit=i Expert Limit
+ (default='200')
+ --threshold=d Listen threshold
+ (default='0.15')
+ --chance=i Chances
+ (default='20')
+ --binaryHash Use binary hash or not
+ (default='false')
+ --offsetType=s Way of update offset/palindrome expert: possible value count, subs
+ (default='counts')
+ --real=s File name of the real compression
+ (default='null')
+ --decode=s File name of the encoded
+ (default='null')
+ --output=s The output file of decoded file
+ (default='decoded')
+ --info=s File name of the infomation content
+ (default='null')
+ --markov=s File name of the markov infomation content
+ (default='null')
+ --optimise Running in optimise mode, just report the entropy,recommended for long sequence
+ (default='false')
+ --checkPoint=i Frequency of check point
+ (default='1000000')
+ --hashType=s Type of Hash table: hash=hashtable, sft=SuffixTree,sfa = SuffixArray
+ (default='hash')
+ --selfRep Propose experts from the sequence to compressed?
+ (default='true')
+ --help Display this usage and exit
+ (default='false')
+
+
+
+
+~~~~~~~~
+Citation
+~~~~~~~~
+
+If you find XM useful for your research, please cite
+
+Cao MD, Dix TI, Allison L, and Mears C,
+*A simple statistical algorithm for biological sequence compression*,
+Data Compression Conference, 2007 (DCC'07), Snowbird, UT, pp43-52.
+
diff --git a/libs/gs-algo-1.3.jar b/libs/gs-algo-1.3.jar
new file mode 100644
index 0000000..1e720ee
Binary files /dev/null and b/libs/gs-algo-1.3.jar differ
diff --git a/libs/gs-core-1.3.jar b/libs/gs-core-1.3.jar
new file mode 100644
index 0000000..c8f50fd
Binary files /dev/null and b/libs/gs-core-1.3.jar differ
diff --git a/libs/gs-ui-1.3.jar b/libs/gs-ui-1.3.jar
new file mode 100644
index 0000000..4661b87
Binary files /dev/null and b/libs/gs-ui-1.3.jar differ
diff --git a/libs/htsjdk-1.126.jar b/libs/htsjdk-1.126.jar
deleted file mode 100644
index cf0c02f..0000000
Binary files a/libs/htsjdk-1.126.jar and /dev/null differ
diff --git a/libs/htsjdk-2.10.1.jar b/libs/htsjdk-2.10.1.jar
new file mode 100644
index 0000000..399634d
Binary files /dev/null and b/libs/htsjdk-2.10.1.jar differ
diff --git a/libs/javax.json-1.0.4.jar b/libs/javax.json-1.0.4.jar
new file mode 100644
index 0000000..09967d8
Binary files /dev/null and b/libs/javax.json-1.0.4.jar differ
diff --git a/libs/javax.json-api-1.0.jar b/libs/javax.json-api-1.0.jar
new file mode 100644
index 0000000..d276c79
Binary files /dev/null and b/libs/javax.json-api-1.0.jar differ
diff --git a/libs/jfxrt.jar b/libs/jfxrt.jar
new file mode 100755
index 0000000..55858aa
Binary files /dev/null and b/libs/jfxrt.jar differ
diff --git a/libs/junit-4.11.jar b/libs/junit-4.11.jar
deleted file mode 100644
index 09db3e7..0000000
Binary files a/libs/junit-4.11.jar and /dev/null differ
diff --git a/libs/slf4j-api-1.7.25.jar b/libs/slf4j-api-1.7.25.jar
new file mode 100644
index 0000000..7e62f13
Binary files /dev/null and b/libs/slf4j-api-1.7.25.jar differ
diff --git a/libs/slf4j-simple-1.7.25.jar b/libs/slf4j-simple-1.7.25.jar
new file mode 100644
index 0000000..b29ca12
Binary files /dev/null and b/libs/slf4j-simple-1.7.25.jar differ
diff --git a/makefile.dev b/makefile.dev
new file mode 100644
index 0000000..efb729c
--- /dev/null
+++ b/makefile.dev
@@ -0,0 +1,180 @@
+#############################################################################
+# Copyright (c) Minh Duc Cao, Monash Uni & UQ, All rights reserved. #
+# #
+# Redistribution and use in source and binary forms, with or without #
+# modification, are permitted provided that the following conditions #
+# are met: #
+# #
+# 1. Redistributions of source code must retain the above copyright notice, #
+# this list of conditions and the following disclaimer. #
+# 2. Redistributions in binary form must reproduce the above copyright #
+# notice, this list of conditions and the following disclaimer in the #
+# documentation and/or other materials provided with the distribution. #
+# 3. Neither the names of the institutions nor the names of the contributors#
+# may be used to endorse or promote products derived from this software #
+# without specific prior written permission. #
+# #
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS #
+# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, #
+# THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR #
+# PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR #
+# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, #
+# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, #
+# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR #
+# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF #
+# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING #
+# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS #
+# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #
+############################################################################/
+
+
+############################################################################
+# A generic makefile for java
+#
+############################################################################
+
+
+###########################################################################
+# Modify the following parameter to suit your setting
+
+#Class directory
+CLASS_DIR=target/dev-classes
+
+#Source directory
+SRC_MAIN_DIR=src/main/java
+SRC_DEV_DIR=src/dev/java
+##SRC_WORK_DIR=src/work/java
+
+#Java compiler
+JCC=javac
+
+#target jar file
+JAR_FILE=japsa-dev.jar
+
+#External library directory
+LIB_DIR=libs
+
+#List of external libraries
+
+#EXT_LIBS=colloquial.jar commons-math3-3.0.jar jhdf5.jar jhdfobj.jar htsjdk-2.10.1.jar guava-18.0.jar jcommon-1.0.23.jar jfreechart-1.0.19.jar JRIEngine.jar JRI.jar gs-core-1.3.jar gs-ui-1.3.jar gs-algo-1.3.jar javax.json-api-1.0.jar jfxrt.jar
+EXT_LIBS=colloquial.jar commons-math3-3.0.jar jhdf5.jar jhdfobj.jar htsjdk-2.10.1.jar guava-18.0.jar jcommon-1.0.23.jar jfreechart-1.0.19.jar JRIEngine.jar JRI.jar gs-core-1.3.jar gs-ui-1.3.jar gs-algo-1.3.jar javax.json-api-1.0.jar javax.json-1.0.4.jar slf4j-api-1.7.25.jar slf4j-simple-1.7.25.jar jfxrt.jar
+
+###########################################################################
+##What this scripts does:
+# 1. Find all packages by searching SRC_DIR subdirectories
+# 2. Get all the java file, and compile them to the CLASS_DIR directory
+# 5. Clean back up files (*.bak, ~, etc)
+# To be implemented:
+# 3. Copy all the resources to class directory
+# 4. Make java file
+##########################################################################
+
+COMMA:= ,
+EMPTY:=
+SPACE:= $(SPACE) $(SPACE)
+LIBS:=$(subst $(SPACE),:, $(addprefix $(LIB_DIR)/, $(EXT_LIBS)))
+
+##Get all the packages in $(SRC_MAIN_DIR).
+PACKAGE_MAIN_DIRS := $(shell echo `cd $(SRC_MAIN_DIR);find . -type d`)
+SRC_MAIN_DIRS := $(addprefix $(SRC_MAIN_DIR)/, $(PACKAGE_MAIN_DIRS))
+CLASS_MAIN_DIRS := $(addprefix $(CLASS_DIR)/, $(PACKAGE_MAIN_DIRS))
+##Get the list of java files
+JAVA_MAIN_FILES := $(foreach dir,$(SRC_MAIN_DIRS),$(wildcard $(dir)/*.java))
+CLASS_MAIN_FILES := $(subst $(SRC_MAIN_DIR)/,$(CLASS_DIR)/,$(JAVA_MAIN_FILES:.java=.class))
+
+PACKAGE_DEV_DIRS := $(shell echo `cd $(SRC_DEV_DIR);find . -type d`)
+SRC_DEV_DIRS := $(addprefix $(SRC_DEV_DIR)/, $(PACKAGE_DEV_DIRS))
+CLASS_DEV_DIRS := $(addprefix $(CLASS_DIR)/, $(PACKAGE_DEV_DIRS))
+##Get the list of java files
+JAVA_DEV_FILES := $(foreach dir,$(SRC_DEV_DIRS),$(wildcard $(dir)/*.java))
+CLASS_DEV_FILES := $(subst $(SRC_DEV_DIR)/,$(CLASS_DIR)/,$(JAVA_DEV_FILES:.java=.class))
+
+#TODO: merge three src dirs to one by foreach
+
+#####################Make targets
+
+VPATH=$(subst ' ',':',$(PACKAGE_MAIN_DIRS))
+$(CLASS_DIR)/%.class: $(SRC_MAIN_DIR)/%.java $(CLASS_DIR)
+ $(JCC) -sourcepath $(SRC_MAIN_DIR) -cp $(CLASS_DIR):$(LIBS) -nowarn -d $(CLASS_DIR) $(JDEBUGFLAGS) $<
+
+VPATH=$(subst ' ',':',$(PACKAGE_DEV_DIRS))
+$(CLASS_DIR)/%.class: $(SRC_DEV_DIR)/%.java $(CLASS_DIR)
+ $(JCC) -sourcepath $(SRC_DEV_DIR):$(SRC_MAIN_DIR) -cp $(CLASS_DIR):$(LIBS) -nowarn -d $(CLASS_DIR) $(JDEBUGFLAGS) $<
+
+all: jar
+
+main-classes: $(CLASS_MAIN_FILES)
+dev-classes: $(CLASS_DEV_FILES)
+
+###Create the class directory if neccessary
+$(CLASS_DIR):
+ mkdir -p $(CLASS_DIR)
+
+images:
+ cp -r $(SRC_MAIN_DIR)/japsa/bio/misc/dnaPlatform/gui/images $(CLASS_DIR)/japsa/bio/misc/dnaPlatform/gui/
+ cp -r $(SRC_MAIN_DIR)/japsa/seq/nanopore/icons $(CLASS_DIR)/japsa/seq/nanopore/
+
+$(JAR_FILE): main-classes dev-classes images
+ jar cf $(JAR_FILE) -C $(CLASS_DIR) .
+
+jar: $(JAR_FILE)
+
+clean:
+ @@for i in $(PACKAGE_MAIN_DIRS); do \
+ echo "rm -f $(SRC_MAIN_DIR)/$$i/*.bak $(SRC_MAIN_DIR)/$$i/*.class $(SRC_MAIN_DIR)/$$i/*~ $(CLASS_DIR)/$$i/*.class"; \
+ rm -f $(SRC_MAIN_DIR)/$$i/*.bak $(SRC_MAIN_DIR)/$$i/*.class $(SRC_MAIN_DIR)/$$i/*~ $(CLASS_DIR)/$$i/*.class; \
+ done
+ @@for i in $(PACKAGE_DEV_DIRS); do \
+ echo "rm -f $(SRC_DEV_DIR)/$$i/*.bak $(SRC_DEV_DIR)/$$i/*.class $(SRC_DEV_DIR)/$$i/*~ $(CLASS_DIR)/$$i/*.class"; \
+ rm -f $(SRC_DEV_DIR)/$$i/*.bak $(SRC_DEV_DIR)/$$i/*.class $(SRC_DEV_DIR)/$$i/*~ $(CLASS_DIR)/$$i/*.class; \
+ done
+ rm -f install.sh install.bat
+
+##############################################################################
+# INSTALLATION
+##############################################################################
+
+ifdef INSTALL_DIR
+O_INS_DIR=-installDir=\"${INSTALL_DIR}\"
+else
+O_INS_DIR=
+endif
+
+
+ifdef JLP
+O_JLP=-jlp=\"${JLP}\"
+else
+O_JLP=
+endif
+
+ifdef MXMEM
+O_MXMEM=-xmx=${MXMEM}
+else
+O_MXMEM=
+endif
+
+
+ifdef SERVER
+O_SERVER=--server=${SERVER}
+else
+O_SERVER=
+endif
+
+
+##=/usr/lib/jni:/usr/lib/R/site-library/rJava/jri
+
+pre-install: jar
+ @@echo "java -cp $(JAR_FILE):$(LIB_DIR)/guava-18.0.jar:$(LIB_DIR)/slf4j-api-1.7.25.jar:$(LIB_DIR)/slf4j-simple-1.7.25.jar japsadev.util.deploy.DevDeploy --mode install --libs $(subst $(SPACE),:, $(EXT_LIBS)) $(O_INS_DIR) $(O_JLP) $(O_MXMEM) $(O_SERVER) --compiler \"`$(JCC) -version 2>&1`\"" > install.sh && \
+ chmod u+x install.sh && \
+ echo "java -cp $(JAR_FILE):$(LIB_DIR)\guava-18.0.jar:$(LIB_DIR)/slf4j-api-1.7.25.jar:$(LIB_DIR)/slf4j-simple-1.7.25.jar japsadev.util.deploy.DevDeploy --mode install --libs $(subst $(SPACE),:, $(EXT_LIBS)) $(O_INS_DIR) $(O_JLP) $(O_MXMEM) $(O_SERVER)--compiler \"`$(JCC) -version 2>&1`\"" > install.bat &&\
+ echo "Installation scripts created"
+
+install: pre-install
+ ./install.sh
+
+uninstall:
+ @@java -cp $(JAR_FILE):$(LIB_DIR)/guava-18.0.jar:$(LIB_DIR)/slf4j-api-1.7.25.jar:$(LIB_DIR)/slf4j-simple-1.7.25.jar japsadev.util.deploy.DevDeploy --mode uninstall --libs $(subst $(SPACE),:, $(EXT_LIBS)) ${O_INS_DIR} && echo "Japsa uninstalled!"
+
+galaxy: jar
+ @@java -cp $(JAR_FILE):$(LIB_DIR)/guava-18.0.jar:$(LIB_DIR)/slf4j-api-1.7.25.jar:$(LIB_DIR)/slf4j-simple-1.7.25.jar japsadev.util.deploy.DevDeploy --mode galaxy && echo "Galaxy wrappers generated!"
+
diff --git a/src/dev/java/japsadev/bio/go/GoTerm.java b/src/dev/java/japsadev/bio/go/GoTerm.java
new file mode 100644
index 0000000..2951407
--- /dev/null
+++ b/src/dev/java/japsadev/bio/go/GoTerm.java
@@ -0,0 +1,80 @@
+/*****************************************************************************
+ * Copyright (c) Minh Duc Cao, Monash Uni & UQ, All rights reserved. *
+ * *
+ * Redistribution and use in source and binary forms, with or without *
+ * modification, are permitted provided that the following conditions *
+ * are met: *
+ * *
+ * 1. Redistributions of source code must retain the above copyright notice, *
+ * this list of conditions and the following disclaimer. *
+ * 2. Redistributions in binary form must reproduce the above copyright *
+ * notice, this list of conditions and the following disclaimer in the *
+ * documentation and/or other materials provided with the distribution. *
+ * 3. Neither the names of the institutions nor the names of the contributors*
+ * may be used to endorse or promote products derived from this software *
+ * without specific prior written permission. *
+ * *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS *
+ * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, *
+ * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR *
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR *
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, *
+ * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, *
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR *
+ * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF *
+ * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING *
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
+ ****************************************************************************/
+
+/*****************************************************************************
+ * Revision History
+ * 20 Jul 2015 - Minh Duc Cao: Created
+ *
+ ****************************************************************************/
+package japsadev.bio.go;
+
+
+
+import java.util.ArrayList;
+
+/**
+ * A simple GO Term class, will be replaced by a standard library
+ * @author minhduc
+ *
+ */
+public class GoTerm {
+ String ID;
+ String name;
+ String desc;
+
+ ArrayList relationship;
+
+ public GoTerm (String anID){
+ ID = anID;
+ relationship = new ArrayList();
+ }
+
+ public void addRelationShip(TypeRelationship type, GoTerm term){
+ relationship.add(new GoRelationship(type, term));
+ }
+
+ public static class GoRelationship{
+ GoTerm relTerm;
+ TypeRelationship relType;
+ GoRelationship(TypeRelationship type, GoTerm term){
+ relType = type;
+ relTerm = term;
+ }
+ }
+
+ /**
+ * @param args
+ */
+ public static void main(String[] args) {
+
+
+
+ }
+
+}
diff --git a/src/dev/java/japsadev/bio/go/ODOParser.java b/src/dev/java/japsadev/bio/go/ODOParser.java
new file mode 100644
index 0000000..88e2833
--- /dev/null
+++ b/src/dev/java/japsadev/bio/go/ODOParser.java
@@ -0,0 +1,590 @@
+/*****************************************************************************
+ * Copyright (c) Minh Duc Cao, Monash Uni & UQ, All rights reserved. *
+ * *
+ * Redistribution and use in source and binary forms, with or without *
+ * modification, are permitted provided that the following conditions *
+ * are met: *
+ * *
+ * 1. Redistributions of source code must retain the above copyright notice, *
+ * this list of conditions and the following disclaimer. *
+ * 2. Redistributions in binary form must reproduce the above copyright *
+ * notice, this list of conditions and the following disclaimer in the *
+ * documentation and/or other materials provided with the distribution. *
+ * 3. Neither the names of the institutions nor the names of the contributors*
+ * may be used to endorse or promote products derived from this software *
+ * without specific prior written permission. *
+ * *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS *
+ * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, *
+ * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR *
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR *
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, *
+ * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, *
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR *
+ * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF *
+ * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING *
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
+ ****************************************************************************/
+
+/*****************************************************************************
+ * Revision History
+ * 20 Jul 2015 - Minh Duc Cao: Created
+ *
+ ****************************************************************************/
+package japsadev.bio.go;
+
+import japsa.seq.Sequence;
+import japsa.seq.SequenceReader;
+import japsadev.bio.go.GoTerm.GoRelationship;
+import org.slf4j.Logger;
+import org.slf4j.LoggerFactory;
+
+import java.io.BufferedReader;
+import java.io.IOException;
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.HashSet;
+
+/**
+ * @author minhduc
+ *
+ */
+public class ODOParser {
+ private static final Logger LOG = LoggerFactory.getLogger(ODOParser.class);
+
+ HashMap terms = new HashMap();
+
+ HashMap relTypes =
+ new HashMap();
+
+ public ODOParser(String fileName) throws IOException{
+
+ //First round to read
+ BufferedReader bf = SequenceReader.openFile(fileName);
+ String line = "";
+ int lineNo = 0;
+
+ while ((line = bf.readLine())!=null) {
+ lineNo ++;
+ line = line.trim();
+ if (line.equals("[Term]")){
+ line = bf.readLine().trim();
+ lineNo ++;
+ if (!line.startsWith("id: "))
+ throw new RuntimeException("Wrong format at line " + lineNo);
+
+ String ID = line.substring(4);
+ line = bf.readLine().trim();
+ lineNo ++;
+
+ if (!line.startsWith("name: "))
+ throw new RuntimeException("Wrong format at line " + lineNo);
+ String name = line.substring(5).trim();
+
+ GoTerm currentTerm = new GoTerm(ID);
+ currentTerm.name = name;
+ terms.put(ID, currentTerm);
+ continue;
+ }
+
+ if (line.equals("[Typedef]")){
+ line = bf.readLine().trim();
+ lineNo ++;
+ if (!line.startsWith("id: "))
+ throw new RuntimeException("Wrong format at line " + lineNo);
+
+ String ID = line.substring(4);
+ line = bf.readLine().trim();
+ lineNo ++;
+
+ if (!line.startsWith("name: "))
+ throw new RuntimeException("Wrong format at line " + lineNo);
+ String name = line.substring(5);
+
+ TypeRelationship relationship = new TypeRelationship (ID, name);
+ relTypes.put(ID, relationship);
+ continue;
+ }
+ }
+ bf.close();
+
+
+ //System.out.println(terms.size());
+ //System.out.println(relTypes.size());
+
+
+ //Second round to build relationship
+ bf = SequenceReader.openFile(fileName);
+ line = "";
+ lineNo = 0;
+
+ while ((line = bf.readLine())!=null) {
+ lineNo ++;
+ line = line.trim();
+
+
+ if (line.equals("[Term]")){
+ line = bf.readLine().trim();
+ lineNo ++;
+ if (!line.startsWith("id: "))
+ throw new RuntimeException("Wrong format at line " + lineNo);
+
+ String ID = line.substring(4);
+
+ line = bf.readLine().trim();
+ lineNo ++;
+
+ GoTerm currentTerm = terms.get(ID);
+
+ // System.out.println(currentTerm + " " + lineNo);
+ while(true){
+ line = bf.readLine();
+ lineNo ++;
+ if (line == null)
+ break;
+
+ line = line.trim();
+ if (line.length() == 0)
+ break;
+
+
+ if (line.startsWith("is_a: ")){
+ String[] toks = line.split(" ",4);
+ TypeRelationship relationship = relTypes.get("is_a");
+ GoTerm term = terms.get(toks[1]);
+ if (term != null)
+ currentTerm.addRelationShip(relationship, term);
+ else
+ LOG.warn("Term " + toks[1] + " not found at line " + lineNo);
+
+ }
+ if (line.startsWith("relationship: ")){
+ String[] toks = line.split(" ",4);
+ TypeRelationship relationship = relTypes.get(toks[1]);
+ GoTerm term = terms.get(toks[2]);
+
+ if (term != null)
+ currentTerm.addRelationShip(relationship, term);
+ else
+ LOG.warn("Term " + toks[2] + " not found at line " + lineNo);
+
+ }//if
+ }//while
+ }//if
+ }//while
+ bf.close();
+ }
+
+ static HashSet conferDrugList = new HashSet();
+ static HashSet conferList = new HashSet();
+ static StringBuilder sb = new StringBuilder();
+
+
+ static void confer(GoTerm term){
+ for (GoRelationship rel:term.relationship){
+ if (rel.relType.ID.equals("confers_resistance_to"))
+ conferList.add(rel.relTerm.name);
+ else if (rel.relType.ID.equals("confers_resistance_to_drug"))
+ conferDrugList.add(rel.relTerm.name);
+ else if (rel.relType.ID.equals("is_a"))
+ confer(rel.relTerm);
+ }
+ }
+
+ static void list(GoTerm term, HashSet relationList){
+ for (GoRelationship rel:term.relationship){
+ if (relationList.contains(rel.relType.ID) && (!blackList.contains(rel.relTerm.ID))){
+ sb.append(term.ID + "->" + rel.relTerm.ID + "[" + rel.relTerm.name + "];");
+ list(rel.relTerm, relationList);
+ }
+ }
+ }
+
+ static ODOParser odo;
+ //static HashSet relationList = new HashSet();
+
+ static HashSet blackList = new HashSet();
+ static{
+ blackList.add("AL513383.1.gene203");
+ blackList.add("DQ303918.1.gene1");
+ blackList.add("AL513383.1.gene214");
+
+ //CONFUSED GENE
+ //blackList.add("KM998962.1.gene1");
+ //
+
+ //DQ464881.1.gene2
+ //AL513383.1.gene203
+ //AY566250.1.gene1
+ //AB091338.1.gene1
+
+
+ //blackList.add("ARO:1000001");//every thing
+ //blackList.add("ARO:3000557");// [antibiotic inactivation enzyme]
+ //blackList.add("ARO:3000000");//determinant of antibiotic resistance
+ /******************************************************************
+ blackList.add("ARO:3000722");
+ blackList.add("ARO:3000717");
+ blackList.add("ARO:3002639");
+ blackList.add("ARO:3000744"); //-same as above
+ blackList.add("ARO:3000730");
+ blackList.add("ARO:3002701");
+ blackList.add("ARO:3002703");
+ blackList.add("ARO:3000816");
+ blackList.add("ARO:3000723");
+ blackList.add("ARO:3000498");
+ blackList.add("ARO:3000744");
+ blackList.add("ARO:3000795");
+ blackList.add("ARO:3000796");
+ blackList.add("ARO:3002675");
+ blackList.add("ARO:3003079");
+ blackList.add("ARO:3000730");
+ blackList.add("ARO:3003078");
+ blackList.add("ARO:3000743");
+ blackList.add("ARO:3000822");
+ blackList.add("ARO:3002986");
+ blackList.add("ARO:3000730");
+ blackList.add("ARO:3000717");
+ blackList.add("ARO:3000735");
+ blackList.add("ARO:3000722");
+ blackList.add("ARO:3000734");
+ blackList.add("ARO:3000730");
+ blackList.add("ARO:3003577");
+ blackList.add("ARO:3000299");
+ blackList.add("ARO:3000717");
+ blackList.add("ARO:3000027");
+ blackList.add("ARO:3000730");
+ blackList.add("ARO:3000734");
+ blackList.add("ARO:3000074");
+ blackList.add("ARO:3000838");
+ blackList.add("ARO:3000735");
+ blackList.add("ARO:3000734");
+ blackList.add("ARO:3003301");
+ blackList.add("ARO:3000723");
+ blackList.add("ARO:3000722");
+ blackList.add("ARO:3000735");
+ blackList.add("ARO:3000733");
+ blackList.add("ARO:3003063");
+ blackList.add("ARO:3000781");
+ blackList.add("ARO:3000780");
+ blackList.add("ARO:3000768");
+ blackList.add("ARO:3000782");
+ blackList.add("ARO:3000722");
+ blackList.add("ARO:3003392");
+ blackList.add("ARO:3000730");
+ blackList.add("ARO:3001214");
+ blackList.add("ARO:3000717");
+ blackList.add("ARO:3000737");
+ blackList.add("ARO:3003064");
+ blackList.add("ARO:3003052");
+ blackList.add("ARO:3003053");
+ blackList.add("ARO:3003051");
+ blackList.add("ARO:3000502");
+ blackList.add("ARO:3000502");
+ blackList.add("ARO:3000502");
+ blackList.add("ARO:3000499");
+ blackList.add("ARO:3000499");
+ blackList.add("ARO:3000499");
+ blackList.add("ARO:3000810");
+ blackList.add("ARO:3000811");
+ blackList.add("ARO:3003046");
+ blackList.add("ARO:3003047");
+ blackList.add("ARO:3000804");
+ blackList.add("ARO:3000805");
+ blackList.add("ARO:3000785");
+ blackList.add("ARO:3000254");
+ blackList.add("ARO:3000254");
+ blackList.add("ARO:3000254");
+ blackList.add("ARO:3000254");
+ blackList.add("ARO:3000237");
+ blackList.add("ARO:3000379");
+ blackList.add("ARO:3000379");
+ blackList.add("ARO:3000379");
+ blackList.add("ARO:3002983");
+ blackList.add("ARO:3002983");
+ blackList.add("ARO:3003033");
+ blackList.add("ARO:3002982");
+ blackList.add("ARO:3002982");
+ blackList.add("ARO:3003034");
+ blackList.add("ARO:3000774");
+ blackList.add("ARO:3000774");
+ blackList.add("ARO:3000812");
+ blackList.add("ARO:3000807");
+ blackList.add("ARO:3000806");
+ blackList.add("ARO:3000377");
+ blackList.add("ARO:3000377");
+ blackList.add("ARO:3000794");
+ blackList.add("ARO:3000794");
+ blackList.add("ARO:3000794");
+ blackList.add("ARO:3000792");
+ blackList.add("ARO:3000792");
+ blackList.add("ARO:3000792");
+ blackList.add("ARO:3000793");
+ blackList.add("ARO:3000793");
+ blackList.add("ARO:3000793");
+ blackList.add("ARO:3000800");
+ blackList.add("ARO:3000800");
+ blackList.add("ARO:3000801");
+ blackList.add("ARO:3000801");
+ blackList.add("ARO:3000790");
+ blackList.add("ARO:3000789");
+ blackList.add("ARO:3000791");
+ blackList.add("ARO:3000779");
+ blackList.add("ARO:3000779");
+ blackList.add("ARO:3000779");
+ blackList.add("ARO:3000207");
+ blackList.add("ARO:3000216");
+ blackList.add("ARO:3000809");
+ blackList.add("ARO:3000808");
+ blackList.add("ARO:3000533");
+ blackList.add("ARO:3000378");
+ blackList.add("ARO:3000803");
+ blackList.add("ARO:3000777");
+ blackList.add("ARO:3000778");
+ blackList.add("ARO:3000535");
+ blackList.add("ARO:3000775");
+ blackList.add("ARO:3000775");
+ blackList.add("ARO:3003057");
+ blackList.add("ARO:3003056");
+ blackList.add("ARO:3003055");
+ blackList.add("ARO:3003039");
+ blackList.add("ARO:3000802");
+ blackList.add("ARO:3000776");
+ blackList.add("ARO:3000783");
+ blackList.add("ARO:3000784");
+ blackList.add("ARO:3003010");
+ blackList.add("ARO:3003009");
+ blackList.add("ARO:3000206");
+/******************************************************************/
+ /*
+ * ARO:3000839
+ * ARO:3000725
+ * ARO:3000815
+ * ARO:3000765
+ * ARO:3000764
+ * ARO:3000508
+ * ARO:3003066
+ * ARO:3003067
+ * ARO:3000725
+ * ARO:3000764
+ * ARO:3000656
+ * ARO:3000502
+ * ARO:3000499
+ * ARO:3000730
+ *
+ *
+ * KM998962.1.gene1
+ *
+ *
+ */
+ }
+
+ public static void resistanceClass(GoTerm term, HashSet allClass, HashSet myClass) throws IOException{
+ if (allClass.contains(term.ID))
+ myClass.add(term.name);
+ else{
+ for (GoRelationship rel:term.relationship){
+ if (rel.relType.ID.equals("is_a"))
+ resistanceClass(rel.relTerm, allClass, myClass);
+ }
+ }
+
+ }
+
+ public static HashMap> checkGenes(ArrayList seqs) throws IOException{
+ //Set up known classes
+ HashSet abrGroups = new HashSet();
+ abrGroups.add("ARO:3000052");// phenicol resistance gene
+ abrGroups.add("ARO:3000102");// fluoroquinolone resistance gene
+ abrGroups.add("ARO:3000104");// aminoglycoside resistance gene
+ abrGroups.add("ARO:3000129");// beta-lactam resistance gene
+ abrGroups.add("ARO:3000240");// streptogramin resistance gene
+ abrGroups.add("ARO:3000241");// lincosamide resistance gene
+ abrGroups.add("ARO:3000267");// linezolid resistance gene
+ abrGroups.add("ARO:3000271");// fosfomycin resistance gene
+ abrGroups.add("ARO:3000315");// macrolide resistance gene
+ abrGroups.add("ARO:3000362");// mosaic antibiotic resistance gene
+ abrGroups.add("ARO:3000383");// rifampin resistance gene
+ abrGroups.add("ARO:3000398");// chloramphenicol resistance gene
+ abrGroups.add("ARO:3000408");// sulfonamide resistance gene
+ abrGroups.add("ARO:3000468");// ethambutol resistance gene
+ abrGroups.add("ARO:3000472");// tetracycline resistance gene
+ abrGroups.add("ARO:3000477");// aminocoumarin resistance gene
+ abrGroups.add("ARO:3000494");// glycopeptide resistance gene
+ abrGroups.add("ARO:3000529");// mupirocin resistance gene
+ abrGroups.add("ARO:3000751");// peptide antibiotic resistance gene
+ abrGroups.add("ARO:3000868");// streptothricin resistance gene
+ abrGroups.add("ARO:3001217");// trimethoprim resistance gene
+ abrGroups.add("ARO:3001311");// elfamycin resistance gene
+ abrGroups.add("ARO:3002984");// polymyxin resistance gene
+ abrGroups.add("ARO:3003024");// fusidic acid resistance gene
+ abrGroups.add("ARO:3003058");// tunicamycin resistance gene
+ abrGroups.add("ARO:3003073");// lipopeptide antibiotic resistance gene
+ abrGroups.add("ARO:3003252");// bacitracin resistance gene
+ abrGroups.add("ARO:3003432");// isoniazid resistance gene
+ abrGroups.add("ARO:3003433");// pyrazinamide resistance gene
+
+
+ abrGroups.add("ARO:3000004");// class B (metallo-) beta-lactamase
+ abrGroups.add("ARO:3000075");// class D beta-lactamase
+ abrGroups.add("ARO:3000076");// class C beta-lactamase
+ abrGroups.add("ARO:3000078");// class A beta-lactamase
+
+ HashMap> ret= new HashMap>();
+
+ for (Sequence seq : seqs){
+ String ID = seq.getName();
+
+ if (blackList.contains(ID))
+ continue;
+
+ String desc = seq.getDesc();
+ HashSet myGroups = new HashSet();
+
+ int indexPos = 0;
+ while (true){
+ indexPos = desc.indexOf("ARO:",indexPos);
+ if (indexPos < 0)
+ break;
+
+ if (indexPos + 11 < desc.length()){
+ String ARO = desc.substring(indexPos, indexPos + 11);
+ indexPos += 11;
+ if (ARO.equals("ARO:1000001"))
+ continue;
+
+ GoTerm term = odo.terms.get(ARO);
+ resistanceClass(term,abrGroups, myGroups);
+ }
+ }//while
+ System.out.print(seq.getName() + " : ");
+ for (String group:myGroups)
+ System.out.print(group + ";" );
+
+ System.out.println();
+
+ ret.put(seq.getName(), myGroups);
+ }
+ return ret;
+ }
+
+ /**
+ * @param args
+ * @throws IOException
+ *
+ */
+ public static void main(String[] args) throws IOException {
+ odo = new ODOParser("aro.obo");
+ //relationList.add("is_a");
+/********************************************************
+ ArrayList seqs = SequenceReader.readAll(args[0],Alphabet.DNA());
+
+ HashMap allele2Sequence = new HashMap();
+
+ for (Sequence seq:seqs){
+ allele2Sequence.put(seq.getName(), seq);
+ }
+
+ HashMap> allele2Groups = checkGenes(seqs);
+
+
+ BufferedReader bf = SequenceReader.openFile("togene");
+ String line = "";
+
+ HashMap> gene2Groups = new HashMap>();
+ HashMap gene2dbGeneID = new HashMap();
+
+ GeneDatabase geneDB = new GeneDatabase();
+
+
+ while ( (line = bf.readLine())!=null){
+ String [] toks = line.trim().split("\t");
+ String alleleID = toks[0];
+
+ if (!allele2Groups.containsKey(alleleID))
+ continue;
+
+ String geneID = ((toks.length > 2)? toks[2]:"None") + "_" + toks[1];
+
+ HashSet agroup = allele2Groups.get(alleleID);
+ if (gene2Groups.containsKey(geneID)){
+ HashSet ggroup = gene2Groups.get(geneID);
+ if (agroup == null){
+ LOG.error(alleleID + " xx " + geneID );
+ }
+
+ if (!ggroup.containsAll(agroup))
+ LOG.error( alleleID + " <> " + geneID );
+
+ if (!agroup.containsAll(ggroup))
+ LOG.error(alleleID + " # " + geneID );
+
+ GeneDatabase.GeneFamily dbFamily = geneDB.getFamily(gene2dbGeneID.get(geneID));
+
+ if (dbFamily == null)
+ LOG.error("Problem finding " + geneID );
+ else{
+ dbFamily.addSequence(allele2Sequence.get(alleleID));
+ }
+
+
+ }else{
+ gene2Groups.put(geneID, allele2Groups.get(toks[0]));
+ Sequence seq = allele2Sequence.get(alleleID);
+ String famID = geneDB.addNewFamily(seq);
+ gene2dbGeneID.put(geneID, famID);
+ String desc = "geneID=" + geneID + ";dg=";
+
+ for (String group:agroup){
+ desc += group + ",";
+ }
+
+ geneDB.getFamily(famID).setDesc(desc);
+ }
+ }
+
+ geneDB.write2File("F.fasta", false);
+ geneDB.write2File("A.fasta", true);
+
+ LOG.info(gene2Groups.size() + " group ");
+
+
+
+ /********************************************************/
+
+ for (String relID:odo.relTypes.keySet()){
+ System.out.println(relID + "\t" + odo.relTypes.get(relID));
+ }
+
+ for (String termID:odo.terms.keySet()){
+ GoTerm term = odo.terms.get(termID);
+ for (GoRelationship rel:term.relationship){
+ if (rel.relType.ID.equals("is_a"))
+ System.out.println(termID + "(" + term.name +")\t" + rel.relTerm.ID + "(" + rel.relTerm.name+")");
+ }
+
+ conferDrugList.clear();
+ conferList.clear();
+ sb = new StringBuilder();
+
+ odo.confer(term);
+
+ System.out.println(" " + term.ID + " :" + term.name + " " + term.desc);
+
+ //sodo.list(term, relationList);
+ System.out.print(" Confer: ");
+ for (String st:conferList){
+ System.out.print(st + ";");
+ }
+
+ System.out.print("\n Confer drug : ");
+ for (String st:conferDrugList){
+ System.out.print(st + ";");
+ }
+ System.out.print("\n-----------------------------------------------------\n");
+
+ }
+ /********************************************************/
+ }
+}
diff --git a/src/dev/java/japsadev/bio/go/RefFinderParser.java b/src/dev/java/japsadev/bio/go/RefFinderParser.java
new file mode 100644
index 0000000..9582f5b
--- /dev/null
+++ b/src/dev/java/japsadev/bio/go/RefFinderParser.java
@@ -0,0 +1,144 @@
+/*****************************************************************************
+ * Copyright (c) Minh Duc Cao, Monash Uni & UQ, All rights reserved. *
+ * *
+ * Redistribution and use in source and binary forms, with or without *
+ * modification, are permitted provided that the following conditions *
+ * are met: *
+ * *
+ * 1. Redistributions of source code must retain the above copyright notice, *
+ * this list of conditions and the following disclaimer. *
+ * 2. Redistributions in binary form must reproduce the above copyright *
+ * notice, this list of conditions and the following disclaimer in the *
+ * documentation and/or other materials provided with the distribution. *
+ * 3. Neither the names of the institutions nor the names of the contributors*
+ * may be used to endorse or promote products derived from this software *
+ * without specific prior written permission. *
+ * *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS *
+ * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, *
+ * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR *
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR *
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, *
+ * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, *
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR *
+ * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF *
+ * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING *
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
+ ****************************************************************************/
+
+/*****************************************************************************
+ * Revision History
+ * 20 Jul 2015 - Minh Duc Cao: Created
+ *
+ ****************************************************************************/
+package japsadev.bio.go;
+
+import japsa.bio.gene.GeneDatabase;
+import japsa.seq.Alphabet;
+import japsa.seq.Sequence;
+import japsa.seq.SequenceReader;
+import org.slf4j.Logger;
+import org.slf4j.LoggerFactory;
+
+
+import java.io.BufferedReader;
+import java.io.IOException;
+import java.util.ArrayList;
+import java.util.HashMap;
+
+/**
+ * @author minhduc
+ *
+ */
+public class RefFinderParser {
+ private static final Logger LOG = LoggerFactory.getLogger(RefFinderParser.class);
+
+ public RefFinderParser(String fileName) throws IOException{
+
+ }
+
+
+ /**
+ * @param args
+ * @throws IOException
+ *
+ */
+ public static void main(String[] args) throws IOException {
+
+ //Read in all sequences
+ ArrayList seqs = SequenceReader.readAll(args[0],Alphabet.DNA());
+
+ HashMap allele2Sequence = new HashMap();
+
+ for (Sequence seq:seqs){
+ allele2Sequence.put(seq.getName(), seq);
+ }
+
+// / HashMap allele2Groups = new HashMap();
+
+ BufferedReader bf = SequenceReader.openFile("togene");
+ String line = "";
+
+ HashMap gene2Groups = new HashMap();
+ HashMap gene2dbGeneID = new HashMap();
+
+ GeneDatabase geneDB = new GeneDatabase();
+
+
+ while ( (line = bf.readLine())!=null){
+ String [] toks = line.trim().split("\t");
+ String alleleID = toks[0];
+ String agroup = toks[1];
+ String geneID = toks[4] + "_" + toks[2];
+
+ if (gene2Groups.containsKey(geneID)){
+ String ggroup = gene2Groups.get(geneID);
+
+ if (agroup == null){
+ LOG.error(alleleID + " xx " + geneID );
+ }
+
+ if (!ggroup.equals(agroup))
+ LOG.error( alleleID + " <> " + geneID );
+
+ GeneDatabase.GeneFamily dbFamily = geneDB.getFamily(gene2dbGeneID.get(geneID));
+
+ if (dbFamily == null)
+ LOG.error("Problem finding " + geneID );
+ else{
+ dbFamily.addSequence(allele2Sequence.get(alleleID));
+ }
+
+
+ }else{
+ gene2Groups.put(geneID, agroup);
+ Sequence seq = allele2Sequence.get(alleleID);
+ String famID = geneDB.addNewFamily(seq);
+ gene2dbGeneID.put(geneID, famID);
+ String desc = "geneID=" + geneID + ";dg="+agroup;
+
+ geneDB.getFamily(famID).setDesc(desc);
+ }
+ }
+
+ geneDB.write2File("F.fasta", false);
+ geneDB.write2File("A.fasta", true);
+
+ LOG.info(gene2Groups.size() + " group ");
+
+ /********************************************************
+ for (String relID:odo.relTypes.keySet()){
+ System.out.println(relID + "\t" + odo.relTypes.get(relID));
+ }
+
+ for (String termID:odo.terms.keySet()){
+ GoTerm term = odo.terms.get(termID);
+ for (GoRelationship rel:term.relationship){
+ if (rel.relType.ID.equals("is_a"))
+ System.out.println(termID + "(" + term.name +")\t" + rel.relTerm.ID + "(" + rel.relTerm.name+")");
+ }
+ }
+ /********************************************************/
+ }
+}
diff --git a/src/dev/java/japsadev/bio/go/TypeRelationship.java b/src/dev/java/japsadev/bio/go/TypeRelationship.java
new file mode 100644
index 0000000..f655f29
--- /dev/null
+++ b/src/dev/java/japsadev/bio/go/TypeRelationship.java
@@ -0,0 +1,43 @@
+/*****************************************************************************
+ * Copyright (c) Minh Duc Cao, Monash Uni & UQ, All rights reserved. *
+ * *
+ * Redistribution and use in source and binary forms, with or without *
+ * modification, are permitted provided that the following conditions *
+ * are met: *
+ * *
+ * 1. Redistributions of source code must retain the above copyright notice, *
+ * this list of conditions and the following disclaimer. *
+ * 2. Redistributions in binary form must reproduce the above copyright *
+ * notice, this list of conditions and the following disclaimer in the *
+ * documentation and/or other materials provided with the distribution. *
+ * 3. Neither the names of the institutions nor the names of the contributors*
+ * may be used to endorse or promote products derived from this software *
+ * without specific prior written permission. *
+ * *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS *
+ * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, *
+ * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR *
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR *
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, *
+ * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, *
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR *
+ * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF *
+ * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING *
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
+ ****************************************************************************/
+
+/*****************************************************************************
+ * Revision History
+ * 20 Jul 2015 - Minh Duc Cao: Created
+ *
+ ****************************************************************************/
+package japsadev.bio.go;
+
+public class TypeRelationship{
+ String ID, name;
+ public TypeRelationship(String id, String aname){
+ ID = id;
+ name = aname;
+ }
+}
\ No newline at end of file
diff --git a/src/dev/java/japsadev/bio/hts/clustering/GettingTreadsFromFasta.java b/src/dev/java/japsadev/bio/hts/clustering/GettingTreadsFromFasta.java
new file mode 100644
index 0000000..f175847
--- /dev/null
+++ b/src/dev/java/japsadev/bio/hts/clustering/GettingTreadsFromFasta.java
@@ -0,0 +1,173 @@
+package japsadev.bio.hts.clustering;
+
+import java.io.*;
+import java.util.*;
+
+
+
+
+/**
+ * @author buvan.suji
+ *
+ */
+
+public class GettingTreadsFromFasta {
+
+ private String[] description;
+ private String[] sequence;
+
+
+ static ArrayList Rname = new ArrayList();
+ static ArrayList TReads = new ArrayList();
+ static String FnamePath;
+ static int sequenceLength;
+ static String fle;
+
+ public GettingTreadsFromFasta(String filename) {
+ readSequenceFromFile(filename);
+ }
+
+ @SuppressWarnings("resource")
+ void readSequenceFromFile(String file) {
+
+ ArrayList desc = new ArrayList();
+ ArrayList seq = new ArrayList();
+ BufferedReader in = null;
+
+ try {
+ in = new BufferedReader(new FileReader(file));
+ StringBuffer buffer = new StringBuffer();
+ String line = in.readLine();
+
+ if (line == null)
+ throw new IOException(file + " is an empty file");
+
+ if (line.charAt(0) != '>')
+ throw new IOException("First line of " + file
+ + " should start with '>'");
+ else
+ desc.add(line);
+
+ for (line = in.readLine().trim(); line != null; line = in
+ .readLine()) {
+ if (line.length() > 0 && line.charAt(0) == '>') {
+ seq.add(buffer.toString());
+ buffer = new StringBuffer();
+ desc.add(line);
+ } else
+ buffer.append(line.trim());
+ }
+ if (buffer.length() != 0)
+ seq.add(buffer.toString());
+ } catch (IOException e) {
+ System.out.println("Error when reading " + file);
+ e.printStackTrace();
+ }
+
+ description = new String[desc.size()];
+ sequence = new String[seq.size()];
+ for (int i = 0; i < seq.size(); i++) {
+ description[i] = (String) desc.get(i);
+ sequence[i] = (String) seq.get(i);
+ }
+ sequenceLength = seq.size();
+
+ }
+
+ // return first sequence as a String
+ public String getSequence() {
+ return sequence[0];
+ }
+
+ // return first xdescription as String
+ public String getDescription() {
+ return description[0];
+ }
+
+ // return sequence as a String
+ public String getSequence(int i) {
+ return sequence[i];
+ }
+
+ // return description as String
+ public String getDescription(int i) {
+ return description[i];
+ }
+
+ public int size() {
+ return sequence.length;
+ }
+
+ public static void FileReading(String filename) {
+ GettingTreadsFromFasta fsf = new GettingTreadsFromFasta(
+ filename);
+ StringBuffer buffer1 = new StringBuffer();
+
+ for (int i = 0; i < fsf.size(); i++) {
+ String temp = ((buffer1.append(fsf.getDescription(i)))
+ .deleteCharAt(0)).toString();
+ Rname.add(temp);
+ TReads.add(fsf.getSequence(i));
+ buffer1 = new StringBuffer();
+ }
+ }
+
+ // call this method-1
+ //public static void DestReads() throws Exception {
+ public static void DestReads(String fle1) throws Exception {
+ /*String fname = "";
+ System.out.print("Enter the name of the FastaFile:");
+ fname = (new BufferedReader(new InputStreamReader(System.in)))
+ .readLine();
+ FnamePath = new File(fname).getName();*/
+
+ // System.out.println(FnamePath);
+ //FileReading(fname);
+ fle = fle1;
+ FnamePath = fle1;
+ //FnamePath = "chr10134169759_134174717.fasta";
+ FileReading(FnamePath);
+ }
+
+ public static void main(String[] args) throws Exception {
+ DestReads(fle);
+ ViewList();
+ }
+
+ // call this method-2
+ public static ArrayList GetRname() {
+ return Rname;
+ }
+
+ // call this method-3
+ public static ArrayList GetTReads() {
+ return TReads;
+ }
+
+ public static int NumberReads() {
+ return TReads.size();
+ }
+
+ // call this method-4
+ public static String GetFileName() {
+ return FnamePath;
+ }
+
+ public static int SeqLength() {
+ return sequenceLength;
+ }
+
+ public static int NelementsClustering() {
+ return sequenceLength * (sequenceLength - 1) / 2;
+ }
+
+ public static void ViewList() {
+ System.out.println(GetRname());
+ // System.out.println(GetTReads());
+ System.out.println(SeqLength());
+ System.out.println(NelementsClustering());
+ }
+
+
+
+}
diff --git a/src/dev/java/japsadev/bio/hts/clustering/KmeanClustering.java b/src/dev/java/japsadev/bio/hts/clustering/KmeanClustering.java
new file mode 100644
index 0000000..ddd55f9
--- /dev/null
+++ b/src/dev/java/japsadev/bio/hts/clustering/KmeanClustering.java
@@ -0,0 +1,300 @@
+package japsadev.bio.hts.clustering;
+
+import java.util.ArrayList;
+import java.io.*;
+
+import japsa.seq.SequenceOutputStream;
+import japsadev.bio.hts.clustering.PairDistance;
+import japsadev.bio.hts.clustering.GettingTreadsFromFasta;
+
+/**
+ * @author buvan.suji
+ *
+ */
+
+public class KmeanClustering {
+
+ public static void Clustering() throws Exception{
+ FileInputStream file1 = new FileInputStream("TRfile.fasta");
+ BufferedReader br = new BufferedReader(new InputStreamReader(file1));
+ String line = null;
+ while((line = br.readLine())!=null){
+ ArrayList descript = new ArrayList();
+ ArrayList reads = new ArrayList();
+ //ArrayList MaxReadLn = new ArrayList();
+ int MaxReadLn = 0;
+ int Nreads;
+ String FileName;
+ int NumberElements;
+ int ClustElements;
+ int n;
+ int d[];
+ int k[][];
+ final int p =2;//number of clusters
+ int tempk[][];
+ int m[];
+ int Nclusters = 0;
+
+ double max = 0;
+ double temp1 = 0;
+ int index1=0;
+ int index2=0;
+ int count1=0,count2=0;
+ long startTime = System.nanoTime();
+
+
+ GettingTreadsFromFasta.DestReads(line);
+ //GettingTreadsFromFasta.DestReads();
+ descript = GettingTreadsFromFasta.GetRname();
+ reads = GettingTreadsFromFasta.GetTReads();
+ Nreads = GettingTreadsFromFasta.NumberReads();
+ FileName = GettingTreadsFromFasta.GetFileName();
+ NumberElements = GettingTreadsFromFasta.SeqLength();
+ ClustElements = GettingTreadsFromFasta.NelementsClustering();
+
+ double[][] table = new double[Nreads][Nreads];
+
+ for (int i=0; im[i]){
+ diff[i] = list[m[i]][a];
+ }
+ else{
+ diff[i]=0;
+ }
+ }
+
+ int val=0;
+ double temp=diff[0];
+ for(int i=0;i<2;++i){
+ if(diff[i] t = new ArrayList() ;
+
+ for(int i=0;i<2;++i){
+ m[i]=0; // initializing means to 0
+ }
+
+ //int cnt=0;
+ for(int i=0;i<2;++i){
+
+ for(int j=0;j s = new ArrayList();
+
+ for(int x1 = 0; x1t.get(x2)){
+ sum = sum+(list[t.get(x2)][t.get(x1)]*list[t.get(x2)][t.get(x1)]);
+ }
+ else{
+ sum = sum+0;
+ }
+ }
+ //System.out.println(sum);
+ s.add(sum/t.size());
+ }
+ //System.out.println(s.size());
+ double min=s.get(0);
+ int d2 = 0;
+ for(int x3=1;x3s.get(x3)){
+ min = s.get(x3);
+ d2 = x3;
+ }
+ }
+
+
+ m[i]=t.get(d2);
+ s.clear();
+ t.clear();
+
+ }
+ }
+
+//This checks if previous k ie. tempk and current k are same.Used as terminating case.
+ static int VerifyEqual(int n, int k[][], int tempk[][]){
+ for(int i=0;i<2;++i){
+ for(int j=0;j>
+ Clustering(ArrayList reads) throws Exception{
+ ArrayList> tempClusterList =
+ new ArrayList>();
+ ArrayList tempList = new ArrayList();
+ if(reads.size()<=2){
+ tempClusterList.add(tempList);
+ tempClusterList.add(tempList);
+ tempClusterList.add(tempList);
+ return tempClusterList;
+ }
+ else{
+ int MaxReadLn = 0;
+ int MinReadLn = 1000000;
+
+ int n;
+ int d[];
+ int k[][];
+ final int p =2;//number of clusters
+ int tempk[][];
+ int m[];
+ int Nclusters = 0;
+
+ double max = 0;
+ double temp1 = 0;
+ int index1=0;
+ int index2=0;
+ int count1=0,count2=0;
+ long startTime = System.nanoTime();
+
+ n = reads.size();
+
+ double[][] table = new double[n][n];
+ ArrayList readsLengthRange = new ArrayList();
+
+ for (int x=0; x> clusterList =
+ new ArrayList>();
+ clusterList.add(readsLengthRange);
+
+ for(int x=0;x tempcluster = new ArrayList();
+ for(int y=0;k[x][y]!=-1 && ym[i]){
+ diff[i] = list[m[i]][a];
+ }
+ else{
+ diff[i]=0;
+ }
+ }
+
+ int val=0;
+ double temp=diff[0];
+ for(int i=0;i<2;++i){
+ if(diff[i] t = new ArrayList() ;
+
+ for(int i=0;i<2;++i){
+ m[i]=0; // initializing means to 0
+ }
+
+ //int cnt=0;
+ for(int i=0;i<2;++i){
+
+ for(int j=0;j1){
+ ArrayList s = new ArrayList();
+
+ for(int x1 = 0; x1t.get(x2)){
+ sum = sum+(list[t.get(x2)][t.get(x1)]*list[t.get(x2)][t.get(x1)]);
+ }
+ else{
+ sum = sum+0;
+ }
+ }
+ //System.out.println(sum);
+ s.add(sum/t.size());
+ }
+ //System.out.println(s.size());
+ double min=s.get(0);
+ int d2 = 0;
+ for(int x3=1;x3s.get(x3)){
+ min = s.get(x3);
+ d2 = x3;
+ }
+ }
+
+
+ m[i]=t.get(d2);
+ s.clear();
+ t.clear();
+ }else{
+ m[i]=t.get(0);
+ }
+
+ }
+ }
+
+//This checks if previous k ie. tempk and current k are same.Used as terminating case.
+ static int VerifyEqual(int n, int k[][], int tempk[][]){
+ for(int i=0;i<2;++i){
+ for(int j=0;j parentMap = new HashMap<>();
+
+ for (int i = 1; i <= m; ++i) {
+ d[i][0] = i;
+ }
+
+ for (int j = 1; j <= n; ++j) {
+ d[0][j] = j;
+ }
+
+
+
+
+ for (int j = 1; j <= n; ++j) {
+ for (int i = 1; i <= m; ++i) {
+ final int delta = (s.charAt(j - 1) == z.charAt(i - 1)) ? 0 : 1;
+
+ double tentativeDistance = d[i - 1][j] +0.5;//gap penalty
+ EditOperation editOperation = EditOperation.INSERT;
+
+ if (tentativeDistance > d[i][j - 1] +0.5) {
+ tentativeDistance = d[i][j - 1] +0.5;
+ editOperation = EditOperation.DELETE;
+ }
+
+ if (tentativeDistance > d[i - 1][j - 1] + delta) {
+ tentativeDistance = d[i - 1][j - 1] + delta;
+ editOperation = EditOperation.SUBSTITUTE;
+ }
+
+ d[i][j] = tentativeDistance;
+
+
+ switch (editOperation) {
+ case SUBSTITUTE:
+ parentMap.put(new Point(i, j), new Point(i - 1, j - 1));
+ break;
+
+ case INSERT:
+ parentMap.put(new Point(i, j), new Point(i - 1, j));
+ break;
+
+ case DELETE:
+ parentMap.put(new Point(i, j), new Point(i, j - 1));
+ break;
+ default:
+ break;
+ }
+ }
+ }
+
+ final StringBuilder topLineBuilder = new StringBuilder(n + m);
+ final StringBuilder bottomLineBuilder = new StringBuilder(n + m);
+ final StringBuilder editSequenceBuilder = new StringBuilder(n + m);
+ Point current = new Point(m, n);
+
+ while (true) {
+ Point predecessor = parentMap.get(current);
+
+ if (predecessor == null) {
+ break;
+ }
+
+ if (current.x != predecessor.x && current.y != predecessor.y) {
+ final char schar = s.charAt(predecessor.y);
+ final char zchar = z.charAt(predecessor.x);
+
+ topLineBuilder.append(schar);
+ bottomLineBuilder.append(zchar);
+ editSequenceBuilder.append(schar != zchar ?
+ EditOperation.SUBSTITUTE :
+ EditOperation.MATCH);
+ } else if (current.x != predecessor.x) {
+ topLineBuilder.append(GAP);
+ bottomLineBuilder.append(z.charAt(predecessor.x));
+ editSequenceBuilder.append(EditOperation.INSERT);
+ } else {
+ topLineBuilder.append(s.charAt(predecessor.y));
+ bottomLineBuilder.append(GAP);
+ editSequenceBuilder.append(EditOperation.DELETE);
+ }
+
+ current = predecessor;
+ }
+
+ // Remove the last characters that correspond to the very beginning
+ // of the alignments and edit sequence (since the path reconstruction
+ // proceeds from the "end" to the "beginning" of the distance matrix.
+ topLineBuilder .deleteCharAt(topLineBuilder.length() - 1);
+ bottomLineBuilder .deleteCharAt(bottomLineBuilder.length() - 1);
+ editSequenceBuilder.deleteCharAt(editSequenceBuilder.length() - 1);
+
+ // Our result data is backwards, reverse them.
+ topLineBuilder .reverse();
+ bottomLineBuilder .reverse();
+ editSequenceBuilder.reverse();
+
+ return new EditDistanceResult(d[m][n],
+ editSequenceBuilder.toString(),
+ topLineBuilder.toString(),
+ bottomLineBuilder.toString());
+ }
+
+
+
+ public static void main(String[] args) throws IOException {
+ EditDistanceResult result = compute("ACTAGGTTA", "TAAGGCTCAAT");
+ System.out.println("Distance: " + result.getDistance());
+ System.out.println("Edit sequence: " + result.getEditSequence());
+ System.out.println("Alignment:");
+ System.out.println(result.getTopAlignmentRow());
+ System.out.println(result.getBottomAlignmentRow());
+
+
+ }
+
+
+
+
+}
diff --git a/src/dev/java/japsadev/bio/hts/clustering/WriteClusterResultOnFile.java b/src/dev/java/japsadev/bio/hts/clustering/WriteClusterResultOnFile.java
new file mode 100644
index 0000000..fe5c98a
--- /dev/null
+++ b/src/dev/java/japsadev/bio/hts/clustering/WriteClusterResultOnFile.java
@@ -0,0 +1,73 @@
+package japsadev.bio.hts.clustering;
+
+import java.util.ArrayList;
+import java.io.*;
+
+import japsa.seq.SequenceOutputStream;
+import japsadev.bio.hts.clustering.PairDistance;
+import japsadev.bio.hts.clustering.GettingTreadsFromFasta;
+import japsa.seq.*;
+
+/**
+ * @author buvan.suji
+ *
+ */
+
+public class WriteClusterResultOnFile {
+
+ public static void writeOnFile(ArrayList>
+ list, Sequence cons1, Sequence cons2, String FileName ) throws Exception{
+
+ File file = new File("ClusterResult_"+FileName+".fasta");
+ FileWriter fw = new FileWriter(file.getAbsoluteFile());
+ BufferedWriter bw = new BufferedWriter(fw);
+
+
+ ArrayList list1 = new ArrayList();
+ list1 = list.get(0);
+
+ bw.write("Minimum Read Length: "+ list1.get(0));
+ bw.newLine();
+ bw.write("Maximum Read Length: "+ list1.get(1));
+ bw.newLine();
+ bw.write("Estimated Time: "+ list1.get(2));
+ bw.newLine();
+ bw.newLine();
+ bw.write("The consensus sequence of cluster1: ");
+ bw.newLine();
+ bw.write(""+cons1);
+ bw.newLine();
+ bw.newLine();
+ bw.write("The C1 members are: ");
+ bw.newLine();
+ for(int x=1;x tempList = new ArrayList();
+ tempList = list.get(x);
+
+ for(int y=0;y vertices = (HashSet) graph.getVertices();
+//
+// //display the initial setup- all vertices adjacent to each other
+// for(Vertex vertex:vertices){
+// System.out.println(vertex);
+//
+// for(int j = 0; j < vertex.getNeighborCount(); j++){
+// System.out.println(vertex.getNeighbor(j));
+// }
+//
+// System.out.println();
+// }
+//
+// Node source=new Node(graph.getVertex("108"),false),
+// dest=new Node(graph.getVertex("201"),true);
+// int distance=2962;
+//
+// ArrayList paths=graph.DFS(source, dest, distance);
+// if(!paths.isEmpty()){
+// System.out.println("Paths found ("+distance+"):");
+// for(Path p:paths)
+// System.out.println(p.toString() + " d=" + p.getDeviation() );
+// }
+// else
+// System.out.println("Path not found ("+distance+")!");
+
+ }catch (Exception e) {
+ // TODO Auto-generated catch block
+ e.printStackTrace();
+ }
+
+// Vertex v1= new Vertex("EDGE_1_length_1000_cov_9"),
+// v2= new Vertex("EDGE_2_length_200_cov_22"),
+// v3= new Vertex("EDGE_3_length_300_cov_12"),
+// v4= new Vertex("EDGE_4_length_500_cov_33"),
+// v5= new Vertex("EDGE_5_length_600_cov_21");
+// graph.addVertex(v1, false);
+// graph.addVertex(v2, false);
+// graph.addVertex(v3, false);
+// graph.addVertex(v4, false);
+// graph.addVertex(v5, false);
+//
+// graph.addEdge(v1, v2, true, true);
+// graph.addEdge(v1, v2, false, false);
+//
+// graph.addEdge(v2, v1, true, true);
+// graph.addEdge(v2, v3, true, true);
+// graph.addEdge(v2, v1, false, false);
+// graph.addEdge(v2, v4, false, false);
+//
+// graph.addEdge(v3, v4, true, true);
+// graph.addEdge(v3, v2, false, false);
+//
+// graph.addEdge(v4, v5, true, true);
+// graph.addEdge(v4, v2, true, true);
+// graph.addEdge(v4, v3, false, false);
+// graph.addEdge(v4, v5, false, false);
+//
+// graph.addEdge(v5, v4, false, false);
+// graph.addEdge(v5, v4, true, true);
+//
+// graph.printStats();
+//
+// Path p1=new Path(graph, "1+,2+,3+"),
+// p2=new Path(graph, "3+,4+");
+// graph.reduce(p1);
+// graph.printStats();
+// graph.reduce(p2);
+// graph.printStats();
+ }
+}
+
+
diff --git a/src/dev/java/japsadev/bio/hts/metagenome/Edge.java b/src/dev/java/japsadev/bio/hts/metagenome/Edge.java
new file mode 100644
index 0000000..99cb9a3
--- /dev/null
+++ b/src/dev/java/japsadev/bio/hts/metagenome/Edge.java
@@ -0,0 +1,180 @@
+package japsadev.bio.hts.metagenome;
+
+/**
+ * This class models an bidirected Edge in my Graph implementation.
+ * An Edge contains two vertices and a weight (distance between them).
+ * A certain edge (v1,v2) can take one among 4 types: ++, --, +- and -+. Each
+ * type corresponds to the way we read the DNA sequence in each read when traversing
+ * this edge.
+ * For example: v1->---<-v2 or (v1,v2)+- spells out (v1 v2') and/or (v2 v1') as in SPAdes output.
+ * This class also deviates from the expectations of the Comparable interface
+ * in that a return value of 0 does not indicate that this.equals(other). The
+ * equals() method only compares the vertices, while the compareTo() method
+ * compares the edge weights. This provides more efficient implementation for
+ * checking uniqueness of edges, as well as the fact that two edges of equal weight
+ * should be considered equitably in a path finding or spanning tree algorithm.
+ *
+ * @author Son Nguyen
+ * @date August 20, 2016
+ */
+public class Edge implements Comparable {
+
+ private Vertex one, two;
+ private boolean dOne, dTwo;
+ private int weight;
+
+ /**
+ *
+ * @param one The first vertex in the Edge
+ * @param two The second vertex in the Edge
+ */
+ public Edge(Vertex one, Vertex two, boolean d1, boolean d2){
+ this(one, two, d1, d2, -127);
+ }
+
+ /**
+ *
+ * @param one The first vertex in the Edge
+ * @param two The second vertex of the Edge
+ * @param weight The weight of this Edge
+ */
+ public Edge(Vertex one, Vertex two, boolean dOne, boolean dTwo, int weight){
+ //this.one = (one.getLabel().compareTo(two.getLabel()) <= 0) ? one : two;
+ //this.two = (this.one == one) ? two : one;
+ this.one=one;
+ this.two=two;
+ this.weight = weight;
+ this.dOne=dOne;
+ this.dTwo=dTwo;
+ }
+
+
+ /**
+ *
+ * @param current
+ * @return The neighbor of current along this Edge
+ */
+ public Vertex getNeighbor(Vertex current){
+ if(!(current.equals(one) || current.equals(two))){
+ return null;
+ }
+
+ return (current.equals(one)) ? two : one;
+ }
+ /**
+ * Return the same Edge but reading the other way around
+ * just swap the order of its vertices upside down
+ * @param
+ * @return the identical Edge
+ */
+ public Edge getReversedRead(){
+ return new Edge(this.two, this.one, !this.dTwo, !this.dOne, this.weight);
+ }
+ /**
+ *
+ * @param current
+ * @return The direction to spell *current* along this Edge
+ */
+ public boolean getDirection(Vertex current){
+ assert (current.equals(one) || current.equals(two)):"Vertex doesn't belong to this Edge!";
+
+ return (current.equals(one)) ? dOne : !dTwo;
+ }
+
+ /**
+ *
+ * @return Vertex this.one
+ */
+ public Vertex getOne(){
+ return this.one;
+ }
+
+ /**
+ *
+ * @return Vertex this.two
+ */
+ public Vertex getTwo(){
+ return this.two;
+ }
+
+ /**
+ *
+ * @return boolean this.dOne
+ */
+ public boolean getDOne(){
+ return this.dOne;
+ }
+
+ /**
+ *
+ * @return boolean this.dTwo
+ */
+ public boolean getDTwo(){
+ return this.dTwo;
+ }
+ /**
+ *
+ * @return int The weight of this Edge
+ */
+ public int getWeight(){
+ return this.weight;
+ }
+
+
+ /**
+ *
+ * @param weight The new weight of this Edge
+ */
+ public void setWeight(int weight){
+ this.weight = weight;
+ }
+
+
+ /**
+ * Note that the compareTo() method deviates from
+ * the specifications in the Comparable interface. A
+ * return value of 0 does not indicate that this.equals(other).
+ * The equals() method checks the Vertex endpoints, while the
+ * compareTo() is used to compare Edge weights
+ *
+ * @param other The Edge to compare against this
+ * @return int this.weight - other.weight
+ */
+ public int compareTo(Edge other){
+ return this.weight - other.weight;
+ }
+
+ /**
+ *
+ * @return String A String representation of this Edge
+ */
+ public String toString(){
+ return "({" + one + (dOne?">":"<") + ", " + two + (dTwo?">":"<") + "}, " + weight + ")";
+ }
+
+ /**
+ *
+ * @return int The hash code for this Edge
+ */
+ public int hashCode(){
+ return (one.getLabel() + (dOne?"":"'") + two.getLabel() + (dTwo?"":"'")).hashCode();
+ }
+
+ /**
+ *
+ * @param other The Object to compare against this
+ * @return true iff other is an Edge with the same Vertices as this
+ */
+ public boolean equals(Object other){
+ if(!(other instanceof Edge)){
+ return false;
+ }
+
+ Edge e = (Edge)other;
+
+ return (e.one.equals(this.one) && e.two.equals(this.two) && (e.getDOne()==this.dOne) && (e.getDTwo()==this.dTwo))
+ || (e.one.equals(this.two) && e.two.equals(this.one) && (e.getDOne()!=this.dOne) && (e.getDTwo()!=this.dTwo));
+ }
+}
+
+
diff --git a/src/dev/java/japsadev/bio/hts/metagenome/Graph.java b/src/dev/java/japsadev/bio/hts/metagenome/Graph.java
new file mode 100644
index 0000000..d419852
--- /dev/null
+++ b/src/dev/java/japsadev/bio/hts/metagenome/Graph.java
@@ -0,0 +1,389 @@
+package japsadev.bio.hts.metagenome;
+
+import java.io.IOException;
+import java.util.*;
+import japsa.seq.Alphabet;
+import japsa.seq.FastaReader;
+import japsa.seq.Sequence;
+import japsa.seq.SequenceReader;
+
+/**
+ * This class models a simple, bidirected graph using an
+ * incidence list representation. Vertices are identified
+ * uniquely by their labels, and only unique vertices are allowed.
+ * At most one unique Edge per vertex pair is allowed in this Graph.
+ *
+ * @author Son Nguyen
+ * @date August 20, 2016
+ */
+public class Graph {
+
+ private HashMap vertices;
+ private HashMap edges;
+ private int kmer;
+
+ static final int TOLERATE=500;
+
+
+ public Graph(){
+ this.vertices = new HashMap();
+ this.edges = new HashMap();
+ setKmerSize(127);//default kmer size used by SPAdes to assembly MiSeq data
+ }
+
+
+ public Graph(String graphFile) throws IOException{
+ this();
+ //1. next iterate over again to read the connections
+ SequenceReader reader = new FastaReader(graphFile);
+ Sequence seq;
+ int shortestLen = 10000;
+ while ((seq = reader.nextSequence(Alphabet.DNA())) != null){
+ if(seq.length() 1){
+ String[] nbList = adjList[1].split(",");
+ for(int i=0; i < nbList.length; i++){
+ // create list of bridges here (distance=-kmer overlapped)
+ String neighbor = nbList[i];
+ boolean dir2=neighbor.contains("'")?false:true;
+ neighbor=neighbor.replaceAll("[^a-zA-Z0-9_.]", "").trim();
+
+ Vertex nbVertex=new Vertex(neighbor);
+ if(getVertex(nbVertex.getLabel())!=null)
+ nbVertex=getVertex(nbVertex.getLabel());
+
+ addVertex(nbVertex, false);
+
+ addEdge(current, nbVertex, dir1, dir2);
+ }
+ }
+
+ }
+ //rough estimation of kmer used
+ if((shortestLen-1) != getKmerSize())
+ setKmerSize(shortestLen-1);
+
+ reader.close();
+ }
+ /**
+ * This constructor accepts an ArrayList and populates
+ * this.vertices. If multiple Vertex objects have the same label,
+ * then the last Vertex with the given label is used.
+ *
+ * @param vertices The initial Vertices to populate this Graph
+ */
+ public Graph(ArrayList vertices){
+ this.vertices = new HashMap();
+ this.edges = new HashMap();
+
+ for(Vertex v: vertices){
+ this.vertices.put(v.getLabel(), v);
+ }
+ setKmerSize(127);//default kmer size used by SPAdes to assembly MiSeq data
+
+ }
+
+ public int getKmerSize(){
+ return this.kmer;
+ }
+ public void setKmerSize(int kmer){
+ this.kmer=kmer;
+ Path.setK(kmer);
+ }
+ /**
+ * This method adds am edge between Vertices one and two
+ * and their corresponding direction of weight kmer,
+ * if no Edge between these Vertices already exists in the Graph.
+ *
+ * @param one The first vertex to add
+ * @param two The second vertex to add
+ * @param d1 The direction on the side of vertex one
+ * @param d2 The direction on the side of vertex two
+ * @return true iff no Edge relating one and two exists in the Graph
+ */
+ public boolean addEdge(Vertex one, Vertex two, boolean d1, boolean d2){
+ return addEdge(one, two, d1, d2, -kmer);
+ }
+
+
+ /**
+ * Accepts two vertices, their directions and a weight, and adds the edge
+ * ({one, two}, {d1, d2}, weight) iff no Edge relating one and two
+ * exists in the Graph.
+ *
+ * @param one The first Vertex of the Edge
+ * @param two The second Vertex of the Edge
+ * @param d1 The direction on the side of vertex one
+ * @param d2 The direction on the side of vertex two
+ * @param weight The weight of the Edge
+ * @return true iff no Edge already exists in the Graph
+ */
+ public boolean addEdge(Vertex one, Vertex two, boolean d1, boolean d2, int weight){
+
+ //ensures the Edge is not in the Graph
+ Edge e = new Edge(one, two, d1, d2, weight);
+ if(edges.containsKey(e.hashCode()) || edges.containsKey(e.getReversedRead().hashCode())){
+ return false;
+ }
+
+ //and that the Edge isn't already incident to one of the vertices
+ else if(one.containsNeighbor(e) || two.containsNeighbor(e.getReversedRead())){
+ return false;
+ }
+
+ edges.put(e.hashCode(), e);
+ one.addNeighbor(e);
+ two.addNeighbor(e.getReversedRead());
+ return true;
+ }
+
+ /**
+ *
+ * @param e The Edge to look up
+ * @return true iff this Graph contains the Edge e
+ */
+ public boolean containsEdge(Edge e){
+ if(e.getOne() == null || e.getTwo() == null){
+ return false;
+ }
+
+ return this.edges.containsKey(e.hashCode())
+ || this.edges.containsKey(e.getReversedRead().hashCode());
+ }
+
+
+ /**
+ * This method removes the specified Edge from the Graph,
+ * including as each vertex's incidence neighborhood.
+ *
+ * @param e The Edge to remove from the Graph
+ * @return Edge The Edge removed from the Graph
+ */
+ public Edge removeEdge(Edge e){
+ e.getOne().removeNeighbor(e);
+ e.getTwo().removeNeighbor(e.getReversedRead());
+ Edge rmEdge = this.edges.remove(e.hashCode());
+ if (rmEdge==null)
+ rmEdge = this.edges.remove(e.getReversedRead().hashCode());
+ return rmEdge;
+ }
+
+ /**
+ *
+ * @param vertex The Vertex to look up
+ * @return true iff this Graph contains vertex
+ */
+ public boolean containsVertex(Vertex vertex){
+ return this.vertices.get(vertex.getLabel()) != null;
+ }
+
+ /**
+ *
+ * @param label The specified Vertex label
+ * @return Vertex The Vertex with the specified label
+ */
+ public Vertex getVertex(String label){
+ return vertices.get(label);
+ }
+
+ /**
+ * This method adds a Vertex to the graph. If a Vertex with the same label
+ * as the parameter exists in the Graph, the existing Vertex is overwritten
+ * only if overwriteExisting is true. If the existing Vertex is overwritten,
+ * the Edges incident to it are all removed from the Graph.
+ *
+ * @param vertex
+ * @param overwriteExisting
+ * @return true iff vertex was added to the Graph
+ */
+ public boolean addVertex(Vertex vertex, boolean overwriteExisting){
+ Vertex current = this.vertices.get(vertex.getLabel());
+ if(current != null){
+ if(!overwriteExisting){
+ return false;
+ }
+
+ while(current.getNeighborCount() > 0){
+ this.removeEdge(current.getNeighbor(0));
+ }
+ }
+
+
+ vertices.put(vertex.getLabel(), vertex);
+ return true;
+ }
+
+ /**
+ *
+ * @param label The label of the Vertex to remove
+ * @return Vertex The removed Vertex object
+ */
+ public Vertex removeVertex(String label){
+ Vertex v = vertices.remove(label);
+
+ while(v.getNeighborCount() > 0){
+ this.removeEdge(v.getNeighbor((0)));
+ }
+
+ return v;
+ }
+
+ /**
+ *
+ * @return Set All Graph's Vertex objects
+ */
+ public Set getVertices(){
+ return new HashSet(this.vertices.values());
+ }
+
+ /**
+ *
+ * @return Set The Edges of this graph
+ */
+ public Set getEdges(){
+ return new HashSet(this.edges.values());
+ }
+
+ /**
+ * Find a path between two nodes within a given distance
+ */
+ public ArrayList DFS(Node source, Node dest, int distance){
+ System.out.println("Looking for path between " + source.toString() + " to " + dest.toString() + " with distance " + distance);
+ Path tmp = new Path();
+ ArrayList retval = new ArrayList();
+ tmp.addNode(source);
+
+ //traverse(tmp, dest, retval, distance+source.getSeq().length()+dest.getSeq().length());
+ traverse(tmp, dest, retval, distance);
+
+ return retval;
+ }
+
+ public void traverse(Path path, Node dest, ArrayList curResult, int distance){
+ Node source=path.getEnd();
+ assert source!=null:"Path null fault!";
+
+ ArrayList nList = source.getVertex().getNeighbors();
+ for(Edge e:nList){
+ if(e.getDOne()==source.getDirection()){
+ path.addNode(e.getTwo(), e.getDTwo());
+
+ if(e.getTwo()==dest.getVertex() && e.getDTwo()==dest.getDirection() && Math.abs(distance+getKmerSize()) < TOLERATE){
+
+ Path curPath=curResult.isEmpty()?new Path():curResult.get(0), //the best path saved among all possible paths from the list curResult
+ tmpPath=new Path();
+ tmpPath.setComp(path.getNodes());
+ tmpPath.setDeviation(Math.abs(distance+getKmerSize()));
+ if( Math.abs(distance+getKmerSize()) < curPath.getDeviation() )
+ curResult.add(0, tmpPath);
+ else
+ curResult.add(tmpPath);
+
+ System.out.println("Hit added: "+path+"(candidate deviation: "+Math.abs(distance+getKmerSize())+")");
+ }else{
+ int newDistance=distance-e.getTwo().getSequence().length()+getKmerSize();
+ if (newDistance+getKmerSize()<-TOLERATE){
+ System.out.println("Stop following path with distance "+newDistance+" already! : "+path);
+ }else
+ traverse(path, dest, curResult, newDistance);
+ }
+ path.removeLast();
+ }
+ }
+ }
+ /**
+ *
+ * @param p Path to be grouped as a virtually vertex
+ */
+ public void reduce(Path p){
+ Vertex comp=new Vertex(p);
+ //add the new composite Vertex to the graph
+ addVertex(comp, true);
+ //remove unique nodes on p
+ ArrayList tobeRemoved=new ArrayList();
+ for(Node n:p.getNodes()){
+ if(n.getVertex().isUnique())
+ tobeRemoved.add(n.getVertex().getLabel());
+ }
+
+ Node start = p.getStart(),
+ end = p.getEnd();
+ //set neighbors of the grouped Vertex
+ for(Edge e:start.getVertex().getNeighbors()){
+ if(e.getDOne()!=start.getDirection()){
+ //comp.addNeighbor(e);
+ addEdge(comp,e.getTwo(),e.getDOne(),e.getDTwo());
+ }
+ }
+ for(Edge e:end.getVertex().getNeighbors()){
+ if(e.getDOne()!=end.getDirection())
+ //comp.addNeighbor(e);
+ addEdge(comp,e.getTwo(),e.getDOne(),e.getDTwo());
+ }
+ for(String lab:tobeRemoved)
+ removeVertex(lab);
+ //TODO: remove bubbles...
+ }
+ /**
+ *
+ * @param v Vertex to be reverted (1-level reverting)
+ */
+ public void revert(Vertex v){
+ //TODO: revert to initial status by extracting a complex vertex into its initial components
+ Path p=v.getSubComps();
+ if(!containsVertex(v)||p==null) return;
+ //add back all vertices first
+ for(Node n:p.getNodes())
+ addVertex(n.getVertex(), false);
+ //then add back all neighbor edges of this composite vertex
+ for(Edge e:v.getNeighbors())
+ addEdge(v,e.getTwo(),e.getDOne(),e.getDTwo());
+ //finally add back all edges from the path
+ Node prev=p.getStart();
+ for(Node cur:p.getNodes()){
+ if(cur==p.getStart())
+ continue;
+ else{
+ addEdge(prev.getVertex(),cur.getVertex(),prev.getDirection(),cur.getDirection());
+ prev=cur;
+ }
+ }
+
+ //remove the original composite vertex
+ removeVertex(v.getLabel());
+ }
+ public void printStats(){
+ System.out.println(vertices.size() + " vertices:");
+ for(String label:vertices.keySet())
+ System.out.print(label+", ");
+ System.out.println();
+ System.out.println(edges.size() + " edges:");
+ for(Edge e:edges.values()){
+ System.out.print(e.toString()+", ");
+ }
+ System.out.println();
+
+ }
+}
+
+
diff --git a/src/dev/java/japsadev/bio/hts/metagenome/Node.java b/src/dev/java/japsadev/bio/hts/metagenome/Node.java
new file mode 100644
index 0000000..d6c1036
--- /dev/null
+++ b/src/dev/java/japsadev/bio/hts/metagenome/Node.java
@@ -0,0 +1,34 @@
+package japsadev.bio.hts.metagenome;
+
+import japsa.seq.Alphabet;
+import japsa.seq.Sequence;
+
+public class Node{
+ Vertex v;
+ boolean dir;
+ Node(Vertex v, boolean dir){
+ this.v=v;
+ this.dir=dir;
+ }
+ public Vertex getVertex(){
+ return v;
+ }
+ public void setVertex(Vertex v){
+ this.v = v;
+ }
+ public boolean getDirection(){
+ return dir;
+ }
+ public void setDirection(boolean dir){
+ this.dir=dir;
+ }
+ public Node getRC(){
+ return new Node(v,!dir);
+ }
+ public Sequence getSeq(){
+ return dir?v.getSequence():Alphabet.DNA.complement(v.getSequence());
+ }
+ public String toString(){
+ return v.getLabel()+ (dir?"+":"-");
+ }
+}
diff --git a/src/dev/java/japsadev/bio/hts/metagenome/Path.java b/src/dev/java/japsadev/bio/hts/metagenome/Path.java
new file mode 100644
index 0000000..5474435
--- /dev/null
+++ b/src/dev/java/japsadev/bio/hts/metagenome/Path.java
@@ -0,0 +1,167 @@
+package japsadev.bio.hts.metagenome;
+
+import java.util.ArrayList;
+import japsa.seq.Alphabet;
+import japsa.seq.Sequence;
+import japsa.seq.SequenceBuilder;
+
+public class Path implements Comparable{
+ public static int kmer=127;
+ ArrayList nodes;
+ //Graph graph;
+ int length, deviation; //how this path differ to long read data (todo: by multiple-alignment??)
+ public Path(){
+ this.nodes=new ArrayList();
+ //graph=new Graph();
+ length=0;
+ deviation=Integer.MAX_VALUE;
+ }
+
+// public Path(Graph graph){
+// this();
+// associate(graph);
+// }
+ public static void setK(int kmer){
+ Path.kmer=kmer;
+ }
+ public Path(Path p){
+ this();
+ //this(p.graph);
+ for(Node node:p.nodes)
+ this.nodes.add(node);
+ this.length=p.length;
+ }
+ /*
+ * @param String: a path as in contigs.paths of SPAdes output
+ * For example: 1+,2-,3+
+ */
+ public Path(Graph graph, String paths){
+ this();
+ //this(graph);
+ paths=paths.replace(";", ""); //optimized it!
+ String[] comps = paths.split(",");
+ for(int i=0; i getNodes(){
+ return nodes;
+ }
+ public void setComp(ArrayList nodes){
+ this.length=0;
+ for(Node node:nodes)
+ this.addNode(node);
+ }
+
+ public Path rc(){
+ //Path retval=new Path(graph);
+ Path retval=new Path();
+ for(Node node:nodes){
+ retval.nodes.add(0, node.getRC());
+ }
+ return retval;
+ }
+
+ public String toString(){
+ return "P"+getID();
+ }
+ public String getID(){
+ String retval="";
+ for(Node node:nodes){
+ retval+=node.toString();
+ }
+ return retval.trim();
+ }
+ public Node removeLast(){
+ Node retval=nodes.remove(nodes.size()-1);
+ length-=retval.getSeq().length()-kmer;
+ return retval;
+ }
+
+ public Sequence spelling(){
+ SequenceBuilder seq = new SequenceBuilder(Alphabet.DNA16(), 1024*1024, this.toString());
+
+ for(int i=0;i neighborhood;
+ private String fullName, label;
+ private double coverage;
+ private Sequence seq=null;
+ //a sub-graph (path) is equivalent to a Vertex recursively
+ private Path components=null;
+ public Vertex(){
+ fullName="";
+ label="";
+ coverage=0;
+ this.neighborhood = new ArrayList();
+ this.seq=new Sequence(Alphabet.DNA5(), 0);
+ }
+ /**
+ *
+ * @param name The unique label associated with this Vertex
+ */
+ public Vertex(String name){
+ this.fullName=name;
+ this.label=extractID(name);
+ this.coverage=extractCoverage(name);
+ this.neighborhood = new ArrayList();
+ this.seq=new Sequence(Alphabet.DNA5(), 0);
+ }
+
+ public Vertex(String name, Sequence seq){
+ this(name);
+ this.seq = seq;
+ }
+
+ public Vertex(Path p){
+ this();
+ fullName=label=p.getID();
+ components=p;
+ coverage=p.averageCov();
+ }
+ /**
+ *
+ * @return If this vertex is unique (only based on its degree)
+ */
+ public boolean isUnique(){
+ return (getNeighborCount() <= 2);
+ }
+ public Path getSubComps(){
+ return components;
+ }
+ /**
+ * Extract ID from name: EDGE_xx_length_yy_cov_zz;
+ * @param name The name of Edge in assembly graph that correspond to this Vertex
+ */
+ private String extractID(String name){
+ return name.split("_")[1];
+ }
+ /**
+ * Extract coverage from name: EDGE_xx_length_yy_cov_zz;
+ * @param name The name of Edge in assembly graph that correspond to this Vertex
+ * @return coverage value
+ */
+ private double extractCoverage(String name){
+ double res=1.0;
+ try{
+ res=Double.parseDouble(name.split("_")[5]);
+ }catch(Exception e){
+ e.printStackTrace();
+ }
+ return res;
+ }
+ public double getCoverage(){
+ return coverage;
+ }
+ /**
+ * This method adds an Edge to the incidence neighborhood of this graph iff
+ * the edge is not already present.
+ *
+ * @param edge The edge to add
+ */
+ public void addNeighbor(Edge edge){
+ if(this.neighborhood.contains(edge)){
+ return;
+ }
+ this.neighborhood.add(edge);
+ }
+
+
+ /**
+ *
+ * @param other The edge for which to search
+ * @return true iff other is contained in this.neighborhood
+ */
+ public boolean containsNeighbor(Edge other){
+ return this.neighborhood.contains(other);
+ }
+
+ /**
+ *
+ * @param index The index of the Edge to retrieve
+ * @return Edge The Edge at the specified index in this.neighborhood
+ */
+ public Edge getNeighbor(int index){
+ return this.neighborhood.get(index);
+ }
+
+
+ /**
+ *
+ * @param index The index of the edge to remove from this.neighborhood
+ * @return Edge The removed Edge
+ */
+ Edge removeNeighbor(int index){
+ return this.neighborhood.remove(index);
+ }
+
+ /**
+ *
+ * @param e The Edge to remove from this.neighborhood
+ */
+ public void removeNeighbor(Edge e){
+ this.neighborhood.remove(e);
+ }
+
+
+ /**
+ *
+ * @return int The number of neighbors of this Vertex
+ */
+ public int getNeighborCount(){
+ return this.neighborhood.size();
+ }
+ /**
+ *
+ * @return String The label of this Vertex
+ */
+ public String getLabel(){
+ return this.label;
+ }
+ /**
+ *
+ * @return String The full name of this Vertex
+ */
+ public String getName(){
+ return this.fullName;
+ }
+ /**
+ *
+ * @param Sequence A sequence
+ */
+ public void setSequence(Sequence seq){
+ this.seq = seq;
+ }
+ /**
+ *
+ * @return Sequence The sequence of this Vertex
+ */
+ public Sequence getSequence(){
+ return this.seq;
+ }
+ /**
+ *
+ * @return String A String representation of this Vertex
+ */
+ public String toString(){
+ return "Vertex " + label;
+ }
+
+ /**
+ *
+ * @return The hash code of this Vertex's label
+ */
+ public int hashCode(){
+ return this.label.hashCode();
+ }
+
+ /**
+ *
+ * @param other The object to compare
+ * @return true iff other instanceof Vertex and the two Vertex objects have the same label
+ */
+ public boolean equals(Object other){
+ if(!(other instanceof Vertex)){
+ return false;
+ }
+
+ Vertex v = (Vertex)other;
+ return this.label.equals(v.label);
+ }
+
+ /**
+ *
+ * @return ArrayList A copy of this.neighborhood. Modifying the returned
+ * ArrayList will not affect the neighborhood of this Vertex
+ */
+ public ArrayList getNeighbors(){
+ return new ArrayList(this.neighborhood);
+ }
+
+}
+
+
diff --git a/src/dev/java/japsadev/bio/hts/newscarf/Alignment.java b/src/dev/java/japsadev/bio/hts/newscarf/Alignment.java
new file mode 100644
index 0000000..41d04fe
--- /dev/null
+++ b/src/dev/java/japsadev/bio/hts/newscarf/Alignment.java
@@ -0,0 +1,157 @@
+package japsadev.bio.hts.newscarf;
+
+import java.util.ArrayList;
+
+import htsjdk.samtools.Cigar;
+import htsjdk.samtools.CigarElement;
+import htsjdk.samtools.SAMRecord;
+import japsa.seq.Sequence;
+
+public class Alignment implements Comparable {
+ public final static int OVERHANG_THRES=1000;
+ public final static int GOOD_QUAL=60;
+
+ public static int MIN_QUAL=30; //TODO: reduce this by doing self-correction
+
+ int alignLength, quality;
+
+ public String readID;
+ BidirectedNode node;
+
+ public int refStart, refEnd; //1-based position on ref of the start and end of the alignment
+
+ //Position on read of the start and end of the alignment (using the direction of read)
+ //readStart map to refStart, readEnd map to refEnd.
+ //readStart < readEnd if strand = true, else readStart > readEnd
+ public int readStart = 0, readEnd = 0;
+
+ //read length
+ public int readLength = 0;
+
+ public boolean strand = true;//positive
+ public boolean prime = true;//primary alignment
+ public boolean goodMargin = false;
+ public boolean useful = false;
+ //SAMRecord mySam;
+
+ ArrayList alignmentCigars = new ArrayList();
+
+
+ //public int readLeft, readRight, readAlign, refLeft, refRight, refAlign;
+ //left and right are in the direction of the reference sequence
+
+ public Alignment(SAMRecord sam, BidirectedNode node) {
+// readID = Integer.parseInt(sam.getReadName().split("_")[0]);
+ readID = sam.getReadName();
+ quality = sam.getMappingQuality();
+ prime=!sam.getNotPrimaryAlignmentFlag();
+ this.node = node;
+
+ refStart = sam.getAlignmentStart();
+ refEnd = sam.getAlignmentEnd();
+
+ Cigar cigar = sam.getCigar();
+ boolean enterAlignment = false;
+ //////////////////////////////////////////////////////////////////////////////////
+
+ for (final CigarElement e : cigar.getCigarElements()) {
+ alignmentCigars.add(e);
+ final int length = e.getLength();
+ switch (e.getOperator()) {
+ case H :
+ case S :
+ case P : //pad is a kind of clipped
+ if (enterAlignment)
+ readEnd = readLength;
+ readLength += length;
+ break; // soft clip read bases
+ case I :
+ case M :
+ case EQ :
+ case X :
+ if (!enterAlignment){
+ readStart = readLength + 1;
+ enterAlignment = true;
+ }
+ readLength += length;
+ break;
+ case D :
+ case N :
+ if (!enterAlignment){
+ readStart = readLength + 1;
+ enterAlignment = true;
+ }
+ break;
+ default : throw new IllegalStateException("Case statement didn't deal with cigar op: " + e.getOperator());
+ }//case
+ }//for
+ if (readEnd == 0)
+ readEnd = readLength;
+ //these temporary variable to determine usefulness
+ int readLeft = readStart -1;
+ int readRight = readLength - readEnd;
+
+ int refLeft = refStart - 1;
+ int refRight = ((Sequence) node.getAttribute("seq")).length() - refEnd;
+
+ alignLength = refEnd + 1 - refStart;
+ if (sam.getReadNegativeStrandFlag()){
+ strand = false;
+ //need to convert the alignment position on read the correct direction
+ readStart = 1 + readLength - readStart;
+ readEnd = 1 + readLength - readEnd;
+ }
+
+ if (
+ (readLeft < OVERHANG_THRES || refLeft < OVERHANG_THRES) &&
+ (readRight < OVERHANG_THRES || refRight < OVERHANG_THRES)
+ )
+ goodMargin=true;
+
+ if ( goodMargin
+ //prime && //TODO: should be separated as another attribute for further consideration??
+ && alignLength > BidirectedGraph.getKmerSize() //FIXME:
+ && quality >= MIN_QUAL
+ )
+ useful = true;
+
+ }
+
+
+ public int readAlignmentStart(){
+ return Math.min(readStart,readEnd);
+
+ }
+
+ public int readAlignmentEnd(){
+ return Math.max(readStart,readEnd);
+ }
+
+ public String toString() {
+ return node.getAttribute("name")
+ + ": " + refStart
+ + " -> " + refEnd
+ + " / " + ((Sequence) node.getAttribute("seq")).length()
+ + " map to "
+ + readID
+ + ": " + readStart
+ + " -> " + readEnd
+ + " / " + readLength
+ + ", strand: " + (strand?"+":"-")
+ + ", prime: " + (prime?"yes":"no")
+ + ", margin: " + (goodMargin?"good":"bad");
+ }
+ public static boolean isOverlap(Alignment alg1, Alignment alg2){
+ boolean retval=false;
+
+
+ return retval;
+ }
+ /* (non-Javadoc)
+ * @see java.lang.Comparable#compareTo(java.lang.Object)
+ */
+ @Override
+ public int compareTo(Alignment o) {
+ return readAlignmentStart() - o.readAlignmentStart();
+ }
+}
diff --git a/src/dev/java/japsadev/bio/hts/newscarf/BidirectedEdge.java b/src/dev/java/japsadev/bio/hts/newscarf/BidirectedEdge.java
new file mode 100644
index 0000000..0e262e5
--- /dev/null
+++ b/src/dev/java/japsadev/bio/hts/newscarf/BidirectedEdge.java
@@ -0,0 +1,176 @@
+package japsadev.bio.hts.newscarf;
+
+import org.graphstream.graph.Edge;
+import org.graphstream.graph.Node;
+import org.graphstream.graph.implementations.AbstractEdge;
+import org.graphstream.graph.implementations.AbstractNode;
+import org.slf4j.Logger;
+import org.slf4j.LoggerFactory;
+
+import japsa.seq.Sequence;
+
+public class BidirectedEdge extends AbstractEdge{
+ protected boolean dir0, dir1;//true: outward, false: inward
+ private int length=-BidirectedGraph.getKmerSize();//length of the edge (distance between tips of seqs represented by 2 nodes)
+ private BidirectedPath path; //path that represents this edge. Should be a list if metagenome/multiploidy
+
+
+ private static final Logger LOG = LoggerFactory.getLogger(BidirectedEdge.class);
+
+ protected BidirectedEdge(String id, AbstractNode src, AbstractNode dst, boolean dir0, boolean dir1) {
+ // id fuck off!!! we'll make one for ourselves
+ this(src,dst,dir0,dir1);
+ }
+ protected BidirectedEdge(AbstractNode src, AbstractNode dst, boolean dir0, boolean dir1) {
+ //this(createID(src,dst,dir0,dir1), src, dst);
+
+ super(createID(src,dst,dir0,dir1),src,dst,false);
+ this.dir0=dir0;
+ this.dir1=dir1;
+
+ path = new BidirectedPath();
+ path.setRoot(src);
+ path.add(this);
+ }
+
+ /* param id must have the form %s[o/i]%s[o/i], e.g. [1+2-]o[3]i
+ * the constructor will translate the id to the direction property
+ * of the bidirected edge
+ */
+// protected BidirectedEdge(String id, AbstractNode source, AbstractNode dest){
+// super(id, source, dest, false);
+// String pattern = "^\\[([0-9\\+\\-]*)\\]([oi])\\[([0-9\\+\\-]*)\\]([oi])$";
+// // Create a Pattern object
+// Pattern r = Pattern.compile(pattern);
+// // Now create matcher object.
+// Matcher m = r.matcher(id);
+// String leftID, rightID,
+// leftDir, rightDir;
+// if(m.find()){
+// leftID=m.group(1);
+// leftDir=m.group(2);
+// rightID=m.group(3);
+// rightDir=m.group(4);
+// if(source.getId().equals(leftID)){
+// dir0=(leftDir.equals("o")?true:false);
+// dir1=(rightDir.equals("o")?true:false);
+// }else if(source.getId().equals(rightID)){
+// dir0=(rightDir.equals("o")?true:false);
+// dir1=(leftDir.equals("o")?true:false);
+// }else{
+// System.err.println("ID does not match");
+// System.exit(1);
+// }
+// } else{
+// System.err.println("Illegal ID for a bidirected edge (id must have the form node_id[o/i]node_id[o/i])");
+// System.exit(1);
+// }
+//
+// }
+
+ /*
+ * To adapt Graph class from GraphStream library (called by newInstance())
+ */
+ protected BidirectedEdge(String id, AbstractNode source, AbstractNode dest){
+ super(id, source, dest, false);
+
+ assert (source.getGraph() == dest.getGraph()):"Nodes come from different graph " + source.getGraph().getId() + " and " + dest.getGraph().getId();
+ BidirectedGraph g = (BidirectedGraph) source.getGraph();
+ path = new BidirectedPath(g,id);
+ BidirectedNode n0 = (BidirectedNode) path.getRoot(),
+ n1 = (BidirectedNode) path.peekNode();
+
+ BidirectedEdge firstEdge = (BidirectedEdge) path.getEdgePath().get(0),
+ lastEdge = (BidirectedEdge) path.peekEdge();
+ if(n0.getId().equals(source.getId())){
+ dir0 = firstEdge.getDir(n0);
+ dir1 = lastEdge.getDir(n1);
+ }else if(n0.getId().equals(dest.getId())){
+ dir1 = firstEdge.getDir(n0);
+ dir0 = lastEdge.getDir(n1);
+ path = path.getReversedComplemented();
+ }else{
+ LOG.error("Path {} conflicts with src={} dst={}!", id, source.getId(), dest.getId());
+ System.exit(1);
+ }
+ }
+
+// public static String createID(AbstractNode source, AbstractNode dst, boolean dir0, boolean dir1){
+// String srcDes = "["+source.getId()+"]"+(dir0 ? "o":"i"),
+// dstDes = "["+dst.getId()+"]"+(dir1 ? "o":"i");
+// if(srcDes.compareTo(dstDes)<0)
+// return String.format("%s%s", srcDes, dstDes);
+// else
+// return String.format("%s%s", dstDes, srcDes);
+// }
+
+
+ //should have smt like this
+ public static BidirectedEdge makeEdge (AbstractNode src, AbstractNode dst, BidirectedPath path) {
+ return new BidirectedEdge(path.getId(), src, dst);
+ }
+
+// public BidirectedEdge getReversedComplemented(){
+// BidirectedEdge retval = makeEdge(getNode1(),getNode0(),path.getReversedComplemented());
+//
+// return retval;
+// }
+
+ public void setPath(BidirectedPath path){
+ this.path = path;
+ //here set the new length due to that path
+ Node curNode = path.getRoot();
+ for(Edge e:path.getEdgePath()){
+ curNode = e.getOpposite(curNode);
+ if(curNode == path.peekNode())
+ break;
+ else{
+ length+=((Sequence)curNode.getAttribute("seq")).length()-BidirectedGraph.getKmerSize();
+ }
+ }
+ }
+ public BidirectedPath getPath(){
+ return path;
+ }
+
+ @Override
+ public String toString() {
+ return String.format("%s:%s-%s-%s-%s", getId(), source, (dir0?">":"<"), (dir1?"<":">"), target);
+ }
+
+ public void setDir0(boolean dir){
+ this.dir0=dir;
+ }
+ public void setDir1(boolean dir){
+ this.dir1=dir;
+ }
+ public boolean getDir0(){
+ return dir0;
+ }
+ public boolean getDir1(){
+ return dir1;
+ }
+ //if kmer!=127 we need to set initial lengths of edges again (easy+tedious way)
+ public void changeKmerSize(int kmer){
+ length=-kmer;
+ }
+ public int getLength(){
+ return length;
+ }
+
+ //TODO: include the case of tandem repeats
+ public boolean getDir(AbstractNode node){
+ assert node==getSourceNode()||node==getTargetNode():"Node " + node.getId() + " does not belong to this edge src=" + getSourceNode().getId() + " dst=" + getTargetNode().getId();
+ return node==getSourceNode()?getDir0():getDir1();
+ }
+
+
+ public static String createID(AbstractNode source, AbstractNode dst, boolean dir0, boolean dir1){
+ String srcDes = source.getId(),
+ dstDes = dst.getId();
+ if(srcDes.compareTo(dstDes)<0)
+ return String.format("%s%s,%s%s", srcDes, (dir0?"+":"-"), dstDes, (dir1?"-":"+"));
+ else
+ return String.format("%s%s,%s%s", dstDes, (dir1?"+":"-"), srcDes, (dir0?"-":"+"));
+ }
+}
diff --git a/src/dev/java/japsadev/bio/hts/newscarf/BidirectedGraph.java b/src/dev/java/japsadev/bio/hts/newscarf/BidirectedGraph.java
new file mode 100644
index 0000000..a3556d3
--- /dev/null
+++ b/src/dev/java/japsadev/bio/hts/newscarf/BidirectedGraph.java
@@ -0,0 +1,586 @@
+package japsadev.bio.hts.newscarf;
+
+import java.io.IOException;
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.Iterator;
+import java.util.List;
+
+import org.graphstream.graph.*;
+import org.graphstream.graph.implementations.*;
+import org.slf4j.Logger;
+import org.slf4j.LoggerFactory;
+
+import japsa.seq.Alphabet;
+import japsa.seq.FastaReader;
+import japsa.seq.Sequence;
+import japsa.seq.SequenceReader;
+
+
+public class BidirectedGraph extends AdjacencyListGraph{
+ static int kmer=127;
+ static final int TOLERATE=500;
+ static final int D_LIMIT=10000;
+ static final int S_LIMIT=50;
+
+ private static final Logger LOG = LoggerFactory.getLogger(BidirectedGraph.class);
+
+ // *** Constructors ***
+ /**
+ * Creates an empty graph.
+ *
+ * @param id
+ * Unique identifier of the graph.
+ * @param strictChecking
+ * If true any non-fatal error throws an exception.
+ * @param autoCreate
+ * If true (and strict checking is false), nodes are
+ * automatically created when referenced when creating a edge,
+ * even if not yet inserted in the graph.
+ * @param initialNodeCapacity
+ * Initial capacity of the node storage data structures. Use this
+ * if you know the approximate maximum number of nodes of the
+ * graph. The graph can grow beyond this limit, but storage
+ * reallocation is expensive operation.
+ * @param initialEdgeCapacity
+ * Initial capacity of the edge storage data structures. Use this
+ * if you know the approximate maximum number of edges of the
+ * graph. The graph can grow beyond this limit, but storage
+ * reallocation is expensive operation.
+ */
+ public BidirectedGraph(String id, boolean strictChecking, boolean autoCreate,
+ int initialNodeCapacity, int initialEdgeCapacity) {
+ super(id, strictChecking, autoCreate);
+ // All we need to do is to change the node & edge factory
+ setNodeFactory(new NodeFactory() {
+ public BidirectedNode newInstance(String id, Graph graph) {
+ return new BidirectedNode((AbstractGraph) graph, id);
+ }
+ });
+
+ setEdgeFactory(new EdgeFactory() {
+ public BidirectedEdge newInstance(String id, Node src, Node dst, boolean directed) { //stupid??
+ return new BidirectedEdge(id, (AbstractNode)src, (AbstractNode)dst);
+ }
+ });
+
+ }
+
+ /**
+ * Creates an empty graph with default edge and node capacity.
+ *
+ * @param id
+ * Unique identifier of the graph.
+ * @param strictChecking
+ * If true any non-fatal error throws an exception.
+ * @param autoCreate
+ * If true (and strict checking is false), nodes are
+ * automatically created when referenced when creating a edge,
+ * even if not yet inserted in the graph.
+ */
+ public BidirectedGraph(String id, boolean strictChecking, boolean autoCreate) {
+ this(id, strictChecking, autoCreate, DEFAULT_NODE_CAPACITY,
+ DEFAULT_EDGE_CAPACITY);
+ }
+
+ /**
+ * Creates an empty graph with strict checking and without auto-creation.
+ *
+ * @param id
+ * Unique identifier of the graph.
+ */
+ public BidirectedGraph(String id) {
+ this(id, true, false, 10000, 100000);
+ }
+
+ public BidirectedGraph(){
+ this("Assembly graph",true,false, 10000, 100000);
+ setKmerSize(127);//default kmer size used by SPAdes to assembly MiSeq data
+ }
+
+ //just to make AbstractGraph.removeEdge(AbstractEdge, boolean, boolean, boolean) visible
+ protected void removeEdgeDup(AbstractEdge edge, boolean graphCallback,
+ boolean sourceCallback, boolean targetCallback) {
+ this.removeEdge(edge, graphCallback, sourceCallback, targetCallback);
+
+ }
+ protected BidirectedEdge addEdge(AbstractNode src, AbstractNode dst, boolean dir0, boolean dir1){
+ BidirectedEdge tmp = addEdge(BidirectedEdge.createID(src, dst, dir0, dir1), src, dst);
+// String s1=tmp.toString();
+// //tmp.setDir0(dir0);
+// //tmp.setDir1(dir1);
+// if(!s1.equals(tmp.toString()))
+// System.out.println(s1 + " ---> " + tmp);
+ return tmp;
+ }
+
+ public String printEdgesOfNode(BidirectedNode node){
+ Iterator ins = getNode(node.getId()).getEnteringEdgeIterator(),
+ outs = getNode(node.getId()).getLeavingEdgeIterator();
+ String retval=node.getId() + ": IN={";
+ while(ins.hasNext())
+ retval += ins.next().getId() + " ";
+ retval+="}; OUT={";
+ while(outs.hasNext())
+ retval += outs.next().getId() + " ";
+ retval+="}";
+ return retval;
+ }
+ /**********************************************************************************
+ * ****************************Algorithms go from here*****************************
+ */
+ //TODO: read from ABySS assembly graph (graph of final contigs, not like SPAdes)
+ private static double aveCov; //TODO: replaced with more accurate method
+
+ public void loadFromFile(String graphFile) throws IOException{
+ setAutoCreate(true);
+ setStrict(false);
+ //1. next iterate over again to read the connections
+ SequenceReader reader = new FastaReader(graphFile);
+ Sequence seq;
+ int shortestLen = 10000;
+ int totReadLen=0, totGenomeLen=0;
+ while ((seq = reader.nextSequence(Alphabet.DNA())) != null){
+ if(seq.length() 1){
+ String[] nbList = adjList[1].split(",");
+ for(int i=0; i < nbList.length; i++){
+ String neighbor = nbList[i];
+ // note that the direction is read reversely in the dest node
+ boolean dir1=neighbor.contains("'")?true:false;
+ neighbor=neighbor.replaceAll("[^a-zA-Z0-9_.]", "").trim();
+
+ String neighborID = neighbor.split("_")[1];
+ AbstractNode nbr = addNode(neighborID);
+ nbr.setAttribute("cov", Double.parseDouble(neighbor.split("_")[5]));
+
+ addEdge(node, nbr, dir0, dir1);
+ //e.addAttribute("ui.label", e.getId());
+ }
+ }
+
+ }
+
+ //rough estimation of kmer used
+ if((shortestLen-1) != getKmerSize()){
+ setKmerSize(shortestLen-1);
+ for(Edge e:getEdgeSet()){
+ ((BidirectedEdge)e).changeKmerSize(kmer);
+ }
+ }
+
+ aveCov = totReadLen/totGenomeLen;
+ reader.close();
+
+ }
+
+
+ public static int getKmerSize(){
+ return BidirectedGraph.kmer;
+ }
+ public static void setKmerSize(int kmer){
+ BidirectedGraph.kmer=kmer;
+ }
+
+
+ /*
+ * Read paths from contigs.path and reduce the graph
+ */
+// public void readPathsFromSpades(String paths) throws IOException{
+//
+// BufferedReader pathReader = new BufferedReader(new FileReader(paths));
+//
+// String s;
+// //Read contigs from contigs.paths and refer themselves to contigs.fasta
+// boolean flag=false;
+// while((s=pathReader.readLine()) != null){
+// if(s.contains("NODE")){
+// flag=s.contains("'")?false:true;
+// continue;
+// }else if(flag){
+// BidirectedPath path=new BidirectedPath(this, s);
+//// System.out.println("Using path to reduce: " + path.getId());
+//// System.out.println("Before reduce => Node: " + getNodeCount() + " Edge: " + getEdgeCount());
+//
+//// AbstractNode comp=
+//
+// this.reduce(path);
+//
+//
+//// if(comp!=null){
+//// System.out.println("Reverting node: " + comp.getId());
+//// revert(comp);
+//// System.out.println("After revert => Node: " + getNodeCount() + " Edge: " + getEdgeCount());
+////
+//// }
+// }
+//
+//
+// }
+// pathReader.close();
+// }
+
+// /**
+// *
+// * @param p Path to be grouped as a virtually vertex
+// */
+// public AbstractNode reduce(BidirectedPath p){
+// //do nothing if the path has only one node
+// if(p==null || p.getEdgeCount()<1)
+// return null;
+//
+// //now only work with path containing more than 2 unique nodes
+// int uniqueCount=0;
+// for(Node n:p.getEachNode()){
+// if(isUnique(n))
+// uniqueCount++;
+// }
+// if(uniqueCount < 2)
+// {
+// System.out.println("ignore path with less than 1 unique contig!");
+// return null;
+// }
+// //add the new composite Node to the graph
+// //compare id from sense & anti-sense to get the unique one
+// AbstractNode comp = addNode(p.getId().compareTo(p.getReversedComplemented().getId())>0?
+// p.getReversedComplemented().getId():p.getId());
+//
+// comp.addAttribute("path", p);
+// comp.addAttribute("seq", p.spelling());
+// comp.addAttribute("ui.label", comp.getId());
+// comp.setAttribute("ui.style", "text-offset: -10;");
+// comp.setAttribute("ui.class", "marked");
+// try { Thread.sleep(100); } catch (Exception e) {}
+//
+// //store unique nodes on p for removing
+// ArrayList tobeRemoved=new ArrayList();
+// for(Node n:p.getEachNode()){
+// if(isUnique(n))
+// tobeRemoved.add(n.getId());
+// }
+// BidirectedNode start = (BidirectedNode) p.getRoot(),
+// end = (BidirectedNode) p.peekNode();
+// boolean startDir = ((BidirectedEdge) p.getEdgePath().get(0)).getDir(start),
+// endDir = ((BidirectedEdge) p.peekEdge()).getDir(end);
+// //set neighbors of the composite Node
+// Iterator startEdges = startDir?start.getEnteringEdgeIterator():start.getLeavingEdgeIterator(),
+// endEdges = endDir?end.getEnteringEdgeIterator():end.getLeavingEdgeIterator();
+// while(startEdges.hasNext()){
+// BidirectedEdge e = (BidirectedEdge) startEdges.next();
+// BidirectedNode opNode = e.getOpposite(start);
+// boolean opDir = e.getDir(opNode);
+// //Edge tmp=
+// addEdge(BidirectedEdge.createID(comp, opNode, false, opDir), comp, opNode);//always into start node
+// //System.out.println("From " + start.getId() + ": " + tmp.getId() + " added!");
+// }
+//
+// while(endEdges.hasNext()){
+// BidirectedEdge e = (BidirectedEdge) endEdges.next();
+// BidirectedNode opNode = e.getOpposite(end);
+// boolean opDir = e.getDir(opNode);
+// //Edge tmp=
+// addEdge(BidirectedEdge.createID(comp, opNode, true, opDir), comp, opNode);//always out of end node
+//
+// //System.out.println("From " + end.getId() + ": " + tmp.getId() + " added!");
+//
+// }
+//
+// for(String nLabel:tobeRemoved){
+// //System.out.println("About to remove " + nLabel);
+// removeNode(nLabel);
+// }
+//
+// //TODO: remove bubbles...
+// return comp;
+// }
+//
+//
+// /**
+// *
+// * @param v Node to be reverted (1-level reverting)
+// */
+// public void revert(AbstractNode v){
+// System.out.println("Reverting...");
+// Path p=v.getAttribute("path");
+// if(p==null) return;
+//
+// BidirectedNode start = (BidirectedNode) p.getRoot(),
+// end = (BidirectedNode) p.peekNode();
+// boolean startDir = ((BidirectedEdge) p.getEdgePath().get(0)).getDir(start),
+// endDir = ((BidirectedEdge) p.peekEdge()).getDir(end);
+//
+// //add back all neighbor edges of this composite vertex
+// Iterator startEdges = v.getEnteringEdgeIterator(),
+// endEdges = v.getLeavingEdgeIterator();
+// //add back all nodes from the path
+// for(Node n:p.getNodeSet()){
+// if(getNode(n.getId())!=null)
+// continue;
+// Node tmp = addNode(n.getId());
+// tmp.addAttribute("seq", (japsa.seq.Sequence)n.getAttribute("seq"));
+// tmp.addAttribute("name", (String)n.getAttribute("name"));
+// tmp.addAttribute("path", (BidirectedPath)n.getAttribute("path"));
+//
+// //System.out.println("Adding back edge "+tmp.getId());
+// }
+// while(startEdges.hasNext()){
+// BidirectedEdge e = (BidirectedEdge) startEdges.next();
+// BidirectedNode opNode = e.getOpposite(v);
+// boolean opDir = e.getDir(opNode);
+// //Edge tmp =
+// addEdge(BidirectedEdge.createID(start, opNode, !startDir, opDir), start, opNode);
+// //System.out.println("Adding back edge "+tmp.getId());
+// }
+//
+// while(endEdges.hasNext()){
+// BidirectedEdge e = (BidirectedEdge) endEdges.next();
+// BidirectedNode opNode = e.getOpposite(v);
+// boolean opDir = e.getDir(opNode);
+// //Edge tmp =
+// addEdge(BidirectedEdge.createID(end, opNode, !endDir, opDir), end, opNode);
+// //System.out.println("Adding back edge "+tmp.getId());
+// }
+//
+// //add back all edges from the path
+// for(Edge e:p.getEdgeSet()){
+// //Edge tmp =
+// addEdge(e.getId(), e.getSourceNode().getId(), e.getTargetNode().getId());
+// //System.out.println("Adding back edge "+tmp.getId());
+// }
+// //finally remove the composite node
+// removeNode(v);
+// }
+
+ /*
+ * This function deduces a full path in this graph between 2 nodes aligned with a long read
+ */
+ protected ArrayList getClosestPath(Alignment from, Alignment to, int distance){
+ BidirectedNode srcNode = from.node,
+ dstNode = to.node;
+ System.out.println("Looking for path between " + srcNode.getId() + " to " + dstNode.getId() + " with distance " + distance);
+ BidirectedPath tmp = new BidirectedPath();
+ ArrayList possiblePaths = new ArrayList(),
+ retval = new ArrayList();
+ tmp.setRoot(srcNode);
+
+ //traverse(tmp, dest, retval, distance+source.getSeq().length()+dest.getSeq().length());
+ traverse(tmp, dstNode, possiblePaths, distance, from.strand, to.strand, 0);
+ //only get the best ones
+ if(possiblePaths.isEmpty()){
+ //if a path couldn't be found between 2 dead-ends but alignments quality are insane high
+ //FIXME: return a pseudo path having an nanopore edge
+ if(isUnique(srcNode) && isUnique(dstNode) && srcNode.getDegree() == 1 && dstNode.getDegree()==1 &&
+ Math.min(from.quality, to.quality) >= Alignment.GOOD_QUAL)
+ {
+ BidirectedEdge pseudoEdge = new BidirectedEdge(srcNode, dstNode, from.strand, to.strand);
+ //TODO: save the corresponding content of long reads to this edge
+ pseudoEdge.setAttribute("pseudo", distance);
+ tmp.add(pseudoEdge);
+ retval.add(tmp);
+ System.out.println("pseudo path from " + srcNode.getId() + " to " + dstNode.getId());
+// HybridAssembler.promptEnterKey();
+ return retval;
+ }else
+ return null;
+ }
+ double bestScore=possiblePaths.get(0).getDeviation();
+ for(int i=0;i curResult,
+ int distance, boolean srcDir, boolean dstDir, int stepCount)
+ {
+ //stop if it's going too far!
+ if(stepCount >= S_LIMIT)
+ return;
+
+ BidirectedNode currentNode=(BidirectedNode) path.peekNode();
+ BidirectedEdge currentEdge;
+ boolean curDir;//direction to the next node, = ! previous'
+
+ Iterator ite;
+ if(path.size() <= 1) //only root
+ curDir=srcDir;//re-check
+ else{
+ currentEdge = (BidirectedEdge) path.peekEdge();
+ curDir = !((BidirectedEdge) currentEdge).getDir(currentNode);
+ }
+ ite=curDir?currentNode.getLeavingEdgeIterator():currentNode.getEnteringEdgeIterator();
+
+ while(ite.hasNext()){
+ BidirectedEdge e = ite.next();
+ path.add(e);
+
+ int toTarget=Math.abs(distance-e.getLength());
+ if(e.getOpposite(currentNode)==dst && e.getDir(dst)!=dstDir && toTarget < TOLERATE){
+ BidirectedPath curPath=curResult.isEmpty()?new BidirectedPath():curResult.get(0), //the best path saved among all possible paths from the list curResult
+ tmpPath=new BidirectedPath(path);
+ tmpPath.setDeviation(toTarget);
+ if( toTarget < curPath.getDeviation() )
+ curResult.add(0, tmpPath);
+ else
+ curResult.add(tmpPath);
+
+ System.out.println("Hit added: "+path.getId()+"(candidate deviation: "+toTarget+")");
+ }else{
+ int newDistance = distance - ((Sequence) e.getOpposite(currentNode).getAttribute("seq")).length() - e.getLength();
+// System.out.println("adding edge: " + e.getId() + " length=" + e.getLength() +" -> distance=" + newDistance);
+ if (newDistance - e.getLength() < -TOLERATE){
+ System.out.println("Stop go to edge " + e.getPath() + " from path with distance "+newDistance+" already! : "+path.getId());
+ }else
+ traverse(path, dst, curResult, newDistance, srcDir, dstDir, stepCount++);
+ }
+ path.popNode();
+
+ }
+ }
+
+ /*
+ * Find a path based on list of Alignments
+ */
+ public BidirectedPath pathFinding(ArrayList alignments) {
+ if(alignments.size()<=1)
+ return null;
+
+ System.out.println("=================================================");
+ for(Alignment alg:alignments)
+ System.out.println("\t"+alg.toString());
+ System.out.println("=================================================");
+ //First bin the alignments into different overlap regions
+ //only considering useful alignments
+ HashMap allAlignments = new HashMap();
+ ArrayList joinPaths = new ArrayList();
+
+ for(Alignment alg:alignments){
+ if(alg.useful){
+ Range range = new Range(alg.readAlignmentStart(),alg.readAlignmentEnd());
+ allAlignments.put(range, alg);
+ }
+ }
+ //now get all the bin in order
+ List baseRanges=new ArrayList(allAlignments.keySet());
+ List> rangeGroups = MetaRange.getOverlappingGroups(baseRanges);
+
+ if(rangeGroups.size() < 2)
+ return null;
+
+ System.out.println("Binning ranges: ");
+ for(List group : rangeGroups){
+ System.out.println(group);
+ }
+
+ //iterate all alignments in adjacent bins to find correct path
+ ArrayList curGroup = new ArrayList(),
+ nextGroup = new ArrayList();
+ List curRanges = rangeGroups.get(0);
+ for(Range r:curRanges)
+ curGroup.add(allAlignments.get(r));
+
+ for(int i=1; i nextRanges = rangeGroups.get(i);
+ nextGroup = new ArrayList();
+
+ for(Range r:nextRanges)
+ nextGroup.add(allAlignments.get(r));
+
+ ArrayList allPaths = new ArrayList();
+
+
+ for(Alignment curAlg:curGroup){
+ for(Alignment nextAlg:nextGroup){
+ int distance = nextAlg.readAlignmentStart()-curAlg.readAlignmentEnd();
+ if(distance bridges = getClosestPath(curAlg, nextAlg, distance);
+ if(bridges!=null)
+ allPaths.addAll(bridges);
+ }
+ }
+ }
+ //join all paths from previous to the new ones
+ //TODO:optimize it
+ if(joinPaths.isEmpty())
+ joinPaths=allPaths;
+ else{
+ System.out.println("=====Current list of paths: " + joinPaths);
+ System.out.println("=====Join to list of paths: " + allPaths);
+
+ for(BidirectedPath p:joinPaths)
+ for(BidirectedPath e:allPaths)
+ p.join(e);
+
+ }
+ curGroup=nextGroup;
+ }
+
+ for(BidirectedPath path:joinPaths)
+ System.out.println("A member Path: " + path.toString() + " deviation: " + path.getDeviation());
+ if(joinPaths.isEmpty())
+ return null;
+ else
+ return joinPaths.get(0);
+ }
+
+ /*
+ * Important function: determine if a node is able to be removed or not
+ * TODO: re-implement it based on statistics of coverage also
+ * 1. pick the least coverage ones among a path as the base
+ * 2. global base
+ */
+ public static boolean isUnique(Node node){
+ boolean res = false;
+
+ if(node.getDegree()<=2){ // not always true, e.g. unique node in a repetitive component
+ Sequence seq = node.getAttribute("seq");
+ if(seq.length() > 10000 || node.getNumber("cov")/aveCov < 1.3)
+ res=true;
+ }
+
+// if(res)
+// LOG.info(node.getAttribute("name") + " with coverage " + node.getNumber("cov") + " is a marker!");
+// else
+// LOG.info(node.getAttribute("name") + " with coverage " + node.getNumber("cov") + " is NOT a marker!");
+
+ return res;
+ }
+ /* need more powerful function:
+ * A-statistics?
+ * Mixture of Poisson distributions??? Kalman filter idea...
+ */
+// public static boolean isUnique(Node node, double cov){
+// boolean res = false;
+// if(node.getDegree()<=2 || Math.abs(node.getAttribute("cov" )) < cov){
+//// if(((Sequence)node.getAttribute("seq")).length() > 5000 || node.getDegree()==0)
+// res=true;
+// }
+//
+// return res;
+// }
+}
diff --git a/src/dev/java/japsadev/bio/hts/newscarf/BidirectedNode.java b/src/dev/java/japsadev/bio/hts/newscarf/BidirectedNode.java
new file mode 100644
index 0000000..0da20a4
--- /dev/null
+++ b/src/dev/java/japsadev/bio/hts/newscarf/BidirectedNode.java
@@ -0,0 +1,301 @@
+package japsadev.bio.hts.newscarf;
+
+import java.security.AccessControlException;
+import java.util.Arrays;
+import java.util.HashMap;
+import java.util.Iterator;
+import java.util.LinkedList;
+import java.util.List;
+import java.util.NoSuchElementException;
+
+import org.graphstream.graph.Edge;
+import org.graphstream.graph.Node;
+import org.graphstream.graph.implementations.*;
+import org.slf4j.Logger;
+import org.slf4j.LoggerFactory;
+/**
+ * Similar to {@link AdjacencyListNode}
+ *
+ */
+public class BidirectedNode extends AbstractNode {
+ protected int COPY=1;
+
+ protected static final int INITIAL_EDGE_CAPACITY;
+ protected static final double GROWTH_FACTOR = 1.1;
+
+ private static final Logger LOG = LoggerFactory.getLogger(BidirectedNode.class);
+
+ static {
+ String p = "org.graphstream.graph.node.initialEdgeCapacity";
+ int initialEdgeCapacity = 32;
+ try {
+ initialEdgeCapacity = Integer.valueOf(System.getProperty(p, "32"));
+ } catch (AccessControlException e) {
+ }
+ INITIAL_EDGE_CAPACITY = initialEdgeCapacity;
+ }
+ //edges are bidirected, here are 4 sub-types (name based on direction of the arrow relative to the corresponding node):
+ //note that neighbor edges here will treat their root node as *left* node, the opposite node as *right* node
+ protected static final byte OO_EDGE = 0b00; // src-->--<--dst
+ protected static final byte OI_EDGE = 0b01; // src-->-->--dst
+ protected static final byte IO_EDGE = 0b10; // src--<--<--dst
+ protected static final byte II_EDGE = 0b11; // src--<-->--dst
+
+ protected BidirectedEdge[] edges;//fast access to edges knowing direction from src
+ protected int oStart, degree;
+
+ protected HashMap> neighborMap; //fast access to edges knowing dst
+ // *** Constructor ***
+
+ protected BidirectedNode(AbstractGraph graph, String id) {
+ super(graph, id);
+ edges = new BidirectedEdge[INITIAL_EDGE_CAPACITY];
+ oStart = degree = 0;
+ neighborMap = new HashMap>(
+ 4 * INITIAL_EDGE_CAPACITY / 3 + 1);
+ }
+
+ // *** Helpers ***
+
+ protected byte edgeType(BidirectedEdge e) {
+ //return (byte) (e.getDir((AbstractNode) e.getOpposite(this))?0:1 + (e.getDir(this)?0:1)<<1); //cool but less efficient
+ BidirectedNode opposite = e.getOpposite(this);
+ if(e.getDir(this))
+ if(e.getDir(opposite))
+ return OO_EDGE;
+ else
+ return OI_EDGE;
+ else
+ if(e.getDir(opposite))
+ return IO_EDGE;
+ else
+ return II_EDGE;
+ }
+
+ @SuppressWarnings("unchecked")
+ protected BidirectedEdge locateEdge(Node opposite, byte type) {
+ List l = neighborMap.get(opposite);
+ if (l == null)
+ return null;
+
+ for (BidirectedEdge e : l) {
+ if(type==edgeType(e))
+ return e;
+ }
+ return null;
+ }
+
+ protected void removeEdge(int i) {
+// System.out.print("From node " + this.getId() + " remove edge number " + i);
+// System.out.println(" from a list of edges: ");
+// for(int j=0;j l = neighborMap.get(opposite);
+ l.remove(edges[i]);
+ if (l.isEmpty())
+ neighborMap.remove(opposite);
+ //remove from the array
+ if (i >= oStart) {
+ edges[i] = edges[--degree];
+ edges[degree] = null;
+ return;
+ }
+
+ edges[i] = edges[--oStart];
+ edges[oStart] = edges[--degree];
+ edges[degree] = null;
+
+ }
+
+ // *** Callbacks ***
+
+ @Override
+ protected boolean addEdgeCallback(AbstractEdge edge) {
+ //LOG.info("Adding edge callback " + edge.getId() + " from graph " + getGraph().getId());
+ AbstractNode opposite = edge.getOpposite(this);
+ List l = neighborMap.get(opposite);
+ if (l == null) {
+ l = new LinkedList();
+ neighborMap.put(opposite, l);
+ }
+ l.add((BidirectedEdge) edge);
+
+ // resize edges if necessary
+ if (edges.length == degree) {
+ BidirectedEdge[] tmp = new BidirectedEdge[(int) (GROWTH_FACTOR * edges.length) + 1];
+ System.arraycopy(edges, 0, tmp, 0, edges.length);
+ Arrays.fill(edges, null);
+ edges = tmp;
+ }
+
+ byte type = edgeType((BidirectedEdge) edge);
+
+ if (type <= OI_EDGE) {
+ edges[degree++] = (BidirectedEdge) edge;
+ return true;
+ }
+
+ edges[degree++] = edges[oStart];
+ edges[oStart++] = (BidirectedEdge) edge;
+ return true;
+ }
+
+ @Override
+ protected void removeEdgeCallback(AbstractEdge edge) {
+ //LOG.info("Removing edge callback " + edge.getId() + " from graph " + getGraph().getId());
+
+ // locate the edge first
+ byte type = edgeType((BidirectedEdge) edge);
+ int i = 0;
+ if (type <= OI_EDGE)
+ i = oStart;
+ while (i <= degree && edges[i] != edge)
+ i++;
+ if(i < degree){ //only remove iff edge is found
+ removeEdge(i);
+ }
+ }
+
+ @Override
+ protected void clearCallback() {
+ Arrays.fill(edges, 0, degree, null);
+ oStart = degree = 0;
+ }
+
+ // *** Access methods ***
+
+ @Override
+ public int getDegree() {
+ return degree;
+ }
+
+ @Override
+ public int getInDegree() {
+ return oStart;
+ }
+
+ @Override
+ public int getOutDegree() {
+ return degree - oStart;
+ }
+
+ @SuppressWarnings("unchecked")
+ @Override
+ public T getEdge(int i) {
+ if (i < 0 || i >= degree)
+ throw new IndexOutOfBoundsException("Node \"" + this + "\""
+ + " has no edge " + i);
+ return (T) edges[i];
+ }
+
+ @SuppressWarnings("unchecked")
+ @Override
+ public T getEnteringEdge(int i) {
+ if (i < 0 || i >= getInDegree())
+ throw new IndexOutOfBoundsException("Node \"" + this + "\""
+ + " has no entering edge " + i);
+ return (T) edges[i];
+ }
+
+ @SuppressWarnings("unchecked")
+ @Override
+ public T getLeavingEdge(int i) {
+ if (i < 0 || i >= getOutDegree())
+ throw new IndexOutOfBoundsException("Node \"" + this + "\""
+ + " has no edge " + i);
+ return (T) edges[oStart + i];
+ }
+
+ // FIXME: I must override these stupid functions, let's just return random edge among 4 types!!!
+ @SuppressWarnings("unchecked")
+ @Override
+ public T getEdgeBetween(Node node) {
+ return (T) locateEdge(node, IO_EDGE);
+ }
+
+ @SuppressWarnings("unchecked")
+ @Override
+ public T getEdgeFrom(Node node) {
+ return (T) locateEdge(node, OO_EDGE);
+ }
+
+ @SuppressWarnings("unchecked")
+ @Override
+ public T getEdgeToward(Node node) {
+ return (T) locateEdge(node, II_EDGE);
+ }
+
+ // *** Iterators ***
+
+ protected class EdgeIterator implements Iterator {
+ protected int iPrev, iNext, iEnd;
+ //0:in, 1:out, other(2):all
+ protected EdgeIterator(int ori) {
+ iPrev = -1;
+ iNext = 0;
+ iEnd = degree;
+ if (ori==0)
+ iEnd = oStart;
+ else if(ori==1)
+ iNext = oStart;
+
+// System.out.println("Iterator " + ori + " of " + getId() + " from " + iNext + " to " + iEnd + " of");
+// for(int i=0;i= iEnd)
+ throw new NoSuchElementException();
+ iPrev = iNext++;
+ return (T) edges[iPrev];
+ }
+
+ public void remove() {
+ if (iPrev == -1)
+ throw new IllegalStateException();
+ AbstractEdge e = edges[iPrev];
+ // do not call the callback because we already know the index
+ //graph.removeEdge(e);
+ ((BidirectedGraph)graph).removeEdgeDup(e, true, e.getSourceNode() != BidirectedNode.this,
+ e.getTargetNode() != BidirectedNode.this);
+ removeEdge(iPrev);
+ iNext = iPrev;
+ iPrev = -1;
+ iEnd--;
+
+ }
+ }
+
+ @Override
+ public Iterator getEdgeIterator() {
+ return new EdgeIterator(2);
+ }
+
+ @Override
+ public Iterator getEnteringEdgeIterator() {
+ return new EdgeIterator(0);
+ }
+
+ @Override
+ public Iterator getLeavingEdgeIterator() {
+ return new EdgeIterator(1);
+ }
+
+}
diff --git a/src/dev/java/japsadev/bio/hts/newscarf/BidirectedPath.java b/src/dev/java/japsadev/bio/hts/newscarf/BidirectedPath.java
new file mode 100644
index 0000000..0ba5f55
--- /dev/null
+++ b/src/dev/java/japsadev/bio/hts/newscarf/BidirectedPath.java
@@ -0,0 +1,191 @@
+package japsadev.bio.hts.newscarf;
+
+import japsa.seq.Alphabet;
+import japsa.seq.Sequence;
+import japsa.seq.SequenceBuilder;
+
+
+import java.util.List;
+
+import org.graphstream.graph.Edge;
+import org.graphstream.graph.Node;
+import org.graphstream.graph.Path;
+import org.graphstream.graph.implementations.AbstractNode;
+import org.slf4j.Logger;
+import org.slf4j.LoggerFactory;
+
+public class BidirectedPath extends Path{
+ int deviation; //how this path differ to long read data (todo: by multiple-alignment??)
+ private double coverage=-1; //representative coverage of this path (basically the lowest cov from its nodes)
+ private static final Logger LOG = LoggerFactory.getLogger(BidirectedPath.class);
+
+ @Override
+ public void add(Edge edge) {
+ super.add(edge);
+
+ double newCoverage = peekNode().getNumber("cov");
+ if(coverage<=0)
+ coverage=newCoverage;
+ else
+ coverage=Math.min(coverage, newCoverage);
+
+ }
+
+ public BidirectedPath(){
+ super();
+ }
+ public BidirectedPath(BidirectedPath p){
+ super();
+ if(p!=null && !p.empty()){
+ setRoot(p.getRoot());
+ for(Edge e:p.getEdgePath())
+ add(e);
+ }
+ deviation=p.deviation;
+ }
+
+ //This constructor is only used to load in contigs.path from SPAdes
+ //So no recursive path here (path contains all primitive edges)
+ public BidirectedPath(BidirectedGraph graph, String paths){
+ super();
+ paths=paths.replace(";", ","); //optimized it!
+ String[] comps = paths.split(",");
+ if(comps.length<1)
+ return;
+ String curID = comps[0], nextID;
+ boolean curDir = curID.contains("+")?true:false,
+ nextDir;
+ BidirectedNode curNode = graph.getNode(curID.substring(0,curID.length()-1)),
+ nextNode;
+ setRoot(curNode);
+ for(int i=1; i edges = this.getEdgePath();
+ for(int i = edges.size()-1; i>=0; i--)
+ rcPath.add(edges.get(i));
+ return rcPath;
+ }
+ //It is not really ID because Path doesn't need an ID
+ public String getId(){
+ //need to make the Id unique for both sense and antisense spelling???
+ BidirectedNode curNode = (BidirectedNode) getRoot();
+ if(getEdgeCount()<1)
+ return curNode.getId();
+
+ String retval=curNode.getId(),
+ curDir=((BidirectedEdge) getEdgePath().get(0)).getDir(curNode)?"+":"-";
+ retval+=curDir;
+ for(Edge e:getEdgePath()){
+ curNode=e.getOpposite(curNode);
+ retval+=","+curNode.getId();
+ curDir=((BidirectedEdge) e).getDir(curNode)?"-":"+"; //note that curNode is target node
+ retval+=curDir;
+ }
+
+ return retval.trim();
+ }
+
+ public String toString(){
+ return "(" + getId() + ")";
+ }
+
+ public Sequence spelling(){
+
+ BidirectedNode curNode = (BidirectedNode) getRoot();
+ Sequence curSeq = curNode.getAttribute("seq");
+
+ if(getEdgeCount()<1)
+ return curSeq;
+
+ SequenceBuilder seq = new SequenceBuilder(Alphabet.DNA16(), 1024*1024, this.toString());
+ boolean curDir=((BidirectedEdge) getEdgePath().get(0)).getDir(curNode);
+ curSeq = curDir?curSeq:Alphabet.DNA.complement(curSeq);
+
+ seq.append(curSeq.subSequence(0, curSeq.length()-BidirectedGraph.getKmerSize()));
+ for(Edge edge:getEdgePath()){
+ for(Edge e:((BidirectedEdge) edge).getPath().getEdgePath()){
+ curNode=e.getOpposite(curNode);
+ curSeq= curNode.getAttribute("seq");
+ curDir=!((BidirectedEdge) e).getDir(curNode);
+ curSeq = curDir?curSeq:Alphabet.DNA.complement(curSeq);
+
+ seq.append(curSeq.subSequence(0, curSeq.length()-(curNode==peekNode()?
+ 0:BidirectedGraph.getKmerSize())));
+ }
+ }
+ return seq.toSequence();
+ }
+ /*
+ * Add a path to the current path. The path to be added must start with the last node
+ * of the current path.
+ */
+ public void join(BidirectedPath bridge) {
+ if(bridge==null || bridge.size() <=1)
+ return;
+ if(bridge.getRoot() != peekNode()){
+ LOG.error("Cannot join path with disagreed first node " + bridge.getRoot().getId());
+ return;
+ }
+ if(((BidirectedEdge) bridge.getEdgePath().get(0)).getDir((AbstractNode) bridge.getRoot())
+ == ((BidirectedEdge) peekEdge()).getDir((AbstractNode) peekNode())){
+ LOG.error("Conflict direction from the first node " + bridge.getRoot().getId());
+ return;
+ }
+ //TODO: need a way to check coverage consistent
+
+
+ for(Edge e:bridge.getEdgePath()){
+ add(e);
+ }
+
+ coverage=Math.min(coverage, bridge.coverage);
+ }
+
+ public int getDeviation(){
+ return this.deviation;
+ }
+ public void setDeviation(int deviation){
+ this.deviation=deviation;
+ }
+
+ public double getCoverage(){
+ return coverage;
+ }
+ /**
+ *
+ * @return average depth of this path
+ */
+// public double averageCov(){
+// int len=0;
+// double res=0;
+// for(Node n:getNodePath()){
+// if(BidirectedGraph.isUnique(n)){
+// Sequence seq = (Sequence) n.getAttribute("seq");
+// double cov = Double.parseDouble(seq.getName().split("_")[5]);
+// len+=(n==getRoot())?seq.length():seq.length()-BidirectedGraph.getKmerSize();
+// res+=seq.length()*cov;
+// }
+// }
+// return res/len;
+// }
+
+ public int length() {
+ int retval = 0;
+ for(Node n:getNodePath()){
+ Sequence seq = (Sequence) n.getAttribute("seq");
+ retval+=(n==getRoot())?seq.length():seq.length()-BidirectedGraph.getKmerSize();
+ }
+ return retval;
+ }
+}
diff --git a/src/dev/java/japsadev/bio/hts/newscarf/GraphExplore.java b/src/dev/java/japsadev/bio/hts/newscarf/GraphExplore.java
new file mode 100644
index 0000000..7f85338
--- /dev/null
+++ b/src/dev/java/japsadev/bio/hts/newscarf/GraphExplore.java
@@ -0,0 +1,127 @@
+package japsadev.bio.hts.newscarf;
+
+import java.io.IOException;
+import java.util.Iterator;
+import org.graphstream.graph.*;
+import japsa.seq.Sequence;
+
+
+public class GraphExplore {
+
+ public static String spadesFolder="/home/sonhoanghguyen/Projects/scaffolding/data/spades_3.7/"; //imb desktop
+// public static String spadesFolder="/home/hoangnguyen/workspace/data/spades/"; //sony
+// public static String spadesFolder="/home/s_hoangnguyen/Projects/scaffolding/test-graph/spades/"; //dell
+
+
+ public static void main(String args[]) {
+ try {
+ new GraphExplore();
+ } catch (IOException e) {
+ // TODO Auto-generated catch block
+ e.printStackTrace();
+ }
+
+
+ }
+
+ public GraphExplore() throws IOException{
+ //System.setProperty("org.graphstream.ui.renderer", "org.graphstream.ui.j2dviewer.J2DGraphRenderer");
+// String sample="EcK12S-careful";
+// String sample="Kp2146-careful";
+// String sample="meta-careful";
+ String sample="cp_S5";
+
+ HybridAssembler ass = new HybridAssembler(spadesFolder+sample+"/assembly_graph.fastg");
+ BidirectedGraph graph= ass.simGraph;
+// graph.addAttribute("ui.quality");
+// graph.addAttribute("ui.antialias");
+ graph.addAttribute("ui.stylesheet", styleSheet);
+ graph.addAttribute("ui.default.title", "New real-time hybrid assembler");
+
+
+ graph.display();
+
+ System.out.println("Node: " + graph.getNodeCount() + " Edge: " + graph.getEdgeCount());
+
+
+ for (Node node : graph) {
+// node.addAttribute("ui.label", node.getId());
+// node.setAttribute("ui.style", "text-offset: -10;");
+ if(BidirectedGraph.isUnique(node))
+ node.setAttribute("ui.class", "marked");
+ }
+
+ //explore(graph.getNode("A"));
+
+ /*
+ * Testing reduce function
+ */
+ try {
+ HybridAssembler.promptEnterKey();
+ ass.reduceFromSPAdesPaths(spadesFolder+sample+"/contigs.paths");
+ HybridAssembler.promptEnterKey();
+ ass.assembly(spadesFolder+sample+"/assembly_graph.sam");
+
+ } catch (IOException e) {
+ // TODO Auto-generated catch block
+ e.printStackTrace();
+ }
+
+ //TODO: thorough cleaning... should have flag dead for each node
+// boolean dead=true;
+// while(dead){
+// dead=false;
+// for (Node node : graph) {
+// if((node.getDegree() < 2)
+//// || (node.getDegree()==0 && ((Sequence)node.getAttribute("seq")).length() < 1000)
+// ){
+// graph.removeNode(node);
+// dead=true;
+// }
+//
+// }
+// }
+ System.out.println("Node: " + graph.getNodeCount() + " Edge: " + graph.getEdgeCount());
+
+ /*
+ * Testing BidirectedEdge id pattern
+ */
+// String pattern = "^\\[([0-9\\+\\-]*)\\]([oi])\\[([0-9\\+\\-]*)\\]([oi])$";
+// // Create a Pattern object
+// Pattern r = Pattern.compile(pattern);
+// // Now create matcher object.
+// String id="[3]i[4+8+]o";
+// Matcher m = r.matcher(id);
+//
+// if(m.find()){
+// System.out.println(m.group(1)+"|"+m.group(2)+"|"+m.group(3)+"|"+m.group(4));
+// } else
+// System.out.println("Fuck");
+ }
+
+ public void explore(Node source) {
+ Iterator k = source.getBreadthFirstIterator();
+
+ while (k.hasNext()) {
+ Node next = k.next();
+ next.setAttribute("ui.class", "marked");
+ sleep();
+ }
+ }
+
+ protected void sleep() {
+ try { Thread.sleep(1000); } catch (Exception e) {}
+ }
+
+ protected String styleSheet =
+ "node {" +
+ " fill-color: black;" +
+ "}" +
+ "node.marked {" +
+ " fill-color: red;" +
+ "}" +
+ "edge.marked {" +
+ " fill-color: red;" +
+ "}";
+
+}
\ No newline at end of file
diff --git a/src/dev/java/japsadev/bio/hts/newscarf/GraphTest.java b/src/dev/java/japsadev/bio/hts/newscarf/GraphTest.java
new file mode 100644
index 0000000..9d231c1
--- /dev/null
+++ b/src/dev/java/japsadev/bio/hts/newscarf/GraphTest.java
@@ -0,0 +1,54 @@
+package japsadev.bio.hts.newscarf;
+
+import java.io.File;
+import java.io.PrintWriter;
+import java.util.Collections;
+import java.util.regex.Matcher;
+import java.util.regex.Pattern;
+
+import org.graphstream.graph.Graph;
+import org.graphstream.graph.implementations.AdjacencyListGraph;
+
+import htsjdk.samtools.*;
+
+public class GraphTest {
+ public static void main(String[] args) throws Exception {
+// System.setProperty("org.graphstream.ui.renderer",
+// "org.graphstream.ui.j2dviewer.J2DGraphRenderer");
+// Graph g = new AdjacencyListGraph("g");
+// g.addNode("A").addAttribute("xyz", new double[] { 0, 0 });
+// g.addNode("B").addAttribute("xyz", new double[] { 10, 10 });
+//
+// g.addEdge("AB", "A", "B", false)
+// .addAttribute(
+// "ui.points",
+// (Object) new double[] { 0, 0, 0, 0, 5, 0, 5, 10, 0, 10,
+// 10, 0 });
+//
+// g.addAttribute("ui.stylesheet", "edge {shape: polyline; }"); // or shape: cubic-curve
+//
+// g.display(false);
+
+// String data="/home/hoangnguyen/workspace/data/spades/EcK12S-careful/assembly_graph.fastg"; //sony
+//
+// final SamFileValidator validator=new SamFileValidator(new PrintWriter(System.out),8000);
+// validator.setIgnoreWarnings(true);
+// validator.setVerbose(true,1000);
+// validator.setErrorsToIgnore(Collections.singletonList(SAMValidationError.Type.MISSING_READ_GROUP));
+// SamReaderFactory factory=SamReaderFactory.makeDefault().validationStringency(ValidationStringency.SILENT);
+// SamReader samReader=factory.open(new File(data));
+// boolean ischeck=false;
+// SAMRecordIterator iter = samReader.iterator();
+// System.out.println(iter.next());
+//
+// samReader.close();
+
+ Pattern versionPattern = Pattern.compile("^Version:\\s(\\d+\\.\\d+\\.\\d+).*");
+ String line="Version: 0.7.5a-r405";
+ Matcher matcher =versionPattern.matcher(line);
+ if (matcher.find()){
+ System.out.println(line);
+ System.out.println(matcher.group(1));
+ }
+ }
+}
diff --git a/src/dev/java/japsadev/bio/hts/newscarf/HybridAssembler.java b/src/dev/java/japsadev/bio/hts/newscarf/HybridAssembler.java
new file mode 100644
index 0000000..f25af64
--- /dev/null
+++ b/src/dev/java/japsadev/bio/hts/newscarf/HybridAssembler.java
@@ -0,0 +1,308 @@
+package japsadev.bio.hts.newscarf;
+
+import java.io.BufferedReader;
+import java.io.File;
+import java.io.FileReader;
+import java.io.IOException;
+import java.util.ArrayList;
+import java.util.Scanner;
+
+import org.graphstream.graph.Edge;
+import org.graphstream.graph.Node;
+import org.slf4j.Logger;
+import org.slf4j.LoggerFactory;
+
+import htsjdk.samtools.SAMRecord;
+import htsjdk.samtools.SAMRecordIterator;
+import htsjdk.samtools.SamInputResource;
+import htsjdk.samtools.SamReader;
+import htsjdk.samtools.SamReaderFactory;
+import htsjdk.samtools.ValidationStringency;
+
+public class HybridAssembler {
+ private static final Logger LOG = LoggerFactory.getLogger(HybridAssembler.class);
+
+// final BidirectedGraph origGraph;
+ public BidirectedGraph simGraph; //original and simplified graph should be separated, no???
+
+ public HybridAssembler(){
+// origGraph=new BidirectedGraph("batch");
+ simGraph=new BidirectedGraph("real");
+ }
+
+
+ public HybridAssembler(String graphFile) throws IOException{
+ this();
+// origGraph.loadFromFile(graphFile);
+ simGraph.loadFromFile(graphFile);
+ }
+
+
+ public void assembly(String bamFile) throws IOException{
+ SamReaderFactory.setDefaultValidationStringency(ValidationStringency.SILENT);
+
+ SamReader reader;
+ if ("-".equals(bamFile))
+ reader = SamReaderFactory.makeDefault().open(SamInputResource.of(System.in));
+ else
+ reader = SamReaderFactory.makeDefault().open(new File(bamFile));
+
+ SAMRecordIterator iter = reader.iterator();
+
+ String readID = "";
+ //ReadFilling readFilling = null;
+ ArrayList samList = new ArrayList();;// alignment record of the same read;
+ BidirectedPath p = new BidirectedPath();
+ while (iter.hasNext()) {
+ SAMRecord rec = iter.next();
+ if (rec.getReadUnmappedFlag())
+ continue;
+// if (rec.getMappingQuality() < qual)
+// continue;
+
+ String refID = rec.getReferenceName().split("_")[1];
+ Alignment myRec = new Alignment(rec, simGraph.getNode(refID)); //FIXME: optimize
+
+ //////////////////////////////////////////////////////////////////
+ // make list of alignments of the same (Nanopore) read.
+
+ //not the first occurrance
+ if (!readID.equals("") && !readID.equals(myRec.readID)) {
+ //Collections.sort(samList);
+ //p=origGraph.pathFinding(samList);
+ p=simGraph.pathFinding(samList); // the graph MUST be the same as from new Alignment(...)
+
+ if(p!=null)
+ System.out.println("Final path found: " + p.getId());
+ reduce(p);
+// reduce2(p);
+ samList = new ArrayList();
+ //readID = myRec.readID;
+ }
+ readID = myRec.readID;
+ samList.add(myRec); // FIXME: (optimize) insert sort here
+
+ }// while
+ iter.close();
+
+ //outOS.close();
+ reader.close();
+
+ }
+ /*
+ * Read paths from contigs.path and reduce the graph
+ */
+ public void reduceFromSPAdesPaths(String paths) throws IOException{
+
+ BufferedReader pathReader = new BufferedReader(new FileReader(paths));
+
+ String s="", curpath="";
+ //Read contigs from contigs.paths and refer themselves to contigs.fasta
+ boolean flag=false;
+ while((s=pathReader.readLine()) != null){
+ if(s.contains("NODE")){
+ if(flag){
+ BidirectedPath path=new BidirectedPath(simGraph, curpath);
+ reduce(path);
+// reduce2(path);
+ }
+ flag=s.contains("'")?false:true;
+ curpath=new String();
+ continue;
+ }else if(flag){
+ curpath+=s;
+ }
+
+
+ }
+ pathReader.close();
+ }
+
+ /**
+ * Another reduce that doesn't remove the unique nodes
+ * Instead redundant edges are removed on a path way
+ * @param p Path to simplify the graph (from origGraph)
+ * @param target Subjected graph for the simplification
+ */
+ private void reduce(BidirectedPath p){
+ //do nothing if the path has only one node
+ if(p==null || p.getEdgeCount()<1)
+ return;
+
+ //loop over the edges of path (like spelling())
+ BidirectedNode markerNode = null,
+ curNodeFromSimGraph = (BidirectedNode) p.getRoot();
+
+ BidirectedPath curPath= null;
+ boolean markerDir=true, curDir;
+
+ if(BidirectedGraph.isUnique(curNodeFromSimGraph)){
+ markerNode=curNodeFromSimGraph;
+ markerDir=((BidirectedEdge) p.getEdgePath().get(0)).getDir(markerNode);
+ curPath = new BidirectedPath();
+ curPath.setRoot(curNodeFromSimGraph);
+ }
+
+
+ //search for an unique node as the marker.
+ ArrayList tobeRemoved = new ArrayList(),
+ tobeAdded = new ArrayList();
+ for(Edge e:p.getEdgePath()){
+
+ curNodeFromSimGraph=e.getOpposite(curNodeFromSimGraph);
+
+// curNodeFromSimGraph = simGraph.getNode(curNodeFromOrigGraph.getId()); //change back to Node belong to simGraph (instead of origGraph)
+ curDir=((BidirectedEdge) e).getDir(curNodeFromSimGraph);
+
+ if(BidirectedGraph.isUnique(curNodeFromSimGraph)){
+
+ if(markerNode!=null){
+ //this is when we have 1 jumping path (both ends are markers)
+ curPath.add(e);
+// LOG.info("Processing path {} with marker {}:{}:{} and curNode {}:{}:{}", curPath.getId(), markerNode.getId(), markerDir?"out":"in", markerNode.getGraph().getId(), curNodeFromSimGraph.getId(), curDir?"out":"in", curNodeFromSimGraph.getGraph().getId());
+ //create an edge connect markerNode to curNode with curPath
+ //Edge reducedEdge = simGraph.addEdge(markerNode, curNodeFromSimGraph, markerDir, curDir);
+ BidirectedEdge reducedEdge = new BidirectedEdge(markerNode, curNodeFromSimGraph, markerDir, curDir);
+
+// if(reducedEdge!=null)
+// reducedEdge.addAttribute("path", new BidirectedPath(curPath));
+
+ tobeAdded.add(reducedEdge);
+
+ //loop over curPath to find out edges needed to be removed
+ Node n0 = curPath.getRoot(),
+ n1 = null;
+ for(Edge ep:curPath.getEdgePath()){
+ n1 = ep.getOpposite(n0);
+ if(!BidirectedGraph.isUnique(n0) == BidirectedGraph.isUnique(n1)){
+ tobeRemoved.add((BidirectedEdge)ep);
+ }
+
+// if(!BidirectedGraph.isUnique(n1)){
+// n1.setAttribute("cov", n1.getNumber("cov")-markerNode.getNumber("cov"));
+// LOG.info("...coverage of " + n1.getAttribute("name") + " now is " + n1.getNumber("cov"));
+// }
+
+ n0=n1;
+ }
+
+ }
+
+
+ markerNode=curNodeFromSimGraph;
+ markerDir=!curDir; //in-out, out-in
+ curPath= new BidirectedPath();
+ curPath.setRoot(curNodeFromSimGraph);
+ }
+ else{
+ if(markerNode!=null){
+ curPath.add(e);
+ }
+ }
+
+ }
+
+ //remove appropriate edges
+ for(BidirectedEdge e:tobeRemoved){
+ LOG.info("REMOVING EDGE " + e.getId() + " from " + e.getNode0().getGraph().getId() + "-" + e.getNode1().getGraph().getId());
+ LOG.info("before: \n\t" + simGraph.printEdgesOfNode(e.getNode0()) + "\n\t" + simGraph.printEdgesOfNode(e.getNode1()));
+ simGraph.removeEdge(e.getId());
+ LOG.info("after: \n\t" + simGraph.printEdgesOfNode(e.getNode0()) + "\n\t" + simGraph.printEdgesOfNode(e.getNode1()));
+ }
+
+ //add appropriate edges
+ for(BidirectedEdge e:tobeAdded){
+ LOG.info("ADDING EDGE " + e.getId()+ " from " + e.getNode0().getGraph().getId() + "-" + e.getNode1().getGraph().getId());
+ LOG.info("before: \n\t" + simGraph.printEdgesOfNode(e.getNode0()) + "\n\t" + simGraph.printEdgesOfNode(e.getNode1()));
+
+ Edge reducedEdge = simGraph.addEdge(e.getSourceNode(),e.getTargetNode(),e.getDir0(),e.getDir1());
+ if(reducedEdge!=null){
+// reducedEdge.addAttribute("ui.label", reducedEdge.getId());
+// reducedEdge.setAttribute("ui.style", "text-offset: -10; text-alignment: along;");
+ reducedEdge.addAttribute("isReducedEdge", true);
+ reducedEdge.setAttribute("ui.class", "marked");
+// reducedEdge.addAttribute("layout.weight", 10);
+ }
+ LOG.info("after: \n\t" + simGraph.printEdgesOfNode(e.getNode0()) + "\n\t" + simGraph.printEdgesOfNode(e.getNode1()));
+
+ }
+
+ }
+ /**
+ * Another reduce that doesn't need to know unique contig
+ * @param p Path to simplify the graph (from origGraph)
+ * @param target Subjected graph for the simplification
+ */
+ private void reduce2(BidirectedPath p){
+ //do nothing if the path has only one node
+ if(p==null || p.getEdgeCount()<1)
+ return;
+ double coverage = p.getCoverage();
+ //loop over the edges of path (like spelling())
+ BidirectedNode firstNode = (BidirectedNode) p.getRoot(),
+ lastNode = (BidirectedNode) p.peekNode();
+
+ boolean firstDir=((BidirectedEdge)p.getEdgePath().get(0)).getDir(firstNode),
+ lastDir=((BidirectedEdge)p.peekEdge()).getDir(lastNode);
+
+ //search for an unique node as the marker.
+ ArrayList tobeRemoved = new ArrayList();
+ BidirectedNode curNode = firstNode;
+ boolean curDir;
+ for(Edge e:p.getEdgePath()){
+
+ double curCoverage=curNode.getNumber("cov"),
+ nextCoverage=e.getOpposite(curNode).getNumber("cov");
+ //if current node has the same coverage as path coverage
+ if(covLeft(curCoverage, coverage)==0 || covLeft(nextCoverage, coverage)==0){
+ tobeRemoved.add((BidirectedEdge) e);
+ }
+
+ curNode=e.getOpposite(curNode);
+ }
+
+ //remove appropriate edges
+ for(BidirectedEdge e:tobeRemoved){
+ LOG.info("REMOVING EDGE " + e.getId() + " from " + e.getNode0().getGraph().getId() + "-" + e.getNode1().getGraph().getId());
+ LOG.info("before: \n\t" + simGraph.printEdgesOfNode(e.getNode0()) + "\n\t" + simGraph.printEdgesOfNode(e.getNode1()));
+ simGraph.removeEdge(e.getId());
+ LOG.info("after: \n\t" + simGraph.printEdgesOfNode(e.getNode0()) + "\n\t" + simGraph.printEdgesOfNode(e.getNode1()));
+ }
+
+ //add appropriate edges
+ Edge reducedEdge = simGraph.addEdge(firstNode,lastNode,firstDir,lastDir);
+
+
+ }
+ private double covLeft(double cov, double pathCov){
+ double retval=0;
+ //TODO: need statistics here...
+ if((cov-pathCov)/pathCov > .2){
+ retval=cov-pathCov;
+ }
+ return retval;
+ }
+
+
+ @SuppressWarnings("resource")
+ public static void promptEnterKey(){
+ System.out.println("Press \"ENTER\" to continue...");
+ Scanner scanner = new Scanner(System.in);
+ scanner.nextLine();
+ }
+
+ protected void sleep() {
+ try { Thread.sleep(1000); } catch (Exception e) {}
+ }
+
+ public static void main(String[] argv) throws IOException{
+ HybridAssembler hbAss = new HybridAssembler(GraphExplore.spadesFolder+"EcK12S-careful/assembly_graph.fastg");
+ //For SAM file, run bwa first on the edited assembly_graph.fastg by running:
+ //awk -F '[:;]' -v q=\' 'BEGIN{flag=0;}/^>/{if(index($1,q)!=0) flag=0; else flag=1;}{if(flag==1) print $1;}' ../EcK12S-careful/assembly_graph.fastg > Eck12-careful.fasta
+ //TODO: need to make this easier
+
+ hbAss.assembly(GraphExplore.spadesFolder+"EcK12S-careful/assembly_graph.sam");
+
+ }
+
+}
diff --git a/src/dev/java/japsadev/bio/hts/newscarf/MetaRange.java b/src/dev/java/japsadev/bio/hts/newscarf/MetaRange.java
new file mode 100644
index 0000000..95259e6
--- /dev/null
+++ b/src/dev/java/japsadev/bio/hts/newscarf/MetaRange.java
@@ -0,0 +1,140 @@
+package japsadev.bio.hts.newscarf;
+
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.List;
+
+/**
+ * A meta range is composed of sub ranges. Its left endpoint is
+ * the leftmost point of its subranges and its right endpoint is
+ * the rightmost point of its subranges.
+ */
+public class MetaRange extends Range{
+
+ /**
+ * meta Ranges are composed of lists of subRanges.
+ */
+ private List subranges;
+
+ public MetaRange(Range initialSubRange){
+ super(initialSubRange.getLeft(),initialSubRange.getRight());
+ this.subranges = new ArrayList();
+ this.subranges.add(initialSubRange);
+ }
+
+ public void join(MetaRange other){
+ verifyHomo(other);
+
+ this.setLeft(Math.min(this.getLeft(),other.getLeft()));
+ this.setRight(Math.max(this.getRight(),other.getRight()));
+ this.subranges.addAll(other.getSubRanges());
+
+ other.invalidate();
+ }
+
+ public void setSubRanges(List subRanges){
+ this.subranges = subRanges;
+ }
+
+ public List getSubRanges(){
+ return this.subranges;
+ }
+
+ private void invalidate(){
+ this.setSubRanges(null);
+ this.setLeft(0);
+ this.setRight(0);
+ }
+
+ public boolean isInvalid(){
+ return this.getSubRanges()==null;
+ }
+
+ /**
+ * This code is retrieved from StackOverFlow: credits to augray :).
+ * From here on, the term "overlap" actually means ranges that
+ * belong to the same homo-group
+ *
+ * Given a list of Ranges, returns a list of Range groups where
+ * the Ranges in the groups all overlap. So if A overlaps with B and
+ * B with C, then A,B,and C will be returned in a group. Supposing Ranges
+ * D and E share nothing with A, B, or C, but do share with each other, D and E
+ * will be returned as a separate group from A,B,C.
+ * @param baseRanges
+ * @return
+ */
+ public static List> getOverlappingGroups(List baseRanges){
+ List baseMetaRanges = toMetaRanges(baseRanges);
+
+ List mergedMetaRanges = getMergedMetaRanges(baseMetaRanges);
+
+ List> RangeGroups = metaRangesToGroups(mergedMetaRanges);
+ return RangeGroups;
+ }
+
+
+
+ private static List getMergedMetaRanges(
+ List metaRanges) {
+ if(metaRanges.isEmpty()){
+ return metaRanges;
+ }
+ //order the MetaRanges by their starting point.
+ Collections.sort(metaRanges);
+
+ //go through and merge the overlapping meta Ranges.
+ //This relies on the logic that if Range i overlaps with
+ //an Range that started before it, then once all the Ranges
+ //before i have been merged, Range i will have a starting point
+ //consecutive to the starting point of the the preceeding Range.
+ for(int i=0; i< metaRanges.size()-1; i++){
+ MetaRange thisRange = metaRanges.get(i);
+ MetaRange nextRange = metaRanges.get(i+1);
+
+ if(thisRange.isHomo(nextRange)){
+ nextRange.join(thisRange);
+ }
+ }
+
+ List resultRanges = new ArrayList();
+
+ //All Ranges from the original list either:
+ //(a) didn't overlap with any others
+ //(b) overlapped with others and were chosen to represent the merged group or
+ //(c) overlapped with others, were represented in the group in another
+ // MetaRange, and then marked as invalid.
+ //Go through and only add the valid Ranges to be returned.
+
+ for(MetaRange i : metaRanges){
+ if(!i.isInvalid()){
+ resultRanges.add(i);
+ }
+ }
+ return resultRanges;
+ }
+
+ /**
+ * Convert a list of MetaRanges into groups of Ranges.
+ * @param mergedMetaRanges
+ * @return
+ */
+ private static List> metaRangesToGroups(
+ List mergedMetaRanges) {
+ List> groups = new ArrayList<>(mergedMetaRanges.size());
+ for(MetaRange metaRange : mergedMetaRanges){
+ groups.add(metaRange.getSubRanges());
+ }
+ return groups;
+ }
+
+ private static List toMetaRanges(
+ List baseRanges) {
+ ArrayList