New Upstream Snapshot - hisat2
Ready changes
Summary
Merged new upstream version: 2.2.1+git20220316.1.5086938 (was: 2.2.1).
Resulting package
Built on 2023-01-25T02:19 (took 11m12s)
The resulting binary packages can be installed (if you have the apt repository enabled) by running one of:
apt install -t fresh-snapshots hisat2-dbgsymapt install -t fresh-snapshots hisat2apt install -t fresh-snapshots python3-hisat2
Lintian Result
- hisat2-dbgsym_2.2.1+git20220316.1.5086938-1~jan+nus4_amd64.deb
- hisat2_2.2.1+git20220316.1.5086938-1~jan+nus4.dsc
- hisat2_2.2.1+git20220316.1.5086938-1~jan+nus4_amd64.buildinfo
- hisat2_2.2.1+git20220316.1.5086938-1~jan+nus4_amd64.deb
- python3-hisat2_2.2.1+git20220316.1.5086938-1~jan+nus4_all.deb
- hisat2_2.2.1+git20220316.1.5086938-1~jan+nus4_amd64.changes
Diff
diff --git a/.gitattributes b/.gitattributes
deleted file mode 100644
index 8edfb9b..0000000
--- a/.gitattributes
+++ /dev/null
@@ -1 +0,0 @@
-*.pbxproj binary merge=union
diff --git a/.gitignore b/.gitignore
deleted file mode 100644
index 74372a8..0000000
--- a/.gitignore
+++ /dev/null
@@ -1,43 +0,0 @@
-*~
-*.dSYM
-.DS_Store
-tags
-*-debug
-*-s
-*-l
-hisat2.xcodeproj/project.xcworkspace
-hisat2.xcodeproj/xcuserdata
-hisat2.xcodeproj/xcshareddata
-*.patch
-
-build_automaton
-build_index
-clean_alignment
-determinize
-gcsa_alignment
-gcsa_test
-hisat2-repeat
-
-hisat2_test/*.bt2
-hisat2_test/*.ht2
-hisat2_test/*.sam
-hisat2_test/paper_example.malignment.automaton
-hisat2_test/paper_example.malignment.backbone
-hisat2_test/paper_example.malignment.gcsa
-hisat2_test/kim_example*.malignment.automaton
-hisat2_test/kim_example*.malignment.backbone
-hisat2_test/kim_example*.malignment.gcsa
-hisat2_test/genome*
-hisat2_test/2*
-hisat2_test/snp142*
-hisat2_test/testset*
-
-.idea
-.vscode
-
-.ht2lib-obj*
-*.a
-*.so
-docs/_site
-docs/*.lock
-docs/.*-cache
diff --git a/HISAT2.png b/HISAT2.png
new file mode 100644
index 0000000..6bf6b3f
Binary files /dev/null and b/HISAT2.png differ
diff --git a/README.md b/README.md
index 44ee862..b05de43 100644
--- a/README.md
+++ b/README.md
@@ -1,202 +1,102 @@
-# Graph-Based Genome Alignment and Genotyping with HISAT2 and HISAT-genotype
+# HISAT2
## Contact
[Daehwan Kim](https://kim-lab.org) (infphilo@gmail.com) and [Chanhee Park](https://www.linkedin.com/in/chanhee-park-97677297/) (parkchanhee@gmail.com)
-## Abstract
-
-Rapid advances in next-generation sequencing technologies have dramatically changed our ability to perform genome-scale analyses. The human reference genome used for most genomic analyses represents only a small number of individuals, limiting its usefulness for genotyping. We designed a novel method, HISAT2, for representing and searching an expanded model of the human reference genome, in which a large catalogue of known genomic variants and haplotypes is incorporated into the data structure used for searching and alignment. This strategy for representing a population of genomes, along with a fast and memory-efficient search algorithm, enables more detailed and accurate variant analyses than previous methods. We demonstrate two initial applications of HISAT2: HLA typing, a critical need in human organ transplantation, and DNA fingerprinting, widely used in forensics. These applications are part of HISAT-genotype, with performance not only surpassing earlier computational methods, but matching or exceeding the accuracy of laboratory-based assays.
-
-
-
-For more information, see the following websites:
-* [HISAT2 website](http://ccb.jhu.edu/software/hisat2)
-* [HISAT-genotype website](http://ccb.jhu.edu/software/hisat-genotype)
-
-## HISAT2
-HISAT2 is a fast and sensitive alignment program for mapping next-generation sequencing reads (whole-genome, transcriptome, and exome sequencing data) to a population of human genomes (as well as to a single reference genome). Based on an extension of BWT for a graph [1], we designed and implemented a graph FM index (GFM), an original approach and its first implementation to the best of our knowledge. In addition to using one global GFM index that represents general population, HISAT2 uses a large set of small GFM indexes that collectively cover the whole genome (each index representing a genomic region of 56 Kbp, with 55,000 indexes needed to cover human population). These small indexes (called local indexes) combined with several alignment strategies enable effective alignment of sequencing reads. This new indexing scheme is called Hierarchical Graph FM index (HGFM). We have developed HISAT2 based on the HISAT [2] and Bowtie 2 [3] implementations. See the [HISAT2 website](http://ccb.jhu.edu/software/hisat2/index.shtml) for
+## Overview
+HISAT2 is a fast and sensitive alignment program for mapping next-generation sequencing reads (whole-genome,
+transcriptome, and exome sequencing data) to a population of human genomes (as well as to a single reference genome).
+Based on an extension of BWT for a graph [1], we designed and implemented a graph FM index (GFM), an original approach
+and its first implementation to the best of our knowledge. In addition to using one global GFM index that represents
+general population, HISAT2 uses a large set of small GFM indexes that collectively cover the whole genome (each index
+representing a genomic region of 56 Kbp, with 55,000 indexes needed to cover human population). These small indexes
+(called local indexes) combined with several alignment strategies enable effective alignment of sequencing reads.
+This new indexing scheme is called Hierarchical Graph FM index (HGFM). We have developed HISAT2 based on the HISAT
+[2] and Bowtie 2 [3] implementations. See the [HISAT2 website](https://daehwankimlab.github.io/hisat2/) for
more information.
-A few notes:
-1) HISAT2's index (HGFM) size for the human reference genome and 12.3 million common SNPs is 6.2GB. The SNPs consist of 11 million single nucleotide polymorphisms, 728,000 deletions, and 555,000 insertions. Insertions and deletions used in this index are small (usually <20bp). We plan to incorporate structural variations (SV) into this index.
+
-2) The memory footprint of HISAT2 is relatively low, 6.7GB.
+Getting started
+============
+HISAT2 requires a 64-bit computer running either Linux or Mac OS X and at least 8 GB of RAM.
+A few notes:
+1) HISAT2's index (HGFM) size for the human reference genome and 12.3 million common SNPs is 6.2 GB. The SNPs consist of 11 million single nucleotide polymorphisms, 728,000 deletions, and 555,000 insertions. Insertions and deletions used in this index are small (usually <20bp). We plan to incorporate structural variations (SV) into this index.
+2) The memory footprint of HISAT2 is relatively low, 6.7 GB.
3) The runtime of HISAT2 is estimated to be slightly slower than HISAT (30–100% slower for some data sets).
-
4) HISAT2 provides greater accuracy for alignment of reads containing SNPs.
+5) Use [HISAT-3N](https://daehwankimlab.github.io/hisat2/hisat-3n/) to align nucleotide converted sequencing reads
+ include [BS-seq], [SLAM-seq], [scBS-seq], [scSLAM-seq], [TAB-seq], [oxBS-seq], [TAPS] and [EM-seq].
+ This alignment process requires about 10 GB of RAM.
+6) HISAT2 repository is seperated with HISAT-genotype repository.
+ Please see the link below for [HISAT-genotype repository](https://github.com/DaehwanKimLab/hisat-genotype)
+ and [HISAT-genotype homepage](https://daehwankimlab.github.io/hisat-genotype/).
-5) We released a first (beta) version of HISAT2 in September 8, 2015.
-
-## License
-
-[GPL-3.0](LICENSE)
-
-# For reviwers, follow the instructions below to reproduce some of the results in the manuscript.
-
-## Code
-
-A specific version of [HISAT2 and HISAT-genotype](http://github.com/infphilo/hisat2) at GitHub is used (a branch name: hisat2_v2.2.0_beta).
-
-## Initial setup
-
-HISAT-genotype requires a 64-bit computer running either Linux or Mac OS X and at least 8 GB of RAM (16 GB of RAM is preferred). All the commands used should be run from the Unix shell prompt within a terminal window and are prefixed with a '$' character.
-
-We refer to <b>hisat-genotype-top</b> as our top directory where all of our programs are located. <b>hisat-genotype-top</b> is a place holder that can be changed to another name according to user preference.
-Run the following commands to install HISAT2 and HISAT-genotype.
-
- $ git clone https://github.com/infphilo/hisat2 hisat-genotype-top
- $ cd hisat-genotype-top
- hisat-genotype-top$ git checkout hisat2_v2.2.0_beta
- hisat-genotype-top$ make hisat2-align-s hisat2-build-s hisat2-inspect-s
-
-To make the binaries built above and other python scripts available everywhere, add the hisat-genotype-top directory to the PATH environment variable (e.g. ~/.bashrc)
-
- export PATH=hisat-genotype-top:hisat-genotype-top/hisatgenotype_scripts:$PATH
- export PYTHONPATH=hisat-genotype-top/hisatgenotype_modules:$PYTHONPATH
-
-To reflect the change, run the following command:
-
- $ source ~/.bashrc
-
-Download real reads, simulated reads, and HISAT2 indexes, then move them into appropriate directories:
-
- hisat-genotype-top$ cd evaluation
- hisat-genotype-top/evaluation$ wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hisat2_20181025.tar.gz
- hisat-genotype-top/evaluation$ tar xvzf hisat2_20181025.tar.gz
- hisat-genotype-top/evaluation$ mkdir aligners aligners/bin; cd aligners/bin; ln -s ../../../hisat2* .; cd ../..
- hisat-genotype-top/evaluation$ mv hisat2/* .
- hisat-genotype-top/evaluation$ cd simulation; ./init.py; cd ../real; ./init.py; cd ..
-
-## Run HISAT2 on the following simulated and real data sets.
-### 10 million simulated read pairs with SNPs and with sequencing errors
-
- hisat-genotype-top/evaluation$ cd simulation/10M_DNA_mismatch_snp_reads_genome
- hisat-genotype-top/evaluation/simulation/10M_DNA_mismatch_snp_reads_genome$ ./calculate_read_cost.py --aligner-list hisat2 --paired-end --fresh
-
-### 10 million simulated read pairs with SNPs and without sequencing errors
-
- hisat-genotype-top/evaluation$ cd simulation/10M_DNA_snp_reads_genome
- hisat-genotype-top/evaluation/simulation/10M_DNA_snp_reads_genome$ ./calculate_read_cost.py --aligner-list hisat2 --paired-end --fresh
-
-### 10 million simulated read pairs without SNPs and with sequencing errors
- hisat-genotype-top/evaluation$ cd simulation/10M_DNA_mismatch_reads_genome
- hisat-genotype-top/evaluation/simulation/10M_DNA_mismatch_reads_genome$ ./calculate_read_cost.py --aligner-list hisat2 --paired-end --fresh
-
-### 10 million simulated read pairs without SNPs and without sequencing errors
- hisat-genotype-top/evaluation$ cd simulation/10M_DNA_reads_genome
- hisat-genotype-top/evaluation/simulation/10M_DNA_reads_genome$ ./calculate_read_cost.py --aligner-list hisat2 --paired-end --fresh
-
-### 10 million real read pairs
- hisat-genotype-top/evaluation$ cd real/DNA/10M
- hisat-genotype-top/evaluation/real/DNA/10M$ ./calculate_read_cost.py --aligner-list hisat2 --paired-end --fresh
-
-### Interpreting output
- Example alignment output for simulated reads
- aligned: 1000000, multi aligned: 2654390
- correctly mapped: 999963 (100.00%)
- uniquely and correctly mapped: 967631 (96.76%)
- 54694 reads per sec (all)
- Memory Usage: 86MB
-
-The above lines show that 1,000,000 read pairs are aligned and the total number of alignments is 2,654,390. 999,963 pairs (100.00%) are correctly aligned (e.g. one of the alignments is correct). 967,631 (96.76%) pairs are uniquely and correctly aligned. HISAT2 aligns 54,594 reads with a peak memory usage of 86 MB of RAM.
+## Install
+ git clone https://github.com/DaehwanKimLab/hisat2.git
+ cd hisat2
+ make
-Each run is expected to take up to several hours mostly due to the comparison of HISAT2’s reported alignments and true alignments and the expansion of repeat alignments.
+Usage
+============
+## Building an index
+`hisat2-build` builds a HISAT2 index from a set of DNA sequences. `hisat2-build` outputs a set of 6 files with
+suffixes `.1.ht2`, `.2.ht2`, `.3.ht2`, `.4.ht2`, `.5.ht2`, `.6.ht2`, .`7.ht2`, and `.8.ht2`.
+In the case of a large index these suffixes will have a `ht2l` termination.
+These files together constitute the index: they are all that is needed to align reads to that reference.
+The original sequence FASTA files are no longer used by HISAT2 once the index is built.
-## Details on HISAT-genotype run for HLA typing and assembly
+Example for HISAT2 index building:
-To create a directory where we perform our analysis for HLA typing and assembly, which here is referred to as hla-analysis but can be changed by the user, execute the following command.
-
- hisat-genotype-top/evaluation$ mkdir hla-analysis
+ hisat2-build genome.fa genome
-The current directory can be changed to hla-analysis as follows:
-
- hisat-genotype-top/evaluation$ cd hla-analysis
+## Alignment with HISAT2
-Additional program requirements: SAMtools (version 1.3 or later)
+Examples alignment with HISAT2:
-### Downloading a Graph Reference and Index
-The graph reference we are going to build incorporates variants of numerous HLA alleles into the linear reference using a graph. The graph reference also includes some known variants of other regions of the genome (e.g. common small variants). To copy the graph reference, type:
+ # for single-end FASTA reads DNA alignment
+ hisat2 -f -x genome -U reads.fa -S output.sam --no-spliced-alignment
- hisat-genotype-top/evaluation/hla-analysis$ mv ../hisat2-genotype/* .
+ # for paired-end FASTQ reads alignment
+ hisat2 -x genome -1 reads_1.fq -2 read2_2.fq -S output.sam
-### Typing and Assembly
-Since whole genome sequencing (WGS) data includes reads that are from the whole genome, the first step is to extract the reads that belong to the HLA genes by aligning them to the graph reference with HISAT2. We provide these extracted reads in hisat-genotype-top/evaluation/hla-analysis/ILMN_20181025.
-HISAT-genotype performs both HLA typing and assembly as follows.
-You can perform HLA typing and assembly for HLA-A gene on sequencing reads from the genome NA12892 (Illumina's HiSeq 2000 platform).
+## For more information, see the following websites:
+* [HISAT2 website](https://daehwankimlab.github.io/hisat2)
+* [HISAT-3N website](https://daehwankimlab.github.io/hisat2/hisat-3n/)
- hisat-genotype-top/evaluation/hla-analysis$ hisatgenotype_locus.py --base hla --locus-list A --assembly -1 ILMN_20181025/NA12892.hla.extracted.1.fq.gz -2 ILMN_20181025/NA12892.hla.extracted.2.fq.gz
-
-### DNA Fingerprinting
-This function can be performed with the same commands used for “Typing and Assembly” and just replacing --base hla with --base codis.
-
-### Interpreting Output
- Typing Output
- Number of reads aligned: 1507
- 1 A*02:01:01:02L (count: 571)
- 2 A*02:01:31 (count: 557)
- 3 A*02:20:02 (count: 557)
- 4 A*02:29 (count: 557)
- 5 A*02:321N (count: 556)
- 6 A*02:372 (count: 556)
- 7 A*02:610:02 (count: 556)
- 8 A*02:249 (count: 555)
- 9 A*02:479 (count: 555)
- 10 A*02:11:01 (count: 554)
-
-The above lines show the top ten alleles that the most number of reads are mapped to or compatible with. For example, the allele first ranked, A\*02:01:01:02L, is compatible with 571 reads. This raw estimate based on the number of reads should not be used to determine the two true alleles because the alleles that resemble both but are not true alleles often tend to be compatible with more reads than either of the true alleles. Thus, we apply a statistical model to identify the two true alleles as described in the main text.
-
- Abundance of alleles
- 1 ranked A*02:01:01:01 (abundance: 54.32%)
- 2 ranked A*11:01:01:01 (abundance: 45.20%)
- 3 ranked A*24:33 (abundance: 0.48%)
-
-The above rankings show the top three alleles that are most abundant in the sample. Normally, the top two alleles in this estimate (e.g. A\*02:01:01:01 and A\*11:01:01:01) are considered as the two alleles that best match a given sequencing data.
-
-Additional tutorials and details are available at the HISAT-genotype website: https://ccb.jhu.edu/hisat-genotype
-
-
-## Data
-
-The Data directory (`/data`) contains all input files for reproducing some of our results such as from the
-evaluation of HISAT2 and other programs using both simulated and real reads, from typing and assembling
-HLA genes of Illumina Platinum Genomes using HISAT-genotype, and from building a HISAT2 graph index.
-
-* **Simulated read pairs**
-
-| Type | Number of pairs | Path |
-| - | - | - |
-| SNPs and sequencing errors included | 10,000,000 | hisat-genotype-top/evaluation/reads/simulation/10M_DNA_mismatch_snp_reads_genome |
-| SNPs included | 10,000,000 | hisat-genotype-top/evaluation/reads/simulation/10M_DNA_snp_reads_genome |
-| Sequencing errors included | 10,000,000 | hisat-genotype-top/evaluation/reads/simulation/10M_DNA_mismatch_reads_genome |
-| No SNPs nor sequencing errors included | 10,000,000 | hisat-genotype-top/evaluation/reads/simulation/10M_DNA_reads_genome |
+## License
-Each directory comes with a true alignment file in SAM format so that users know where the reads were generated in the human reference genome.
+[GPL-3.0](LICENSE)
-* **Real read pairs**
+## References
-| Number of read pairs | Path |
-| - | - |
-| 10,000,000 | hisat-genotype-top/evaluation/reads/real/DNA |
+[1] Sirén J, Välimäki N, Mäkinen V (2014) Indexing graphs for path queries with applications in genome research. IEEE/ACM Transactions on Computational Biology and Bioinformatics 11: 375–388. doi: 10.1109/tcbb.2013.2297101
-* **Human reference genome, SNPs, haplotypes, and HISAT2's indexes**
+[2] Kim D, Langmead B, and Salzberg SL HISAT: a fast spliced aligner with low memory requirements, Nature methods, 2015
-| Type | Path |
-| - | - |
-| GRCh38 reference | hisat-genotype-top/evaluation/data/genome.fa |
-| SNPs | hisat-genotype-top/evaluation/data/genome.snp |
-| Haplotypes | hisat-genotype-top/evaluation/data/genome.haplotype |
-| HISAT2's prebuilt graph index for comparison with other aligners | hisat-genotype-top/evaluation/indexes/HISAT2/genome.[1-8].ht2 |
-| HISAT2's prebuilt graph index for genotyping | hisat-genotype-top/evaluation/hla-analysis/genotype_genome.[1-8].ht2 |
+[3] Langmead B, Salzberg SL: Fast gapped-read alignment with Bowtie 2. Nat Methods 2012, 9:357-359
-## References
+[HISAT2]:https://github.com/DaehwanKimLab/hisat2
+[BS-seq]: https://en.wikipedia.org/wiki/Bisulfite_sequencing
+[SLAM-seq]: https://www.nature.com/articles/nmeth.4435
+[scBS-seq]: https://www.nature.com/articles/nmeth.3035
+[scSLAM-seq]: https://www.nature.com/articles/s41586-019-1369-y
+[TAPS]: https://www.nature.com/articles/s41587-019-0041-2
+[oxBS-seq]: https://science.sciencemag.org/content/336/6083/934
+[TAB-seq]: https://www.cell.com/fulltext/S0092-8674%2812%2900534-X
+[EM-seq]: https://genome.cshlp.org/cgi/content/long/gr.266551.120
-[1] Sirén J, Välimäki N, Mäkinen V (2014) Indexing graphs for path queries with applications in genome research. IEEE/ACM Transactions on Computational Biology and Bioinformatics 11: 375–388. doi: 10.1109/tcbb.2013.2297101
+Publication
+============
+* ### HISAT2 and HISAT-genotype
+ Kim, D., Paggi, J.M., Park, C. et al. [Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype](https://www.nature.com/articles/s41587-019-0201-4). Nat Biotechnol 37, 907–915 (2019)
-[2] Kim D, Langmead B, and Salzberg SL HISAT: a fast spliced aligner with low memory requirements, Nature methods, 2015
+* ### HISAT-3N
+ Zhang, Y., Park, C., Bennett, C., Thornton, M. and Kim, D [Rapid and accurate alignment of nucleotide conversion sequencing reads with HISAT-3N](https://doi.org/10.1101/gr.275193.120) Genome Research 31(7): 1290-1295 (2021)
+
-[3] Langmead B, Salzberg SL: Fast gapped-read alignment with Bowtie 2. Nat Methods 2012, 9:357-359
diff --git a/debian/changelog b/debian/changelog
index 585e483..11b8041 100644
--- a/debian/changelog
+++ b/debian/changelog
@@ -1,3 +1,9 @@
+hisat2 (2.2.1+git20220316.1.5086938-1) UNRELEASED; urgency=low
+
+ * New upstream snapshot.
+
+ -- Debian Janitor <janitor@jelmer.uk> Wed, 25 Jan 2023 02:11:06 -0000
+
hisat2 (2.2.1-4) unstable; urgency=medium
* Team upload
diff --git a/debian/patches/compilerflags b/debian/patches/compilerflags
index 7f3cc75..607b2e5 100644
--- a/debian/patches/compilerflags
+++ b/debian/patches/compilerflags
@@ -1,9 +1,11 @@
Description: remove restrictive compiler flags
Author: Sascha Steinbiss <satta@debian.org>
Bug-Debian: https://bugs.debian.org/cgi-bin/bugreport.cgi?bug=834862
---- hisat2.orig/Makefile
-+++ hisat2/Makefile
-@@ -165,14 +165,14 @@
+Index: hisat2.git/Makefile
+===================================================================
+--- hisat2.git.orig/Makefile
++++ hisat2.git/Makefile
+@@ -165,14 +165,14 @@ ifneq (,$(findstring AMD64,$(PROCESSOR_A
endif
BITS_FLAG =
diff --git a/debian/patches/fix-freebsd b/debian/patches/fix-freebsd
index 7f80bc0..20b884b 100644
--- a/debian/patches/fix-freebsd
+++ b/debian/patches/fix-freebsd
@@ -3,9 +3,11 @@ Description: Identify *freebsd as a UNIX-like OS
build on *FreeBSD.
Author: Steven Chamberlain <stevenc@debian.org>
Bug-Debian: https://bugs.debian.org/cgi-bin/bugreport.cgi?bug=834864
---- a/hisat2
-+++ b/hisat2
-@@ -45,7 +45,7 @@
+Index: hisat2.git/hisat2
+===================================================================
+--- hisat2.git.orig/hisat2
++++ hisat2.git/hisat2
+@@ -45,7 +45,7 @@ while (-f $prog && -l $prog){
($vol,$script_path,$prog)
= File::Spec->splitpath($prog);
diff --git a/debian/patches/fix-manual b/debian/patches/fix-manual
index b827240..de9e488 100644
--- a/debian/patches/fix-manual
+++ b/debian/patches/fix-manual
@@ -1,8 +1,10 @@
From: Michael R. Crusoe <michael.crusoe@gmail.com>
Subject: Update manual path
---- hisat2.orig/Makefile
-+++ hisat2/Makefile
-@@ -548,16 +548,16 @@
+Index: hisat2.git/Makefile
+===================================================================
+--- hisat2.git.orig/Makefile
++++ hisat2.git/Makefile
+@@ -548,16 +548,16 @@ hisat2-bin: $(BIN_PKG_LIST) $(HISAT2_BIN
rm -rf .bin.tmp
.PHONY: doc
diff --git a/debian/patches/hardening b/debian/patches/hardening
index 0f60c53..cb61049 100644
--- a/debian/patches/hardening
+++ b/debian/patches/hardening
@@ -1,9 +1,11 @@
From: Michael R. Crusoe <crusoe@debian.org>
Subject: Enable the passing of extra compilation flags, for hardening
Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/249
---- hisat2.orig/Makefile
-+++ hisat2/Makefile
-@@ -176,7 +176,7 @@
+Index: hisat2.git/Makefile
+===================================================================
+--- hisat2.git.orig/Makefile
++++ hisat2.git/Makefile
+@@ -176,7 +176,7 @@ SSE_FLAG=-msse2
DEBUG_FLAGS = -O0 -g3 $(BITS_FLAG) $(SSE_FLAG)
DEBUG_DEFS = -DCOMPILER_OPTIONS="\"$(DEBUG_FLAGS) $(EXTRA_FLAGS)\""
@@ -12,7 +14,7 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/249
RELEASE_DEFS = -DCOMPILER_OPTIONS="\"$(RELEASE_FLAGS) $(EXTRA_FLAGS)\""
NOASSERT_FLAGS = -DNDEBUG
FILE_FLAGS = -D_LARGEFILE_SOURCE -D_FILE_OFFSET_BITS=64 -D_GNU_SOURCE
-@@ -306,7 +306,7 @@
+@@ -306,7 +306,7 @@ hisat-bp-bin: hisat_bp.cpp $(SEARCH_CPPS
$(INC) \
-o $@ $< \
$(SHARED_CPPS) $(HISAT_CPPS_MAIN) \
@@ -21,7 +23,7 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/249
hisat-bp-bin-debug: hisat_bp.cpp $(SEARCH_CPPS) $(SHARED_CPPS) $(HEADERS) $(SEARCH_FRAGMENTS)
$(CXX) $(DEBUG_FLAGS) \
-@@ -315,7 +315,7 @@
+@@ -315,7 +315,7 @@ hisat-bp-bin-debug: hisat_bp.cpp $(SEARC
$(INC) \
-o $@ $< \
$(SHARED_CPPS) $(HISAT_CPPS_MAIN) \
@@ -30,7 +32,7 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/249
#
# hisat2-repeat targets
-@@ -327,7 +327,7 @@
+@@ -327,7 +327,7 @@ hisat2-repeat: hisat2_repeat.cpp $(REPEA
$(INC) \
-o $@ $< \
$(SHARED_CPPS) $(HISAT2_REPEAT_CPPS_MAIN) \
@@ -39,7 +41,7 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/249
hisat2-repeat-debug: hisat2_repeat.cpp $(REPEAT_CPPS) $(SHARED_CPPS) $(HEADERS)
$(CXX) $(DEBUG_FLAGS) $(DEBUG_DEFS) $(EXTRA_FLAGS) \
-@@ -335,7 +335,7 @@
+@@ -335,7 +335,7 @@ hisat2-repeat-debug: hisat2_repeat.cpp $
$(INC) \
-o $@ $< \
$(SHARED_CPPS) $(HISAT2_REPEAT_CPPS_MAIN) \
@@ -48,7 +50,7 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/249
#
-@@ -348,7 +348,7 @@
+@@ -348,7 +348,7 @@ hisat2-build-s: hisat2_build.cpp $(SHARE
$(INC) \
-o $@ $< \
$(SHARED_CPPS) $(HISAT2_BUILD_CPPS_MAIN) \
@@ -57,7 +59,7 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/249
hisat2-build-l: hisat2_build.cpp $(SHARED_CPPS) $(HEADERS)
$(CXX) $(RELEASE_FLAGS) $(RELEASE_DEFS) $(EXTRA_FLAGS) \
-@@ -356,7 +356,7 @@
+@@ -356,7 +356,7 @@ hisat2-build-l: hisat2_build.cpp $(SHARE
$(INC) \
-o $@ $< \
$(SHARED_CPPS) $(HISAT2_BUILD_CPPS_MAIN) \
@@ -66,7 +68,7 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/249
hisat2-build-s-debug: hisat2_build.cpp $(SHARED_CPPS) $(HEADERS)
$(CXX) $(DEBUG_FLAGS) $(DEBUG_DEFS) $(EXTRA_FLAGS) \
-@@ -364,7 +364,7 @@
+@@ -364,7 +364,7 @@ hisat2-build-s-debug: hisat2_build.cpp $
$(INC) \
-o $@ $< \
$(SHARED_CPPS) $(HISAT2_BUILD_CPPS_MAIN) \
@@ -75,7 +77,7 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/249
hisat2-build-l-debug: hisat2_build.cpp $(SHARED_CPPS) $(HEADERS)
$(CXX) $(DEBUG_FLAGS) $(DEBUG_DEFS) $(EXTRA_FLAGS) \
-@@ -372,7 +372,7 @@
+@@ -372,7 +372,7 @@ hisat2-build-l-debug: hisat2_build.cpp $
$(INC) \
-o $@ $< \
$(SHARED_CPPS) $(HISAT2_BUILD_CPPS_MAIN) \
@@ -84,7 +86,7 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/249
#
# hisat2 targets
-@@ -384,7 +384,7 @@
+@@ -384,7 +384,7 @@ hisat2-align-s: hisat2.cpp $(SEARCH_CPPS
$(INC) $(SEARCH_INC) \
-o $@ $< \
$(SHARED_CPPS) $(HISAT2_CPPS_MAIN) \
@@ -93,7 +95,7 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/249
hisat2-align-l: hisat2.cpp $(SEARCH_CPPS) $(SHARED_CPPS) $(HEADERS) $(SEARCH_FRAGMENTS)
$(CXX) $(RELEASE_FLAGS) $(RELEASE_DEFS) $(EXTRA_FLAGS) \
-@@ -392,7 +392,7 @@
+@@ -392,7 +392,7 @@ hisat2-align-l: hisat2.cpp $(SEARCH_CPPS
$(INC) $(SEARCH_INC) \
-o $@ $< \
$(SHARED_CPPS) $(HISAT2_CPPS_MAIN) \
@@ -102,7 +104,7 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/249
hisat2-align-s-debug: hisat2.cpp $(SEARCH_CPPS) $(SHARED_CPPS) $(HEADERS) $(SEARCH_FRAGMENTS)
$(CXX) $(DEBUG_FLAGS) \
-@@ -401,7 +401,7 @@
+@@ -401,7 +401,7 @@ hisat2-align-s-debug: hisat2.cpp $(SEARC
$(INC) $(SEARCH_INC) \
-o $@ $< \
$(SHARED_CPPS) $(HISAT2_CPPS_MAIN) \
@@ -111,7 +113,7 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/249
hisat2-align-l-debug: hisat2.cpp $(SEARCH_CPPS) $(SHARED_CPPS) $(HEADERS) $(SEARCH_FRAGMENTS)
$(CXX) $(DEBUG_FLAGS) \
-@@ -410,7 +410,7 @@
+@@ -410,7 +410,7 @@ hisat2-align-l-debug: hisat2.cpp $(SEARC
$(INC) $(SEARCH_INC) \
-o $@ $< \
$(SHARED_CPPS) $(HISAT2_CPPS_MAIN) \
@@ -120,7 +122,7 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/249
#
# hisat2-inspect targets
-@@ -423,7 +423,7 @@
+@@ -423,7 +423,7 @@ hisat2-inspect-s: hisat2_inspect.cpp $(H
$(INC) -I . \
-o $@ $< \
$(SHARED_CPPS) \
@@ -129,7 +131,7 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/249
hisat2-inspect-l: hisat2_inspect.cpp $(HEADERS) $(SHARED_CPPS)
$(CXX) $(RELEASE_FLAGS) \
-@@ -432,7 +432,7 @@
+@@ -432,7 +432,7 @@ hisat2-inspect-l: hisat2_inspect.cpp $(H
$(INC) -I . \
-o $@ $< \
$(SHARED_CPPS) \
@@ -138,7 +140,7 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/249
hisat2-inspect-s-debug: hisat2_inspect.cpp $(HEADERS) $(SHARED_CPPS)
$(CXX) $(DEBUG_FLAGS) \
-@@ -441,7 +441,7 @@
+@@ -441,7 +441,7 @@ hisat2-inspect-s-debug: hisat2_inspect.c
$(INC) -I . \
-o $@ $< \
$(SHARED_CPPS) \
@@ -147,7 +149,7 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/249
hisat2-inspect-l-debug: hisat2_inspect.cpp $(HEADERS) $(SHARED_CPPS)
$(CXX) $(DEBUG_FLAGS) \
-@@ -450,7 +450,7 @@
+@@ -450,7 +450,7 @@ hisat2-inspect-l-debug: hisat2_inspect.c
$(INC) -I . \
-o $@ $< \
$(SHARED_CPPS) \
diff --git a/debian/patches/python3 b/debian/patches/python3
index 34a9eea..a920bff 100644
--- a/debian/patches/python3
+++ b/debian/patches/python3
@@ -1,14 +1,16 @@
Author: Michael R. Crusoe <michael.crusoe@gmail.com>
Description: 2to3
---- hisat2.orig/evaluation/build_indexes.py
-+++ hisat2/evaluation/build_indexes.py
+Index: hisat2.git/evaluation/build_indexes.py
+===================================================================
+--- hisat2.git.orig/evaluation/build_indexes.py
++++ hisat2.git/evaluation/build_indexes.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
import sys, os
use_message = '''
-@@ -47,7 +47,7 @@
+@@ -47,7 +47,7 @@ def build_indexes():
assert False
else:
assert False
@@ -17,15 +19,17 @@ Description: 2to3
os.system(cmd)
os.chdir("..")
---- hisat2.orig/evaluation/generate_reads.py
-+++ hisat2/evaluation/generate_reads.py
+Index: hisat2.git/evaluation/generate_reads.py
+===================================================================
+--- hisat2.git.orig/evaluation/generate_reads.py
++++ hisat2.git/evaluation/generate_reads.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
import sys, os, random
from argparse import ArgumentParser, FileType
-@@ -18,7 +18,7 @@
+@@ -18,7 +18,7 @@ def shuffle_reads(read_fname, random_lis
assert len(random_list) == len(reads)
for i in random_list:
read = reads[random_list[i]]
@@ -34,7 +38,7 @@ Description: 2to3
read_file_out.close()
-@@ -73,11 +73,11 @@
+@@ -73,11 +73,11 @@ def simulate_reads():
data_dir_base = "../../../data"
def generate_reads(cmd):
@@ -48,15 +52,17 @@ Description: 2to3
shuffle_pairs("sim_1.fa", "sim_2.fa")
shuffle_reads_cmd = " mv sim_1.fa.shuffle sim_1.fa"
shuffle_reads_cmd += "; mv sim_2.fa.shuffle sim_2.fa"
---- hisat2.orig/evaluation/get_data.py
-+++ hisat2/evaluation/get_data.py
+Index: hisat2.git/evaluation/get_data.py
+===================================================================
+--- hisat2.git.orig/evaluation/get_data.py
++++ hisat2.git/evaluation/get_data.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
import sys, os
from argparse import ArgumentParser, FileType
-@@ -21,7 +21,7 @@
+@@ -21,7 +21,7 @@ def get_data(small = False):
if os.path.exists(file):
continue
wget_cmd = "wget %s/data/%s" % (data_root, file)
@@ -65,7 +71,7 @@ Description: 2to3
os.system(wget_cmd)
os.chdir("..")
-@@ -42,7 +42,7 @@
+@@ -42,7 +42,7 @@ def get_data(small = False):
continue
cmd = "wget %s/indexes/%s.tar.gz; tar xvzf %s.tar.gz; rm %s.tar.gz" % \
(data_root, aligner_dir, aligner_dir, aligner_dir)
@@ -74,7 +80,7 @@ Description: 2to3
os.system(cmd)
os.chdir("..")
-@@ -83,7 +83,7 @@
+@@ -83,7 +83,7 @@ def get_data(small = False):
continue
cmd = "wget %s/reads/%s/%s.tar.gz; tar xvzf %s.tar.gz; rm %s.tar.gz" % \
(data_root, type, file, file, file)
@@ -83,15 +89,17 @@ Description: 2to3
os.system(cmd)
os.chdir("..")
---- hisat2.orig/evaluation/get_programs.py
-+++ hisat2/evaluation/get_programs.py
+Index: hisat2.git/evaluation/get_programs.py
+===================================================================
+--- hisat2.git.orig/evaluation/get_programs.py
++++ hisat2.git/evaluation/get_programs.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
import sys, os
use_message = '''
-@@ -112,7 +112,7 @@
+@@ -112,7 +112,7 @@ def get_aligners():
cmd += "; source ./source_me.sh; make; cp bin/%s ../bin; cd .." % (program)
else:
assert False
@@ -100,15 +108,17 @@ Description: 2to3
os.system(cmd)
files = ["hisat2", "hisat2-align-s", "hisat2-build", "hisat2-build-s", "hisat2-inspect", "hisat2-inspect-s", "extract_splice_sites.py", "hisat2_extract_snps_haplotypes_UCSC.py", "hisat2_simulate_reads.py"]
---- hisat2.orig/evaluation/real/calculate_read_cost.py
-+++ hisat2/evaluation/real/calculate_read_cost.py
+Index: hisat2.git/evaluation/real/calculate_read_cost.py
+===================================================================
+--- hisat2.git.orig/evaluation/real/calculate_read_cost.py
++++ hisat2.git/evaluation/real/calculate_read_cost.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
import sys, os, subprocess, signal
import multiprocessing
-@@ -86,7 +86,7 @@
+@@ -86,7 +86,7 @@ def read_genome(genome_filename):
genome_file.close()
@@ -117,7 +127,7 @@ Description: 2to3
return chr_dic
-@@ -111,7 +111,7 @@
+@@ -111,7 +111,7 @@ def read_snp(snp_filename):
data = int(data)
snps[chr].append([snpID, type, int(pos), data])
@@ -126,7 +136,7 @@ Description: 2to3
return snps
-@@ -158,7 +158,7 @@
+@@ -158,7 +158,7 @@ def extract_splice_sites(gtf_fname):
gtf_file.close()
# Sort exons and merge where separating introns are <=5 bps
@@ -135,7 +145,7 @@ Description: 2to3
exons.sort()
tmp_exons = [exons[0]]
for i in range(1, len(exons)):
-@@ -170,7 +170,7 @@
+@@ -170,7 +170,7 @@ def extract_splice_sites(gtf_fname):
# Calculate and print the unique junctions
junctions = set()
@@ -144,7 +154,7 @@ Description: 2to3
for i in range(1, len(exons)):
junctions.add(to_junction_str([chrom, exons[i-1][1], exons[i][0]]))
-@@ -411,7 +411,7 @@
+@@ -411,7 +411,7 @@ def getSNPs(chr_snps, left, right):
high = mid - 1
snps = []
@@ -153,7 +163,7 @@ Description: 2to3
snp = chr_snps[i]
snpID, type, pos, data = snp
pos2 = pos
-@@ -539,7 +539,7 @@
+@@ -539,7 +539,7 @@ def extract_reads_and_pairs(chr_dic, sam
if cigar_op == "S":
if i != 0 and i != len(cigars) - 1:
@@ -162,7 +172,7 @@ Description: 2to3
if cigar_op in "MS":
ref_pos = right_pos
-@@ -596,8 +596,8 @@
+@@ -596,8 +596,8 @@ def extract_reads_and_pairs(chr_dic, sam
if hisat2:
XM, NM = HISAT2_XM, HISAT2_NM
if NM < MAX_EDIT:
@@ -173,7 +183,7 @@ Description: 2to3
found = False
me = "%s\t%s\t%d" % (read_id, chr, pos)
-@@ -611,8 +611,8 @@
+@@ -611,8 +611,8 @@ def extract_reads_and_pairs(chr_dic, sam
flag, chr, pos, cigar_str, XM, NM, mate_flag, mate_chr_str, mate_pos, mate_cigar_str, mate_XM, mate_NM = \
mate_flag, mate_chr, mate_pos, mate_cigar_str, mate_XM, mate_NM, flag, chr, pos, cigar_str, XM, NM
@@ -184,7 +194,7 @@ Description: 2to3
found = True
break
-@@ -764,7 +764,7 @@
+@@ -764,7 +764,7 @@ def cal_read_len(cigar_str):
def is_concordantly(read_id, flag, chr, pos, cigar_str, XM, NM, mate_flag, mate_chr, mate_pos, mate_cigar_str, mate_XM, mate_NM):
concord_length = 1000
@@ -193,7 +203,7 @@ Description: 2to3
pairs = {}
pairs[0] = [flag, chr, pos, cigar_str, XM, NM]
-@@ -819,7 +819,7 @@
+@@ -819,7 +819,7 @@ def pair_stat(pair_filename, gtf_junctio
# check concordantly
concord_align, segment_len = is_concordantly(read_id, flag, chr, pos, cigar_str, XM, NM, mate_flag, mate_chr, mate_pos, mate_cigar_str, mate_XM, mate_NM)
@@ -202,7 +212,7 @@ Description: 2to3
if junction_pair:
for junction_str, is_gtf_junction in pair_junctions:
-@@ -911,7 +911,7 @@
+@@ -911,7 +911,7 @@ def sql_execute(sql_db, sql_query):
def create_sql_db(sql_db):
if os.path.exists(sql_db):
@@ -211,7 +221,7 @@ Description: 2to3
return
columns = [
-@@ -961,14 +961,14 @@
+@@ -961,14 +961,14 @@ def write_analysis_data(sql_db, database
database_fname = database_name + "_" + paired + ".analysis"
database_file = open(database_fname, "w")
@@ -228,7 +238,7 @@ Description: 2to3
database_file.close()
-@@ -1072,7 +1072,7 @@
+@@ -1072,7 +1072,7 @@ def calculate_read_cost(single_end,
gtf_junctions.append(junction)
gtf_junctions = sorted(gtf_junctions, cmp=junction_cmp)
@@ -237,7 +247,7 @@ Description: 2to3
for paired in [False, True]:
if not paired and not single_end:
-@@ -1219,7 +1219,7 @@
+@@ -1219,7 +1219,7 @@ def calculate_read_cost(single_end,
if version != "":
version = int(version)
else:
@@ -246,7 +256,7 @@ Description: 2to3
if not RNA:
cmd += ["--no-spliced-alignment"]
-@@ -1441,7 +1441,7 @@
+@@ -1441,7 +1441,7 @@ def calculate_read_cost(single_end,
dummy_cmd = get_aligner_cmd(RNA, aligner, type, index_type, version, options, "../one.fq", "../two.fq", "/dev/null")
start_time = datetime.now()
if verbose:
@@ -255,7 +265,7 @@ Description: 2to3
if aligner in ["hisat2", "hisat", "bowtie", "bowtie2", "gsnap", "bwa"]:
proc = subprocess.Popen(dummy_cmd, stdout=open("/dev/null", "w"), stderr=subprocess.PIPE)
else:
-@@ -1451,7 +1451,7 @@
+@@ -1451,7 +1451,7 @@ def calculate_read_cost(single_end,
duration = finish_time - start_time
duration = duration.total_seconds()
if verbose:
@@ -264,7 +274,7 @@ Description: 2to3
loading_time = duration
# align all reads
-@@ -1459,7 +1459,7 @@
+@@ -1459,7 +1459,7 @@ def calculate_read_cost(single_end,
sweep_read_cmd = "cat ../%s ../%s > /dev/null" % (type_read1_fname, type_read2_fname)
else:
sweep_read_cmd = "cat ../%s > /dev/null" % (type_read1_fname)
@@ -273,7 +283,7 @@ Description: 2to3
os.system(sweep_read_cmd)
skip_alignment = False
-@@ -1470,7 +1470,7 @@
+@@ -1470,7 +1470,7 @@ def calculate_read_cost(single_end,
aligner_cmd = get_aligner_cmd(RNA, aligner, type, index_type, version, options, "../" + type_read1_fname, "../" + type_read2_fname, out_fname)
start_time = datetime.now()
if verbose:
@@ -282,7 +292,7 @@ Description: 2to3
if aligner in ["hisat2", "hisat", "bowtie", "bowtie2", "gsnap", "bwa", "vg", "minimap2"]:
proc = subprocess.Popen(aligner_cmd, stdout=open(out_fname, "w"), stderr=subprocess.PIPE)
else:
-@@ -1483,14 +1483,14 @@
+@@ -1483,14 +1483,14 @@ def calculate_read_cost(single_end,
if duration < 0.1:
duration = 0.1
if verbose:
@@ -300,7 +310,7 @@ Description: 2to3
# break
if aligner == "star" and type in ["", "gtf"]:
-@@ -1498,7 +1498,7 @@
+@@ -1498,7 +1498,7 @@ def calculate_read_cost(single_end,
elif aligner in ["hisat2", "hisat"] and type == "x2":
aligner_cmd = get_aligner_cmd(RNA, aligner, type, index_type, version, options, "../" + type_read1_fname, "../" + type_read2_fname, out_fname, 1)
if verbose:
@@ -309,7 +319,7 @@ Description: 2to3
start_time = datetime.now()
proc = subprocess.Popen(aligner_cmd, stdout=open(out_fname, "w"), stderr=subprocess.PIPE)
proc.communicate()
-@@ -1508,7 +1508,7 @@
+@@ -1508,7 +1508,7 @@ def calculate_read_cost(single_end,
if duration < 0.1:
duration = 0.1
if verbose:
@@ -318,7 +328,7 @@ Description: 2to3
elif aligner == "star" and type == "x2":
assert os.path.exists("SJ.out.tab")
os.system("awk 'BEGIN {OFS=\"\t\"; strChar[0]=\".\"; strChar[1]=\"+\"; strChar[2]=\"-\";} {if($5>0){print $1,$2,$3,strChar[$4]}}' SJ.out.tab > SJ.out.tab.Pass1.sjdb")
-@@ -1517,18 +1517,18 @@
+@@ -1517,18 +1517,18 @@ def calculate_read_cost(single_end,
continue
os.remove(file)
star_index_cmd = "STAR --genomeDir ./ --runMode genomeGenerate --genomeFastaFiles ../../../../data/genome.fa --sjdbFileChrStartEnd SJ.out.tab.Pass1.sjdb --sjdbOverhang 100 --runThreadN %d" % (num_threads)
@@ -341,7 +351,7 @@ Description: 2to3
proc = subprocess.Popen(aligner_cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
proc.communicate()
finish_time = datetime.now()
-@@ -1537,7 +1537,7 @@
+@@ -1537,7 +1537,7 @@ def calculate_read_cost(single_end,
if duration < 0.1:
duration = 0.1
if verbose:
@@ -350,7 +360,7 @@ Description: 2to3
os.system("mv Aligned.out.sam %s" % out_fname)
elif aligner == "tophat2":
os.system("samtools sort -n tophat_out/accepted_hits.bam accepted_hits; samtools view -h accepted_hits.bam > %s" % out_fname)
-@@ -1597,8 +1597,8 @@
+@@ -1597,8 +1597,8 @@ def calculate_read_cost(single_end,
sql_execute("../" + sql_db_name, sql_insert)
@@ -361,7 +371,7 @@ Description: 2to3
else:
sum = [0, 0, 0, 0, 0]
stat = read_stat(read_sam, gtf_junctions, chr_dic)
-@@ -1616,8 +1616,8 @@
+@@ -1616,8 +1616,8 @@ def calculate_read_cost(single_end,
(workdir, genome, "single", aligner_name, get_aligner_version(aligner), "no", i, mapped_reads, junction_reads, gtf_junction_reads, num_junctions, num_gtf_junctions, duration, platform.node(), " ".join(aligner_cmd))
sql_execute("../" + sql_db_name, sql_insert)
@@ -372,15 +382,17 @@ Description: 2to3
done_file.close()
---- hisat2.orig/evaluation/real/init.py
-+++ hisat2/evaluation/real/init.py
+Index: hisat2.git/evaluation/real/init.py
+===================================================================
+--- hisat2.git.orig/evaluation/real/init.py
++++ hisat2.git/evaluation/real/init.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
import sys, os, signal
import string, re
-@@ -51,7 +51,7 @@
+@@ -51,7 +51,7 @@ def init():
not os.path.exists(read_dir_base + read_dir + "/1.fq"):
continue
@@ -389,7 +401,7 @@ Description: 2to3
os.mkdir(read_dir)
os.chdir(read_dir)
-@@ -74,18 +74,18 @@
+@@ -74,18 +74,18 @@ def init():
if dir_name == "whole":
ln_cmd = "ln -s ../../%s%s/%s ." % (read_dir_base, read_dir, fq_1_name)
@@ -412,15 +424,17 @@ Description: 2to3
os.system(cmd)
os.system("ln -s ../../calculate_read_cost.py .")
---- hisat2.orig/evaluation/simulation/calculate_read_cost.py
-+++ hisat2/evaluation/simulation/calculate_read_cost.py
+Index: hisat2.git/evaluation/simulation/calculate_read_cost.py
+===================================================================
+--- hisat2.git.orig/evaluation/simulation/calculate_read_cost.py
++++ hisat2.git/evaluation/simulation/calculate_read_cost.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
import sys, os, subprocess
import multiprocessing
-@@ -170,18 +170,18 @@
+@@ -170,18 +170,18 @@ def read_repeatdb(repeat_filename):
repeat_db[rptRefName][rptName].allele[allele_idx].add_position(chr, pos, strand)
else:
@@ -444,7 +458,7 @@ Description: 2to3
repeat_right = 0
#for allele_id, repeatAllele in repeat.allele.items():
-@@ -242,11 +242,11 @@
+@@ -242,11 +242,11 @@ def repeat_to_genome_pos(repeat_db, repe
i = find_leftmost_pos(rmap, left)
if i >= len(rmap):
@@ -458,7 +472,7 @@ Description: 2to3
return
repeat = repeats[rmap[i][1]]
-@@ -260,7 +260,7 @@
+@@ -260,7 +260,7 @@ def repeat_to_genome_pos(repeat_db, repe
if (left >= rpos) and (right <= rpos + rlen):
offset = left - rpos
for genome_pos in allele.positions:
@@ -467,7 +481,7 @@ Description: 2to3
"""
"""
-@@ -284,7 +284,7 @@
+@@ -284,7 +284,7 @@ def read_genome(genome_filename):
genome_file.close()
@@ -476,7 +490,7 @@ Description: 2to3
return chr_dic
-@@ -331,7 +331,7 @@
+@@ -331,7 +331,7 @@ def extract_splice_sites(gtf_fname):
gtf_file.close()
# Sort exons and merge where separating introns are <=5 bps
@@ -485,7 +499,7 @@ Description: 2to3
exons.sort()
tmp_exons = [exons[0]]
for i in range(1, len(exons)):
-@@ -343,7 +343,7 @@
+@@ -343,7 +343,7 @@ def extract_splice_sites(gtf_fname):
# Calculate and print the unique junctions
junctions = set()
@@ -494,7 +508,7 @@ Description: 2to3
for i in range(1, len(exons)):
junctions.add(to_junction_str([chrom, exons[i-1][1], exons[i][0]]))
-@@ -381,7 +381,7 @@
+@@ -381,7 +381,7 @@ def read_repeat_info(repeat_filename):
repeat_dic[rep][allele].append([chr, pos, strand])
repeat_pos[rep][allele].add(pos)
@@ -503,7 +517,7 @@ Description: 2to3
def my_cmp(a, b):
if a[1] < b[1]:
return -1
-@@ -687,11 +687,11 @@
+@@ -687,11 +687,11 @@ def repeat_to_genome_alignment(repeat_db
i = find_leftmost_pos(rmap, left)
if i >= len(rmap):
@@ -517,7 +531,7 @@ Description: 2to3
return alignments
repeat = repeats[rmap[i][1]]
-@@ -775,7 +775,7 @@
+@@ -775,7 +775,7 @@ def extract_single(infilename,
prev_read_id = read_id
continue
@@ -526,7 +540,7 @@ Description: 2to3
for i in range(11, len(cols)):
col = cols[i]
# "nM" from STAR
-@@ -867,7 +867,7 @@
+@@ -867,7 +867,7 @@ def extract_single(infilename,
chr, pos, cigar_str = alignment
pos, cigar_str, NM_real = adjust_alignment(chr, pos, cigar_str)
p_str = "%s\t%s\t%d\t%s\tNM:i:%d" % (read_id, chr, pos, cigar_str, NM_real)
@@ -535,7 +549,7 @@ Description: 2to3
if aligner == "hisat2":
if prev_read_id != read_id:
-@@ -940,7 +940,7 @@
+@@ -940,7 +940,7 @@ def extract_pair(infilename,
num_aligned_reads, num_ualigned_reads = 0, 0
prev_read_id, pair_list = "", set()
@@ -544,7 +558,7 @@ Description: 2to3
prev_NH1, prev_NH2 = 0, 0
NH1_real, NH2_real = 0, 0
-@@ -969,7 +969,7 @@
+@@ -969,7 +969,7 @@ def extract_pair(infilename,
if read_id != prev_read_id:
num_pairs += 1
pair_list = set()
@@ -553,7 +567,7 @@ Description: 2to3
flag = int(flag)
canonical_pos1, canonical_pos2 = int(pos1), int(pos2)
-@@ -982,7 +982,7 @@
+@@ -982,7 +982,7 @@ def extract_pair(infilename,
continue
concordant = (flag & 0x2 != 0)
@@ -562,7 +576,7 @@ Description: 2to3
for i in range(11, len(cols)):
col = cols[i]
# "nM" from STAR
-@@ -1151,7 +1151,7 @@
+@@ -1151,7 +1151,7 @@ def extract_pair(infilename,
if p_str not in pair_reported:
pair_reported.add(p_str)
@@ -571,7 +585,7 @@ Description: 2to3
if not me in read_dic:
read_dic[me] = []
-@@ -1348,7 +1348,7 @@
+@@ -1348,7 +1348,7 @@ def compare_single_sam(RNA,
file.close()
temp_junctions, temp_gtf_junctions = set(), set()
@@ -580,7 +594,7 @@ Description: 2to3
if len(can_junctions) <= 0:
continue
-@@ -1492,7 +1492,7 @@
+@@ -1492,7 +1492,7 @@ def compare_single_sam(RNA,
break
if found:
@@ -589,7 +603,7 @@ Description: 2to3
mapped += 1
if snp_included:
snp_mapped += 1
-@@ -1501,7 +1501,7 @@
+@@ -1501,7 +1501,7 @@ def compare_single_sam(RNA,
if snp_included:
snp_unique_mapped += 1
if found_at_first:
@@ -598,7 +612,7 @@ Description: 2to3
first_mapped += 1
if snp_included:
snp_first_mapped += 1
-@@ -1527,7 +1527,7 @@
+@@ -1527,7 +1527,7 @@ def compare_single_sam(RNA,
false_can_junctions += 1
else:
false_noncan_junctions += 1
@@ -607,7 +621,7 @@ Description: 2to3
return mapped, unique_mapped, first_mapped, unmapped, aligned, multi_aligned, \
snp_mapped, snp_unique_mapped, snp_first_mapped, snp_unmapped, \
-@@ -1587,7 +1587,7 @@
+@@ -1587,7 +1587,7 @@ def compare_paired_sam(RNA,
file.close()
temp_junctions, temp_gtf_junctions = set(), set()
@@ -616,7 +630,7 @@ Description: 2to3
if len(can_junctions) <= 0:
continue
-@@ -1767,17 +1767,17 @@
+@@ -1767,17 +1767,17 @@ def compare_paired_sam(RNA,
break
if found:
@@ -637,7 +651,7 @@ Description: 2to3
first_mapped += 1
if snp_included:
snp_first_mapped += 1
-@@ -1802,7 +1802,7 @@
+@@ -1802,7 +1802,7 @@ def compare_paired_sam(RNA,
false_can_junctions += 1
else:
false_noncan_junctions += 1
@@ -646,7 +660,7 @@ Description: 2to3
return mapped, unique_mapped, first_mapped, unmapped, aligned, multi_aligned, \
-@@ -1832,9 +1832,9 @@
+@@ -1832,9 +1832,9 @@ def extract_mapped_unmapped(read_fname,
write = read_id in mapped_ids
if write:
@@ -658,7 +672,7 @@ Description: 2to3
read_file.close()
mapped_file.close()
-@@ -1865,7 +1865,7 @@
+@@ -1865,7 +1865,7 @@ def sql_execute(sql_db, sql_query):
"""
def create_sql_db(sql_db):
if os.path.exists(sql_db):
@@ -667,7 +681,7 @@ Description: 2to3
return
columns = [
-@@ -1936,14 +1936,14 @@
+@@ -1936,14 +1936,14 @@ def write_analysis_data(sql_db, genome_n
database_fname += "_single"
database_fname += ".analysis"
database_file = open(database_fname, "w")
@@ -684,7 +698,7 @@ Description: 2to3
database_file.close()
-@@ -2069,7 +2069,7 @@
+@@ -2069,7 +2069,7 @@ def calculate_read_cost(single_end,
type_sam_file.close()
if numreads <= 0:
continue
@@ -693,7 +707,7 @@ Description: 2to3
junctions, junctions_set = [], set()
type_junction_file = open(type_junction_fname)
-@@ -2208,7 +2208,7 @@
+@@ -2208,7 +2208,7 @@ def calculate_read_cost(single_end,
if version != "":
version = int(version)
else:
@@ -702,7 +716,7 @@ Description: 2to3
if not RNA:
cmd += ["--no-spliced-alignment"]
-@@ -2417,7 +2417,7 @@
+@@ -2417,7 +2417,7 @@ def calculate_read_cost(single_end,
if not RNA and readtype != "all":
continue
@@ -711,7 +725,7 @@ Description: 2to3
if options != "":
option_name = options.replace(' ', '').replace('-', '').replace(',', '')
aligner_name = aligner_name + '_' + option_name
-@@ -2457,7 +2457,7 @@
+@@ -2457,7 +2457,7 @@ def calculate_read_cost(single_end,
dummy_cmd = get_aligner_cmd(RNA, aligner, type, index_type, version, options, "../one.fa", "../two.fa", "/dev/null")
start_time = datetime.now()
if verbose:
@@ -720,7 +734,7 @@ Description: 2to3
if aligner in ["hisat2", "hisat", "bowtie", "bowtie2", "gsnap", "bwa"]:
proc = subprocess.Popen(dummy_cmd, stdout=open("/dev/null", "w"), stderr=subprocess.PIPE)
else:
-@@ -2467,14 +2467,14 @@
+@@ -2467,14 +2467,14 @@ def calculate_read_cost(single_end,
duration = finish_time - start_time
duration = duration.total_seconds()
if verbose:
@@ -737,7 +751,7 @@ Description: 2to3
if aligner in ["hisat2", "hisat", "bowtie", "bowtie2", "gsnap", "bwa", "vg", "minimap2"]:
proc = subprocess.Popen(aligner_cmd, stdout=open(out_fname, "w"), stderr=subprocess.PIPE)
else:
-@@ -2487,11 +2487,11 @@
+@@ -2487,11 +2487,11 @@ def calculate_read_cost(single_end,
if duration < 0.1:
duration = 0.1
if verbose:
@@ -751,7 +765,7 @@ Description: 2to3
if aligner == "star" and type in ["", "gtf"]:
os.system("mv Aligned.out.sam %s" % out_fname)
-@@ -2499,7 +2499,7 @@
+@@ -2499,7 +2499,7 @@ def calculate_read_cost(single_end,
aligner_cmd = get_aligner_cmd(RNA, aligner, type, index_type, version, options, "../" + type_read1_fname, "../" + type_read2_fname, out_fname, 1)
start_time = datetime.now()
if verbose:
@@ -760,7 +774,7 @@ Description: 2to3
proc = subprocess.Popen(aligner_cmd, stdout=open(out_fname, "w"), stderr=subprocess.PIPE)
proc.communicate()
finish_time = datetime.now()
-@@ -2508,7 +2508,7 @@
+@@ -2508,7 +2508,7 @@ def calculate_read_cost(single_end,
if duration < 0.1:
duration = 0.1
if verbose:
@@ -769,7 +783,7 @@ Description: 2to3
elif aligner == "star" and type == "x2":
assert os.path.exists("SJ.out.tab")
os.system("awk 'BEGIN {OFS=\"\t\"; strChar[0]=\".\"; strChar[1]=\"+\"; strChar[2]=\"-\";} {if($5>0){print $1,$2,$3,strChar[$4]}}' SJ.out.tab > SJ.out.tab.Pass1.sjdb")
-@@ -2518,18 +2518,18 @@
+@@ -2518,18 +2518,18 @@ def calculate_read_cost(single_end,
os.remove(file)
star_index_cmd = "%s/STAR --genomeDir ./ --runMode genomeGenerate --genomeFastaFiles ../../../data/%s.fa --sjdbFileChrStartEnd SJ.out.tab.Pass1.sjdb --sjdbOverhang 99 --runThreadN %d" % (aligner_bin_base, genome, num_threads)
if verbose:
@@ -792,7 +806,7 @@ Description: 2to3
proc = subprocess.Popen(aligner_cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
proc.communicate()
finish_time = datetime.now()
-@@ -2538,7 +2538,7 @@
+@@ -2538,7 +2538,7 @@ def calculate_read_cost(single_end,
if duration < 0.1:
duration = 0.1
if verbose:
@@ -801,7 +815,7 @@ Description: 2to3
os.system("mv Aligned.out.sam %s" % out_fname)
elif aligner == "tophat2":
os.system("samtools sort -n tophat_out/accepted_hits.bam accepted_hits; samtools view -h accepted_hits.bam > %s" % out_fname)
-@@ -2551,7 +2551,7 @@
+@@ -2551,7 +2551,7 @@ def calculate_read_cost(single_end,
os.system("tar cvzf %s.tar.gz %s &> /dev/null" % (out_fname, out_fname))
if runtime_only:
@@ -810,7 +824,7 @@ Description: 2to3
os.chdir("..")
continue
-@@ -2560,7 +2560,7 @@
+@@ -2560,7 +2560,7 @@ def calculate_read_cost(single_end,
pid_list = []
if paired:
if mp_mode:
@@ -819,7 +833,7 @@ Description: 2to3
p = Process(target=extract_pair, args=(out_fname, out_fname2, chr_dic, RNA, aligner, version, repeat_db, repeat_map, debug_dic, i))
pid_list.append(p)
p.start()
-@@ -2570,7 +2570,7 @@
+@@ -2570,7 +2570,7 @@ def calculate_read_cost(single_end,
# merge
os.system("mv %s %s" % (out_fname2 + ".0", out_fname2))
@@ -828,7 +842,7 @@ Description: 2to3
os.system("cat %s >> %s" % (out_fname2 + "." + str(i), out_fname2))
os.system("rm %s" % (out_fname2 + "." + str(i)))
-@@ -2581,7 +2581,7 @@
+@@ -2581,7 +2581,7 @@ def calculate_read_cost(single_end,
else:
if mp_mode:
# Prepare queues
@@ -837,7 +851,7 @@ Description: 2to3
p = Process(target=extract_single, args=(out_fname, out_fname2, chr_dic, aligner, version, repeat_db, repeat_map, debug_dic, i))
pid_list.append(p)
p.start()
-@@ -2592,7 +2592,7 @@
+@@ -2592,7 +2592,7 @@ def calculate_read_cost(single_end,
# merge
os.system("mv %s %s" % (out_fname2 + ".0", out_fname2))
@@ -846,7 +860,7 @@ Description: 2to3
os.system("cat %s >> %s" % (out_fname2 + "." + str(i), out_fname2))
os.system("rm %s" % (out_fname2 + "." + str(i)))
-@@ -2628,25 +2628,25 @@
+@@ -2628,25 +2628,25 @@ def calculate_read_cost(single_end,
assert mapped + unmapped == numreads
if two_step:
@@ -885,7 +899,7 @@ Description: 2to3
if duration > 0.0:
if sql_write and os.path.exists("../" + sql_db_name):
-@@ -2669,14 +2669,14 @@
+@@ -2669,14 +2669,14 @@ def calculate_read_cost(single_end,
os.chdir("..")
@@ -902,15 +916,17 @@ Description: 2to3
if os.path.exists(sql_db_name):
write_analysis_data(sql_db_name, genome, data_base)
---- hisat2.orig/evaluation/simulation/init.py
-+++ hisat2/evaluation/simulation/init.py
+Index: hisat2.git/evaluation/simulation/init.py
+===================================================================
+--- hisat2.git.orig/evaluation/simulation/init.py
++++ hisat2.git/evaluation/simulation/init.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
import sys, os
import string, re
-@@ -50,7 +50,7 @@
+@@ -50,7 +50,7 @@ def extract_pair(RNA):
if p_str not in pair_reported:
pair_reported.add(p_str)
@@ -919,7 +935,7 @@ Description: 2to3
if not me in read_dic:
read_dic[me] = []
-@@ -83,8 +83,7 @@
+@@ -83,8 +83,7 @@ def init_reads(read_dir, RNA):
elif field.startswith("Zs"):
Zs = "\t" + field
@@ -929,7 +945,7 @@ Description: 2to3
sim1_file.close()
sim2_file.close()
-@@ -194,7 +193,7 @@
+@@ -194,7 +193,7 @@ def classify_reads(RNA):
break
assert left_anchor > 0
right_anchor = 0
@@ -938,7 +954,7 @@ Description: 2to3
c = cigars[ci][-1]
c_len = int(cigars[ci][:-1])
if c in "MI":
-@@ -262,7 +261,7 @@
+@@ -262,7 +261,7 @@ def classify_reads(RNA):
if readtype == "all" or readtype == readtype2:
read_ids.add(read_id)
@@ -947,7 +963,7 @@ Description: 2to3
junctions += get_junctions(chr, int(pos), cigar)
if paired:
junctions += get_junctions(chr2, int(pos2), cigar2)
-@@ -285,7 +284,7 @@
+@@ -285,7 +284,7 @@ def classify_reads(RNA):
# write the junctions into a file
type_junction_file = open(type_junction_fname, "w")
for junction in junctions:
@@ -956,7 +972,7 @@ Description: 2to3
type_junction_file.close()
def write_reads(read_fname, type_read_fname):
-@@ -299,7 +298,7 @@
+@@ -299,7 +298,7 @@ def classify_reads(RNA):
write = read_id in read_ids
if write:
@@ -965,7 +981,7 @@ Description: 2to3
read_file.close()
type_read_file.close()
-@@ -323,7 +322,7 @@
+@@ -323,7 +322,7 @@ def init():
not os.path.exists(read_dir_base + read_dir + "/sim_2.fa"):
continue
@@ -974,15 +990,17 @@ Description: 2to3
os.mkdir(read_dir)
os.chdir(read_dir)
---- hisat2.orig/evaluation/tests/HLA_novel/hisatgenotype_locus_prev.py
-+++ hisat2/evaluation/tests/HLA_novel/hisatgenotype_locus_prev.py
+Index: hisat2.git/evaluation/tests/HLA_novel/hisatgenotype_locus_prev.py
+===================================================================
+--- hisat2.git.orig/evaluation/tests/HLA_novel/hisatgenotype_locus_prev.py
++++ hisat2.git/evaluation/tests/HLA_novel/hisatgenotype_locus_prev.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
#
# Copyright 2015, Daehwan Kim <infphilo@gmail.com>
#
-@@ -79,7 +79,7 @@
+@@ -79,7 +79,7 @@ def align_reads(ex_path,
aligner_cmd += ["-1", "%s" % read_fname[0],
"-2", "%s" % read_fname[1]]
if verbose >= 1:
@@ -991,7 +1009,7 @@ Description: 2to3
align_proc = subprocess.Popen(aligner_cmd,
stdout=subprocess.PIPE,
stderr=open("/dev/null", 'w'))
-@@ -116,7 +116,7 @@
+@@ -116,7 +116,7 @@ def align_reads(ex_path,
"""
def normalize(prob):
total = sum(prob.values())
@@ -1000,7 +1018,7 @@ Description: 2to3
prob[allele] = mass / total
-@@ -124,7 +124,7 @@
+@@ -124,7 +124,7 @@ def normalize(prob):
"""
def prob_diff(prob1, prob2):
diff = 0.0
@@ -1009,7 +1027,7 @@ Description: 2to3
if allele in prob2:
diff += abs(prob1[allele] - prob2[allele])
else:
-@@ -153,15 +153,15 @@
+@@ -153,15 +153,15 @@ def single_abundance(Gene_cmpt,
Gene_length):
def normalize2(prob, length):
total = 0
@@ -1028,7 +1046,7 @@ Description: 2to3
alleles = cmpt.split('-')
for allele in alleles:
if allele not in Gene_prob:
-@@ -171,7 +171,7 @@
+@@ -171,7 +171,7 @@ def single_abundance(Gene_cmpt,
normalize(Gene_prob)
def next_prob(Gene_cmpt, Gene_prob, Gene_length):
Gene_prob_next = {}
@@ -1037,7 +1055,7 @@ Description: 2to3
alleles = cmpt.split('-')
alleles_prob = 0.0
for allele in alleles:
-@@ -201,14 +201,14 @@
+@@ -201,14 +201,14 @@ def single_abundance(Gene_cmpt,
Gene_prob_next2 = next_prob(Gene_cmpt, Gene_prob_next, Gene_length)
sum_squared_r, sum_squared_v = 0.0, 0.0
p_r, p_v = {}, {}
@@ -1054,7 +1072,7 @@ Description: 2to3
Gene_prob_next2[a] = max(0.0, Gene_prob[a] - 2 * gamma * p_r[a] + gamma * gamma * p_v[a]);
Gene_prob_next = next_prob(Gene_cmpt, Gene_prob_next2, Gene_length)
-@@ -219,17 +219,17 @@
+@@ -219,17 +219,17 @@ def single_abundance(Gene_cmpt,
if iter >= 10:
Gene_prob2 = {}
avg_prob = sum(Gene_prob.values()) / len(Gene_prob)
@@ -1076,7 +1094,7 @@ Description: 2to3
iter += 1
-@@ -241,7 +241,7 @@
+@@ -241,7 +241,7 @@ def single_abundance(Gene_cmpt,
# normalize(Gene_prob)
normalize2(Gene_prob, Gene_length)
@@ -1085,7 +1103,7 @@ Description: 2to3
Gene_prob = sorted(Gene_prob, cmp=Gene_prob_cmp)
return Gene_prob
-@@ -288,7 +288,7 @@
+@@ -288,7 +288,7 @@ Report variant IDs whose var is within e
"""
def get_exonic_vars(Vars, exons):
vars = set()
@@ -1094,7 +1112,7 @@ Description: 2to3
var_type, var_left, var_data = var
var_right = var_left
if var_type == "deletion":
-@@ -305,7 +305,7 @@
+@@ -305,7 +305,7 @@ Get representative alleles among those t
"""
def get_rep_alleles(Links, exon_vars):
allele_vars = {}
@@ -1103,7 +1121,7 @@ Description: 2to3
if var not in exon_vars:
continue
for allele in alleles:
-@@ -314,7 +314,7 @@
+@@ -314,7 +314,7 @@ def get_rep_alleles(Links, exon_vars):
allele_vars[allele].add(var)
allele_groups = {}
@@ -1112,7 +1130,7 @@ Description: 2to3
vars = '-'.join(vars)
if vars not in allele_groups:
allele_groups[vars] = []
-@@ -322,7 +322,7 @@
+@@ -322,7 +322,7 @@ def get_rep_alleles(Links, exon_vars):
allele_reps = {} # allele representatives
allele_rep_groups = {} # allele groups by allele representatives
@@ -1121,7 +1139,7 @@ Description: 2to3
assert len(allele_members) > 0
allele_rep = allele_members[0]
allele_rep_groups[allele_rep] = allele_members
-@@ -402,9 +402,9 @@
+@@ -402,9 +402,9 @@ def get_alternatives(ref_seq, Vars, Var_
if j_type == "single" and j_pos == latest_pos - 1:
j_cmp_pos = j_pos - del_len + other_del_len
if debug:
@@ -1134,7 +1152,7 @@ Description: 2to3
if j_data == ref_seq[j_cmp_pos]:
add_alt(Alts, alt_list, var_id, j_id)
latest_pos = j_pos
-@@ -465,9 +465,9 @@
+@@ -465,9 +465,9 @@ def get_alternatives(ref_seq, Vars, Var_
# DK - debugging purposes
if debug:
pos2_ = latest_pos + 1 + del_len - other_del_len
@@ -1147,7 +1165,7 @@ Description: 2to3
latest_pos += 1
add_alt(Alts, alt_list, var_id, str(latest_pos))
-@@ -479,9 +479,9 @@
+@@ -479,9 +479,9 @@ def get_alternatives(ref_seq, Vars, Var_
if j_type == "single" and j_pos == latest_pos + 1:
j_cmp_pos = j_pos + del_len - other_del_len
if debug:
@@ -1160,7 +1178,7 @@ Description: 2to3
if j_data == ref_seq[j_cmp_pos]:
add_alt(Alts, alt_list, var_id, j_id)
-@@ -542,7 +542,7 @@
+@@ -542,7 +542,7 @@ def get_alternatives(ref_seq, Vars, Var_
continue
debug = (var_id == "hv454a")
if debug:
@@ -1169,7 +1187,7 @@ Description: 2to3
alt_list = []
var_j = lower_bound(Var_list, var_pos + del_len - 1)
-@@ -578,7 +578,7 @@
+@@ -578,7 +578,7 @@ def get_alternatives(ref_seq, Vars, Var_
debug)
if debug:
@@ -1178,7 +1196,7 @@ Description: 2to3
sys.exit(1)
def assert_print_alts(Alts, dir):
-@@ -640,8 +640,8 @@
+@@ -640,8 +640,8 @@ def get_alternatives(ref_seq, Vars, Var_
return seq, seq_left, seq_right
@@ -1189,7 +1207,7 @@ Description: 2to3
out_str = "\t\t"
alt_list1 = alt_list1.split('-')
for i in range(len(alt_list1)):
-@@ -664,7 +664,7 @@
+@@ -664,7 +664,7 @@ def get_alternatives(ref_seq, Vars, Var_
if j + 1 < len(alt_list3):
out_str += ", "
out_str += "]"
@@ -1198,7 +1216,7 @@ Description: 2to3
for i in range(len(alt_list2)):
alt_list3 = alt_list2[i]
-@@ -702,9 +702,9 @@
+@@ -702,9 +702,9 @@ def get_alternatives(ref_seq, Vars, Var_
else:
seq2 += ref_seq[seq2_right:seq2_right+len_diff]
if verbose >= 3:
@@ -1211,7 +1229,7 @@ Description: 2to3
assert seq1 == seq2
assert_print_alts(Alts_left, "left")
-@@ -739,7 +739,7 @@
+@@ -739,7 +739,7 @@ def identify_ambigious_diffs(Vars,
# Left direction
id_str = var_id
total_del_len = length if type == "deletion" else 0
@@ -1220,7 +1238,7 @@ Description: 2to3
cmp_j = cmp_list[j]
j_type, j_pos, j_len = cmp_j[:3]
if j_type != "match":
-@@ -792,16 +792,16 @@
+@@ -792,16 +792,16 @@ def identify_ambigious_diffs(Vars,
alt_left_pos += alt_total_del_len
if left_pos >= alt_left_pos:
if verbose >= 2:
@@ -1242,7 +1260,7 @@ Description: 2to3
# sys.exit(1)
# Right direction
-@@ -824,7 +824,7 @@
+@@ -824,7 +824,7 @@ def identify_ambigious_diffs(Vars,
# DK - debugging purposes
if debug:
@@ -1251,7 +1269,7 @@ Description: 2to3
if id_str in Alts_right:
orig_alts = id_str.split('-')
-@@ -877,8 +877,8 @@
+@@ -877,8 +877,8 @@ def identify_ambigious_diffs(Vars,
if right_pos <= alt_right_pos:
if verbose >= 2:
@@ -1262,7 +1280,7 @@ Description: 2to3
cmp_right = i - 1
break
i += 1
-@@ -971,7 +971,7 @@
+@@ -971,7 +971,7 @@ def get_mpileup(alignview_cmd,
num_nt = sum(nt_dic.values())
nt_set = []
if num_nt >= 20:
@@ -1271,7 +1289,7 @@ Description: 2to3
if nt not in "ACGT":
continue
if count >= num_nt * 0.2 or count >= 7:
-@@ -980,7 +980,7 @@
+@@ -980,7 +980,7 @@ def get_mpileup(alignview_cmd,
# Sort variants
var_list = [[] for i in range(len(mpileup))]
@@ -1280,7 +1298,7 @@ Description: 2to3
var_type, var_pos, var_data = value
assert var_pos < len(var_list)
var_list[var_pos].append([var_id, var_type, var_data])
-@@ -991,7 +991,7 @@
+@@ -991,7 +991,7 @@ def get_mpileup(alignview_cmd,
nt_dic = mpileup[i][1]
ref_nt = ref_seq[i]
new_nt_dic = {}
@@ -1289,7 +1307,7 @@ Description: 2to3
var_id = ""
if nt == 'D':
if i <= skip_i:
-@@ -1041,8 +1041,8 @@
+@@ -1041,8 +1041,8 @@ def error_correct(ref_seq,
cmp_list,
debug = False):
if debug:
@@ -1300,7 +1318,7 @@ Description: 2to3
i = 0
while i < len(cmp_list):
-@@ -1089,7 +1089,7 @@
+@@ -1089,7 +1089,7 @@ def error_correct(ref_seq,
nt_set = mpileup[left][0]
if debug:
@@ -1309,7 +1327,7 @@ Description: 2to3
if len(nt_set) > 0 and read_bp not in nt_set:
read_bp = 'N' if len(nt_set) > 1 else nt_set[0]
-@@ -1113,8 +1113,8 @@
+@@ -1113,8 +1113,8 @@ def error_correct(ref_seq,
var_idx += 1
if debug:
@@ -1320,7 +1338,7 @@ Description: 2to3
read_pos += length
i += 1
-@@ -1132,8 +1132,8 @@
+@@ -1132,8 +1132,8 @@ def error_correct(ref_seq,
i += 1
if debug:
@@ -1331,7 +1349,7 @@ Description: 2to3
return cmp_list, read_seq
-@@ -1267,9 +1267,9 @@
+@@ -1267,9 +1267,9 @@ def typing(ex_path,
for aligner, index_type in aligners:
for f_ in [sys.stderr, report_file]:
if index_type == "graph":
@@ -1343,7 +1361,7 @@ Description: 2to3
remove_alignment_file = False
if alignment_fname == "":
-@@ -1364,7 +1364,7 @@
+@@ -1364,7 +1364,7 @@ def typing(ex_path,
# List of nodes that represent alleles
allele_vars = {}
@@ -1352,7 +1370,7 @@ Description: 2to3
for allele_id in allele_list:
if allele_id not in Genes[gene]:
continue
-@@ -1480,8 +1480,8 @@
+@@ -1480,8 +1480,8 @@ def typing(ex_path,
# Unalined?
if flag & 0x4 != 0:
if simulation and verbose >= 2:
@@ -1363,7 +1381,7 @@ Description: 2to3
continue
# Concordantly mapped?
-@@ -1686,7 +1686,7 @@
+@@ -1686,7 +1686,7 @@ def typing(ex_path,
# Check if this deletion is artificial alignment
assert right_pos < mpileup
del_count, nt_count = 0, 0
@@ -1372,7 +1390,7 @@ Description: 2to3
count = value[0]
if nt == 'D':
del_count += count
-@@ -1775,7 +1775,7 @@
+@@ -1775,7 +1775,7 @@ def typing(ex_path,
def add_stat(Gene_cmpt, Gene_counts, Gene_count_per_read, include_alleles = set()):
max_count = max(Gene_count_per_read.values())
cur_cmpt = set()
@@ -1381,7 +1399,7 @@ Description: 2to3
if count < max_count:
continue
-@@ -1796,7 +1796,7 @@
+@@ -1796,7 +1796,7 @@ def typing(ex_path,
# alleles = ["A*24:36N", "A*24:359N"]
allele1_found, allele2_found = False, False
if alleles[0] != "":
@@ -1390,7 +1408,7 @@ Description: 2to3
if count < max_count:
continue
if allele == alleles[0]:
-@@ -1804,13 +1804,13 @@
+@@ -1804,13 +1804,13 @@ def typing(ex_path,
elif allele == alleles[1]:
allele2_found = True
if allele1_found != allele2_found:
@@ -1409,7 +1427,7 @@ Description: 2to3
cur_cmpt = sorted(list(cur_cmpt))
cur_cmpt = '-'.join(cur_cmpt)
-@@ -1832,7 +1832,7 @@
+@@ -1832,7 +1832,7 @@ def typing(ex_path,
if debug_read_count == debug_max_read_count and \
Gene_count_per_read["A*24:02:01:02L"] < debug_max_read_count and \
Gene_count_per_read["A*01:01:01:01"] < debug_max_read_count:
@@ -1418,7 +1436,7 @@ Description: 2to3
None
if prev_read_id == "HSQ1008:175:C0JVFACXX:7:1208:5604:41201":
-@@ -1861,9 +1861,9 @@
+@@ -1861,9 +1861,9 @@ def typing(ex_path,
base_fname == "hla":
cur_cmpt = cur_cmpt.split('-')
if not(set(cur_cmpt) & set(test_Gene_names)):
@@ -1430,7 +1448,7 @@ Description: 2to3
prev_lines = []
-@@ -1880,18 +1880,18 @@
+@@ -1880,18 +1880,18 @@ def typing(ex_path,
alleles = Links[var_id]
if verbose >= 2:
if add > 0 and not (set(alleles) & debug_allele_names):
@@ -1455,7 +1473,7 @@ Description: 2to3
if var_id == "unknown" or var_id.startswith("nv"):
continue
var_type, var_pos, var_data = data
-@@ -1924,14 +1924,14 @@
+@@ -1924,14 +1924,14 @@ def typing(ex_path,
DK_debug = False
if orig_read_id == "a46|L_451_88M12D12M_88|D|hv2":
DK_debug = True
@@ -1476,7 +1494,7 @@ Description: 2to3
var_type, var_pos, var_data = data
if var_type != "deletion":
continue
-@@ -2002,8 +2002,8 @@
+@@ -2002,8 +2002,8 @@ def typing(ex_path,
positive_vars.add(var_id)
if read_id == "HSQ1008:175:C0JVFACXX:6:2207:13481:60924" and False:
@@ -1487,7 +1505,7 @@ Description: 2to3
cmp_MD += ("%d%s" % (MD_match_len, ref_seq[ref_pos]))
-@@ -2069,15 +2069,15 @@
+@@ -2069,15 +2069,15 @@ def typing(ex_path,
if read_pos != len(read_seq) or \
cmp_cigar_str != cigar_str:
# cmp_MD != MD: # Disabled due to error correction
@@ -1508,7 +1526,7 @@ Description: 2to3
# Node
if assembly:
-@@ -2100,8 +2100,8 @@
+@@ -2100,8 +2100,8 @@ def typing(ex_path,
add_count(Gene_gen_count_per_read, positive_var, 1)
if read_id == "HSQ1008:175:C0JVFACXX:6:2207:13481:60924" and False:
@@ -1519,7 +1537,7 @@ Description: 2to3
for negative_var in negative_vars:
-@@ -2118,7 +2118,7 @@
+@@ -2118,7 +2118,7 @@ def typing(ex_path,
continue
for f_ in [sys.stderr, report_file]:
@@ -1528,7 +1546,7 @@ Description: 2to3
if prev_read_id != None:
if base_fname == "hla":
-@@ -2182,7 +2182,7 @@
+@@ -2182,7 +2182,7 @@ def typing(ex_path,
if base_fname != "hla":
Gene_counts = Gene_gen_counts
@@ -1537,7 +1555,7 @@ Description: 2to3
def Gene_count_cmp(a, b):
if a[1] != b[1]:
return b[1] - a[1]
-@@ -2199,7 +2199,7 @@
+@@ -2199,7 +2199,7 @@ def typing(ex_path,
for test_Gene_name in test_Gene_names:
if count[0] == test_Gene_name:
for f_ in [sys.stderr, report_file]:
@@ -1546,7 +1564,7 @@ Description: 2to3
found = True
"""
if count_i > 0 and Gene_counts[0][1] > count[1]:
-@@ -2210,14 +2210,14 @@
+@@ -2210,14 +2210,14 @@ def typing(ex_path,
"""
if count_i < 5 and not found:
for f_ in [sys.stderr, report_file]:
@@ -1564,7 +1582,7 @@ Description: 2to3
# Calculate the abundance of representative alleles on exonic sequences
if base_fname == "hla":
-@@ -2238,7 +2238,7 @@
+@@ -2238,7 +2238,7 @@ def typing(ex_path,
if len(gen_alleles) > 0:
Gene_gen_cmpt2 = {}
@@ -1573,7 +1591,7 @@ Description: 2to3
cmpt2 = []
for allele in cmpt.split('-'):
if allele in gen_alleles:
-@@ -2262,7 +2262,7 @@
+@@ -2262,7 +2262,7 @@ def typing(ex_path,
Gene_combined_prob[allele] = prob
for allele, prob in Gene_gen_prob:
Gene_combined_prob[allele] = prob * gen_prob_sum
@@ -1582,7 +1600,7 @@ Description: 2to3
Gene_prob = sorted(Gene_prob, cmp=Gene_prob_cmp)
else:
Gene_prob = single_abundance(Gene_gen_cmpt, Gene_lengths[gene])
-@@ -2302,8 +2302,8 @@
+@@ -2302,8 +2302,8 @@ def typing(ex_path,
ref_allele)
def get_best_alleles(left, right, vars):
@@ -1593,7 +1611,7 @@ Description: 2to3
tmp_vars = allele_node.get_var_ids(left, right)
tmp_common = len(set(vars) & set(tmp_vars))
tmp_common -= len(set(vars) | set(tmp_vars))
-@@ -2315,13 +2315,13 @@
+@@ -2315,13 +2315,13 @@ def typing(ex_path,
return max_alleles
for run, plus, minus in run_alignments:
@@ -1612,7 +1630,7 @@ Description: 2to3
assert False
-@@ -2342,7 +2342,7 @@
+@@ -2342,7 +2342,7 @@ def typing(ex_path,
# Compare two alleles
if simulation and len(test_Gene_names) == 2:
allele_name1, allele_name2 = test_Gene_names
@@ -1621,7 +1639,7 @@ Description: 2to3
asm_graph.print_node_comparison(asm_graph.true_allele_nodes)
def compare_alleles(vars1, vars2, print_output = True):
-@@ -2374,8 +2374,8 @@
+@@ -2374,8 +2374,8 @@ def typing(ex_path,
skip = False
if print_output:
if cmp_var_in_exon:
@@ -1632,7 +1650,7 @@ Description: 2to3
var_i += 1; var_j += 1
continue
if cmp_var[1] <= node_var[1]:
-@@ -2384,42 +2384,42 @@
+@@ -2384,42 +2384,42 @@ def typing(ex_path,
if print_output:
if cmp_var_in_exon:
for f_ in [sys.stderr, report_file]:
@@ -1684,7 +1702,7 @@ Description: 2to3
for cmp_Gene_name in cmp_Gene_names:
tmp_vars = allele_nodes[cmp_Gene_name].get_var_ids(node.left, node.right)
tmp_common = len(set(node_vars) & set(tmp_vars))
-@@ -2432,19 +2432,19 @@
+@@ -2432,19 +2432,19 @@ def typing(ex_path,
for allele_name, cmp_vars in alleles:
for f_ in [sys.stderr, report_file]:
@@ -1710,7 +1728,7 @@ Description: 2to3
vars2 = set(vars2)
tmp_common = len(vars & vars2) - len(vars | vars2)
if tmp_common > max_common:
-@@ -2454,21 +2454,21 @@
+@@ -2454,21 +2454,21 @@ def typing(ex_path,
max_allele_names.append(allele_name)
for f_ in [sys.stderr, report_file]:
@@ -1737,7 +1755,7 @@ Description: 2to3
"""
allele_exon_vars = {}
-@@ -2520,7 +2520,7 @@
+@@ -2520,7 +2520,7 @@ def typing(ex_path,
else:
break
for f_ in [sys.stderr, report_file]:
@@ -1746,7 +1764,7 @@ Description: 2to3
if rank_i < len(success):
success[rank_i] = True
found_list[name_i] = True
-@@ -2530,15 +2530,15 @@
+@@ -2530,15 +2530,15 @@ def typing(ex_path,
break
if not found:
for f_ in [sys.stderr, report_file]:
@@ -1765,7 +1783,7 @@ Description: 2to3
if simulation and not False in success:
aligner_type = "%s %s" % (aligner, index_type)
-@@ -2589,7 +2589,7 @@
+@@ -2589,7 +2589,7 @@ def read_Gene_vars(fname):
Vars[gene][var_id] = [var_type, pos - left, data]
Var_list[gene].append([pos - left, var_id])
@@ -1774,7 +1792,7 @@ Description: 2to3
Var_list[gene] = sorted(in_var_list)
return Vars, Var_list
-@@ -2707,7 +2707,7 @@
+@@ -2707,7 +2707,7 @@ def test_Gene_genotyping(base_fname,
base_fname + ".link"]
if verbose >= 1:
@@ -1783,7 +1801,7 @@ Description: 2to3
if not typing_common.check_files(Gene_fnames):
extract_hla_script = os.path.join(ex_path, "hisatgenotype_extract_vars.py")
-@@ -2731,12 +2731,12 @@
+@@ -2731,12 +2731,12 @@ def test_Gene_genotyping(base_fname,
# DK - debugging purposes
# extract_cmd += ["--ext-seq", "300"]
if verbose >= 1:
@@ -1798,7 +1816,7 @@ Description: 2to3
sys.exit(1)
for aligner, index_type in aligners:
-@@ -2753,11 +2753,11 @@
+@@ -2753,11 +2753,11 @@ def test_Gene_genotyping(base_fname,
"%s_backbone.fa" % base_fname,
"%s.graph" % base_fname]
if verbose >= 1:
@@ -1812,7 +1830,7 @@ Description: 2to3
sys.exit(1)
# Build HISAT2 linear indexes based on the above information
else:
-@@ -2771,7 +2771,7 @@
+@@ -2771,7 +2771,7 @@ def test_Gene_genotyping(base_fname,
proc = subprocess.Popen(build_cmd, stdout=open("/dev/null", 'w'), stderr=open("/dev/null", 'w'))
proc.communicate()
if not typing_common.check_files(Gene_hisat2_linear_index_fnames):
@@ -1821,7 +1839,7 @@ Description: 2to3
sys.exit(1)
else:
assert aligner == "bowtie2" and index_type == "linear"
-@@ -2785,7 +2785,7 @@
+@@ -2785,7 +2785,7 @@ def test_Gene_genotyping(base_fname,
proc = subprocess.Popen(build_cmd, stdout=open("/dev/null", 'w'))
proc.communicate()
if not typing_common.check_files(Gene_bowtie2_index_fnames):
@@ -1830,7 +1848,7 @@ Description: 2to3
sys.exit(1)
# Read partial alleles from hla.data (temporary)
-@@ -2815,21 +2815,21 @@
+@@ -2815,21 +2815,21 @@ def test_Gene_genotyping(base_fname,
refGene_loci[Gene_gene] = [Gene_name, chr, left, right, exons]
Genes = {}
if len(locus_list) == 0:
@@ -1856,7 +1874,7 @@ Description: 2to3
Gene_lengths[Gene_gene][allele_name] = len(seq)
# Read HLA variants, and link information
-@@ -2886,7 +2886,7 @@
+@@ -2886,7 +2886,7 @@ def test_Gene_genotyping(base_fname,
if str(test_i + 1) not in test_ids:
continue
@@ -1865,7 +1883,7 @@ Description: 2to3
test_locus_list = test_list[test_i]
num_frag_list = typing_common.simulate_reads(Genes,
base_fname,
-@@ -2910,7 +2910,7 @@
+@@ -2910,7 +2910,7 @@ def test_Gene_genotyping(base_fname,
gene = test_Gene_name.split('*')[0]
test_Gene_seq = Genes[gene][test_Gene_name]
seq_type = "partial" if test_Gene_name in partial_alleles else "full"
@@ -1874,7 +1892,7 @@ Description: 2to3
if "single-end" in debug_instr:
read_fname = ["%s_input_1.fa" % base_fname]
-@@ -2948,20 +2948,20 @@
+@@ -2948,20 +2948,20 @@ def test_Gene_genotyping(base_fname,
best_alleles,
verbose)
@@ -1900,7 +1918,7 @@ Description: 2to3
fastq = True
typing(ex_path,
simulation,
-@@ -3125,7 +3125,7 @@
+@@ -3125,7 +3125,7 @@ if __name__ == '__main__':
else:
locus_list = args.locus_list.split(',')
if args.aligners == "":
@@ -1909,7 +1927,7 @@ Description: 2to3
sys.exit(1)
args.aligners = args.aligners.split(',')
for i in range(len(args.aligners)):
-@@ -3136,7 +3136,7 @@
+@@ -3136,7 +3136,7 @@ if __name__ == '__main__':
args.read_fname = []
if args.alignment_fname != "" and \
not os.path.exists(args.alignment_fname):
@@ -1918,7 +1936,7 @@ Description: 2to3
sys.exit(1)
if args.verbose and args.verbose_level == 0:
-@@ -3154,10 +3154,10 @@
+@@ -3154,10 +3154,10 @@ if __name__ == '__main__':
debug[item] = 1
if not args.partial:
@@ -1931,7 +1949,7 @@ Description: 2to3
skip_fragment_regions = []
if args.skip_fragment_regions != "":
-@@ -3178,7 +3178,7 @@
+@@ -3178,7 +3178,7 @@ if __name__ == '__main__':
if args.stranded_seq != "":
stranded_seq = args.stranded_seq.split(',')
if len(stranded_seq) != 2:
@@ -1940,15 +1958,17 @@ Description: 2to3
sys.exit(1)
else:
stranded_seq = []
---- hisat2.orig/evaluation/tests/genotype_genome/hisatgenotype_prev.py
-+++ hisat2/evaluation/tests/genotype_genome/hisatgenotype_prev.py
+Index: hisat2.git/evaluation/tests/genotype_genome/hisatgenotype_prev.py
+===================================================================
+--- hisat2.git.orig/evaluation/tests/genotype_genome/hisatgenotype_prev.py
++++ hisat2.git/evaluation/tests/genotype_genome/hisatgenotype_prev.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
#
# Copyright 2016, Daehwan Kim <infphilo@gmail.com>
-@@ -51,9 +51,9 @@
+@@ -51,9 +51,9 @@ def align_reads(base_fname,
aligner_cmd += ["-1", read_fnames[0],
"-2", read_fnames[1]]
@@ -1960,7 +1980,7 @@ Description: 2to3
align_proc = subprocess.Popen(aligner_cmd,
stdout=subprocess.PIPE,
-@@ -69,7 +69,7 @@
+@@ -69,7 +69,7 @@ def align_reads(base_fname,
stderr=open("/dev/null", 'w'))
sambam_proc.communicate()
@@ -1969,7 +1989,7 @@ Description: 2to3
bamsort_cmd = ["samtools",
"sort",
"--threads", str(threads),
-@@ -79,7 +79,7 @@
+@@ -79,7 +79,7 @@ def align_reads(base_fname,
stderr=open("/dev/null", 'w'))
bamsort_proc.communicate()
@@ -1978,7 +1998,7 @@ Description: 2to3
bamindex_cmd = ["samtools",
"index",
-@@ -116,9 +116,9 @@
+@@ -116,9 +116,9 @@ def genotype(base_fname,
# hisat2 graph index files
genotype_fnames += ["%s.%d.ht2" % (base_fname, i+1) for i in range(8)]
if not typing_common.check_files(genotype_fnames):
@@ -1990,7 +2010,7 @@ Description: 2to3
sys.exit(1)
# Align reads, and sort the alignments into a BAM file
-@@ -165,11 +165,11 @@
+@@ -165,11 +165,11 @@ def genotype(base_fname,
# gene alleles
allele_names = {}
@@ -2004,7 +2024,7 @@ Description: 2to3
allele_name1 = allele_name.split('*')[0]
if gene_name2 == allele_name1:
allele_names[gene_name].append(allele_name)
-@@ -196,7 +196,7 @@
+@@ -196,7 +196,7 @@ def genotype(base_fname,
Vars[gene_name][var_id] = [var_type, pos, data]
Var_list[gene_name].append([pos, var_id])
@@ -2013,7 +2033,7 @@ Description: 2to3
Var_list[gene_name] = sorted(in_var_list)
def lower_bound(Var_list, pos):
low, high = 0, len(Var_list)
-@@ -233,7 +233,7 @@
+@@ -233,7 +233,7 @@ def genotype(base_fname,
for test_i in range(len(test_list)):
test_HLA_list = test_list[test_i]
for test_HLA_names in test_HLA_list:
@@ -2022,7 +2042,7 @@ Description: 2to3
for gene in test_HLA_names:
ref_allele = genes[gene]
ref_seq = gene_seqs[gene]
-@@ -303,8 +303,8 @@
+@@ -303,8 +303,8 @@ def genotype(base_fname,
debug = False
if read_id in ["2339"] and False:
debug = True
@@ -2033,7 +2053,7 @@ Description: 2to3
vars = []
if Zs:
-@@ -394,7 +394,7 @@
+@@ -394,7 +394,7 @@ def genotype(base_fname,
def add_stat(HLA_cmpt, HLA_counts, HLA_count_per_read, exon = True):
max_count = max(HLA_count_per_read.values())
cur_cmpt = set()
@@ -2042,7 +2062,7 @@ Description: 2to3
if count < max_count:
continue
"""
-@@ -414,7 +414,7 @@
+@@ -414,7 +414,7 @@ def genotype(base_fname,
alleles = ["", ""]
# alleles = ["B*40:304", "B*40:02:01"]
allele1_found, allele2_found = False, False
@@ -2051,7 +2071,7 @@ Description: 2to3
if count < max_count:
continue
if allele == alleles[0]:
-@@ -422,13 +422,13 @@
+@@ -422,13 +422,13 @@ def genotype(base_fname,
elif allele == alleles[1]:
allele2_found = True
if allele1_found != allele2_found:
@@ -2070,7 +2090,7 @@ Description: 2to3
cur_cmpt = sorted(list(cur_cmpt))
cur_cmpt = '-'.join(cur_cmpt)
-@@ -462,11 +462,11 @@
+@@ -462,11 +462,11 @@ def genotype(base_fname,
# daehwan - for debugging purposes
if debug:
if allele in ["DQA1*05:05:01:01", "DQA1*05:05:01:02"]:
@@ -2084,7 +2104,7 @@ Description: 2to3
var_type, var_pos, var_data = data
if var_type != "deletion":
continue
-@@ -490,13 +490,13 @@
+@@ -490,13 +490,13 @@ def genotype(base_fname,
add_count(var_id, -1)
# daehwan - for debugging purposes
if debug:
@@ -2100,7 +2120,7 @@ Description: 2to3
# Check if this might be one of the two tandem repeats (the same left coordinate)
cmp_left, cmp_right = cmp[1], cmp[1] + cmp[2]
test1_seq1 = ref_seq[cmp_left-base_locus:cmp_right-base_locus]
-@@ -510,7 +510,7 @@
+@@ -510,7 +510,7 @@ def genotype(base_fname,
add_count(var_id, -1)
else:
if debug:
@@ -2109,7 +2129,7 @@ Description: 2to3
add_count(var_id, -1)
var_idx += 1
-@@ -531,7 +531,7 @@
+@@ -531,7 +531,7 @@ def genotype(base_fname,
if var_data == read_base:
# daehwan - for debugging purposes
if debug:
@@ -2118,7 +2138,7 @@ Description: 2to3
# daehwan - for debugging purposes
if False:
-@@ -553,8 +553,8 @@
+@@ -553,8 +553,8 @@ def genotype(base_fname,
var_idx = lower_bound(Var_list[gene], ref_pos)
# daehwan - for debugging purposes
if debug:
@@ -2129,7 +2149,7 @@ Description: 2to3
# sys.exit(1)
while var_idx < len(Var_list[gene]):
var_pos, var_id = Var_list[gene][var_idx]
-@@ -566,7 +566,7 @@
+@@ -566,7 +566,7 @@ def genotype(base_fname,
if var_data == ins_seq:
# daehwan - for debugging purposes
if debug:
@@ -2138,7 +2158,7 @@ Description: 2to3
add_count(var_id, 1)
var_idx += 1
-@@ -602,8 +602,8 @@
+@@ -602,8 +602,8 @@ def genotype(base_fname,
var_len = int(var_data)
if var_len == length:
if debug:
@@ -2149,7 +2169,7 @@ Description: 2to3
add_count(var_id, 1)
var_idx += 1
-@@ -634,9 +634,9 @@
+@@ -634,9 +634,9 @@ def genotype(base_fname,
if read_pos != len(read_seq) or \
cmp_cigar_str != cigar_str or \
cmp_MD != MD:
@@ -2162,7 +2182,7 @@ Description: 2to3
assert False
prev_read_id = read_id
-@@ -648,7 +648,7 @@
+@@ -648,7 +648,7 @@ def genotype(base_fname,
if prev_read_id != None:
add_stat(HLA_cmpt, HLA_counts, HLA_count_per_read)
@@ -2171,7 +2191,7 @@ Description: 2to3
def HLA_count_cmp(a, b):
if a[1] != b[1]:
return b[1] - a[1]
-@@ -660,28 +660,28 @@
+@@ -660,28 +660,28 @@ def genotype(base_fname,
HLA_counts = sorted(HLA_counts, cmp=HLA_count_cmp)
for count_i in range(len(HLA_counts)):
count = HLA_counts[count_i]
@@ -2206,7 +2226,7 @@ Description: 2to3
if allele in prob2:
diff += abs(prob1[allele] - prob2[allele])
else:
-@@ -701,7 +701,7 @@
+@@ -701,7 +701,7 @@ def genotype(base_fname,
return 1
HLA_prob, HLA_prob_next = {}, {}
@@ -2215,7 +2235,7 @@ Description: 2to3
alleles = cmpt.split('-')
for allele in alleles:
if allele not in HLA_prob:
-@@ -718,7 +718,7 @@
+@@ -718,7 +718,7 @@ def genotype(base_fname,
normalize(HLA_prob)
def next_prob(HLA_cmpt, HLA_prob, HLA_length):
HLA_prob_next = {}
@@ -2224,7 +2244,7 @@ Description: 2to3
alleles = cmpt.split('-')
alleles_prob = 0.0
for allele in alleles:
-@@ -745,17 +745,17 @@
+@@ -745,17 +745,17 @@ def genotype(base_fname,
HLA_prob[allele] /= float(allele_len)
normalize(HLA_prob)
"""
@@ -2245,7 +2265,7 @@ Description: 2to3
"""
if len(test_HLA_names) == 2:
-@@ -963,15 +963,15 @@
+@@ -963,15 +963,15 @@ def genotype(base_fname,
if cigar_op in "MIS":
read_pos += length
@@ -2264,7 +2284,7 @@ Description: 2to3
-@@ -1038,7 +1038,7 @@
+@@ -1038,7 +1038,7 @@ if __name__ == '__main__':
read_fnames = [args.read_fname_U]
else:
if args.read_fname_1 == "" or args.read_fname_2 == "":
@@ -2273,15 +2293,17 @@ Description: 2to3
sys.exit(1)
read_fnames = [args.read_fname_1, args.read_fname_2]
---- hisat2.orig/evaluation/tests/genotype_genome/paper_sensitivity/sensitivity.py
-+++ hisat2/evaluation/tests/genotype_genome/paper_sensitivity/sensitivity.py
+Index: hisat2.git/evaluation/tests/genotype_genome/paper_sensitivity/sensitivity.py
+===================================================================
+--- hisat2.git.orig/evaluation/tests/genotype_genome/paper_sensitivity/sensitivity.py
++++ hisat2.git/evaluation/tests/genotype_genome/paper_sensitivity/sensitivity.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
#
# Copyright 2017, Daehwan Kim <infphilo@gmail.com>
#
-@@ -69,7 +69,7 @@
+@@ -69,7 +69,7 @@ if __name__ == '__main__':
aligner_cmd += ["-1", fq_fnames[0],
"-2", fq_fnames[1]]
# print >> sys.stderr, "Running:", ' '.join(aligner_cmd)
@@ -2290,7 +2312,7 @@ Description: 2to3
align_proc = subprocess.Popen(aligner_cmd,
stdout=subprocess.PIPE,
stderr=open("/dev/null", 'w'))
-@@ -113,7 +113,7 @@
+@@ -113,7 +113,7 @@ if __name__ == '__main__':
if ((aligner == "hisat2" and NH == 1) or (aligner == "bowtie2" and AS > XS and read1_first if flag & 0x40 else read2_first)):
if chr in region_loci:
@@ -2299,22 +2321,24 @@ Description: 2to3
_, _, loci_left, loci_right = loci
# there might be a different candidate region for each of left and right reads
if pos >= loci_left and pos < loci_right:
-@@ -136,5 +136,5 @@
+@@ -136,5 +136,5 @@ if __name__ == '__main__':
gene = "HLA-" + gene
if gene not in region_count:
continue
- print "\t%s pair: %d, left+right: %d" % (gene, region_count[gene], region_read1_count[gene] + region_read2_count[gene])
+ print("\t%s pair: %d, left+right: %d" % (gene, region_count[gene], region_read1_count[gene] + region_read2_count[gene]))
---- hisat2.orig/evaluation/tests/one_snp_test/evaluate_one_snp_reads.py
-+++ hisat2/evaluation/tests/one_snp_test/evaluate_one_snp_reads.py
+Index: hisat2.git/evaluation/tests/one_snp_test/evaluate_one_snp_reads.py
+===================================================================
+--- hisat2.git.orig/evaluation/tests/one_snp_test/evaluate_one_snp_reads.py
++++ hisat2.git/evaluation/tests/one_snp_test/evaluate_one_snp_reads.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
import sys, os, subprocess
import multiprocessing
-@@ -154,7 +154,7 @@
+@@ -154,7 +154,7 @@ def evaluate(read_fname,
if aligner == "hisat2" and index_type != "":
aligner_name += ("_" + index_type)
two_step = (aligner == "tophat2" or type == "x2" or (aligner in ["hisat2", "hisat"] and type == ""))
@@ -2323,7 +2347,7 @@ Description: 2to3
if paired:
aligner_dir = aligner_name + "_paired"
else:
-@@ -170,7 +170,7 @@
+@@ -170,7 +170,7 @@ def evaluate(read_fname,
aligner_cmd = get_aligner_cmd(aligner, type, index_type, version, "../" + type_read1_fname, "../" + type_read2_fname, out_fname)
start_time = datetime.now()
if verbose:
@@ -2332,7 +2356,7 @@ Description: 2to3
if aligner in ["hisat2", "hisat", "bowtie", "bowtie2", "gsnap", "bwa"]:
proc = subprocess.Popen(aligner_cmd, stdout=open(out_fname, "w"), stderr=subprocess.PIPE)
else:
-@@ -180,7 +180,7 @@
+@@ -180,7 +180,7 @@ def evaluate(read_fname,
duration = finish_time - start_time
duration = duration.total_seconds()
if verbose:
@@ -2341,7 +2365,7 @@ Description: 2to3
assert os.path.exists(out_fname)
correct_reads, correct_multi_reads, num_reads = 0, 0, 0
-@@ -205,8 +205,8 @@
+@@ -205,8 +205,8 @@ def evaluate(read_fname,
prev_read_id = read_id
@@ -2352,15 +2376,17 @@ Description: 2to3
os.chdir("..")
---- hisat2.orig/evaluation/tests/one_snp_test/simulate_one_snp_reads.py
-+++ hisat2/evaluation/tests/one_snp_test/simulate_one_snp_reads.py
+Index: hisat2.git/evaluation/tests/one_snp_test/simulate_one_snp_reads.py
+===================================================================
+--- hisat2.git.orig/evaluation/tests/one_snp_test/simulate_one_snp_reads.py
++++ hisat2.git/evaluation/tests/one_snp_test/simulate_one_snp_reads.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
#
# Copyright 2015, Daehwan Kim <infphilo@gmail.com>
#
-@@ -125,8 +125,8 @@
+@@ -125,8 +125,8 @@ def getSamAlignment(chr_seq, read_len, s
Zs = "%d|I|%s" % (left_read_len, snp_id)
if len(read_seq) != read_len:
@@ -2371,7 +2397,7 @@ Description: 2to3
assert False
ref_read_seq = chr_seq[pos:pos+read_len]
-@@ -257,8 +257,8 @@
+@@ -257,8 +257,8 @@ def samRepOk(genome_seq, read_seq, chr,
tMD += ("{}".format(match_len))
if tMD != MD:
@@ -2382,7 +2408,7 @@ Description: 2to3
assert False
-@@ -272,14 +272,14 @@
+@@ -272,14 +272,14 @@ def simulate_reads(genome_file, snp_file
genome_seq = read_genome(genome_file)
snps = read_snp(snp_file)
@@ -2401,7 +2427,7 @@ Description: 2to3
read_file = open(base_fname + "_snp_1.fa", "w")
ref_read_file = open(base_fname + "_ref_1.fa", "w")
-@@ -312,12 +312,12 @@
+@@ -312,12 +312,12 @@ def simulate_reads(genome_file, snp_file
# Zs2 = ("\tZs:Z:{}".format(Zs2))
read_id_str = "{}_{}_{}_{}".format(cur_read_id, chr, pos, cigar_str)
@@ -2419,15 +2445,17 @@ Description: 2to3
"""
if paired_end:
print >> read2_file, ">{}".format(cur_read_id)
---- hisat2.orig/evaluation/tests/repeat/generate_repeats.py
-+++ hisat2/evaluation/tests/repeat/generate_repeats.py
+Index: hisat2.git/evaluation/tests/repeat/generate_repeats.py
+===================================================================
+--- hisat2.git.orig/evaluation/tests/repeat/generate_repeats.py
++++ hisat2.git/evaluation/tests/repeat/generate_repeats.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
import sys
import struct
-@@ -25,7 +25,7 @@
+@@ -25,7 +25,7 @@ while True:
chr_sa.append(num)
if len(chr_sa) % 5000000 == 0:
@@ -2436,7 +2464,7 @@ Description: 2to3
f.close()
assert chr_sa[0] + 1 == len(chr_sa)
-@@ -41,7 +41,7 @@
+@@ -41,7 +41,7 @@ repeats = []
while i < len(chr_sa) - 1:
pos = chr_sa[i]
base_seq = chr_seq[pos:pos+seq_len]
@@ -2445,7 +2473,7 @@ Description: 2to3
pos2 = chr_sa[j]
cmp_seq = chr_seq[pos2:pos2+seq_len]
if base_seq != cmp_seq:
-@@ -53,13 +53,13 @@
+@@ -53,13 +53,13 @@ while i < len(chr_sa) - 1:
i = j
if i % 5000000 == 0:
@@ -2463,7 +2491,7 @@ Description: 2to3
if j in deleted:
continue
-@@ -84,17 +84,17 @@
+@@ -84,17 +84,17 @@ for i in xrange(len(repeats) - 1):
if num_close == 1 and num_close2 < 5:
found = True
@@ -2486,7 +2514,7 @@ Description: 2to3
nt = pos_seq2_rc[k]
if nt == 'A':
nt = 'T'
-@@ -108,11 +108,11 @@
+@@ -108,11 +108,11 @@ for i in xrange(len(repeats) - 1):
pos_seq2_rc[k] = nt
pos_seq2_rc = ''.join(pos_seq2_rc)
@@ -2503,7 +2531,7 @@ Description: 2to3
file1.close()
file2.close()
-@@ -121,7 +121,7 @@
+@@ -121,7 +121,7 @@ for i in xrange(len(repeats) - 1):
if found:
break
@@ -2512,7 +2540,7 @@ Description: 2to3
chr_seq = ""
for line in open("%s.fa" % chr_name):
-@@ -132,7 +132,7 @@
+@@ -132,7 +132,7 @@ for line in open("%s.fa" % chr_name):
N_ranges = []
prev_nt = None
@@ -2521,7 +2549,7 @@ Description: 2to3
nt = chr_seq[i]
if nt == 'N':
if prev_nt != 'N':
-@@ -158,7 +158,7 @@
+@@ -158,7 +158,7 @@ for N_start, N_end in N_ranges:
to_genome_list = [[y, x] for x, y in to_joined_list]
N_ranges_tmp = []
@@ -2530,7 +2558,7 @@ Description: 2to3
to_genome = to_genome_list[i]
if i == 0:
if to_genome[1] > 0:
-@@ -171,8 +171,8 @@
+@@ -171,8 +171,8 @@ assert N_ranges == N_ranges_tmp
file = open("%s_rep.info" % chr_name, "w")
def print_rep_info(rep_name, rep_pos, rep_len, pos_set, pos_seq):
@@ -2541,7 +2569,7 @@ Description: 2to3
output = ""
for j in range(i, i + 10):
if j >= len(pos_set):
-@@ -181,7 +181,7 @@
+@@ -181,7 +181,7 @@ def print_rep_info(rep_name, rep_pos, re
output += " "
def convert(pos):
@@ -2550,7 +2578,7 @@ Description: 2to3
if i + 1 == len(to_genome_list) or (pos >= to_genome_list[i][0] and pos < to_genome_list[i+1][0]):
return pos - to_genome_list[i][0] + to_genome_list[i][1]
-@@ -190,7 +190,7 @@
+@@ -190,7 +190,7 @@ def print_rep_info(rep_name, rep_pos, re
pos = convert(pos_set[j])
assert chr_seq[pos:pos+seq_len] == pos_seq
output += ("%s:%d:+" % (chr_name, pos))
@@ -2559,7 +2587,7 @@ Description: 2to3
print_rep_info("rep1", 0, seq_len, pos_set, pos_seq)
print_rep_info("rep2", seq_len, seq_len, pos_set2, pos_seq2)
file.close()
-@@ -198,15 +198,15 @@
+@@ -198,15 +198,15 @@ file.close()
chr_seq = chr_seq.replace(pos_seq, 'N' * seq_len)
chr_seq = chr_seq.replace(pos_seq2, 'N' * seq_len)
file = open("%s_mask.fa" % chr_name, "w")
@@ -2581,15 +2609,17 @@ Description: 2to3
+ print(rep_seq[i:i+60], file=file)
file.close()
---- hisat2.orig/evaluation/tests/repeat/test_repeat.py
-+++ hisat2/evaluation/tests/repeat/test_repeat.py
+Index: hisat2.git/evaluation/tests/repeat/test_repeat.py
+===================================================================
+--- hisat2.git.orig/evaluation/tests/repeat/test_repeat.py
++++ hisat2.git/evaluation/tests/repeat/test_repeat.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
import sys, os, subprocess, random
from argparse import ArgumentParser, FileType
-@@ -52,7 +52,7 @@
+@@ -52,7 +52,7 @@ def read_genome(genome_filename):
genome_file.close()
@@ -2598,7 +2628,7 @@ Description: 2to3
return chr_dic
-@@ -62,7 +62,7 @@
+@@ -62,7 +62,7 @@ def read_genome(genome_filename):
def generate_random_seq(seq_len):
assert seq_len > 0
random_seq = ""
@@ -2607,7 +2637,7 @@ Description: 2to3
random_seq += "ACGT"[random.randint(0, 3)]
return random_seq
-@@ -109,10 +109,10 @@
+@@ -109,10 +109,10 @@ def test_repeat(verbose):
]
for id, seq in seqs:
@@ -2622,33 +2652,41 @@ Description: 2to3
"""
---- hisat2.orig/hisat2-build
-+++ hisat2/hisat2-build
+Index: hisat2.git/hisat2-build
+===================================================================
+--- hisat2.git.orig/hisat2-build
++++ hisat2.git/hisat2-build
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
"""
Copyright 2015, Daehwan Kim <infphilo@gmail.com>
---- hisat2.orig/hisat2-build-new
-+++ hisat2/hisat2-build-new
+Index: hisat2.git/hisat2-build-new
+===================================================================
+--- hisat2.git.orig/hisat2-build-new
++++ hisat2.git/hisat2-build-new
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
"""
Copyright 2018, Chanhee Park <parkchanhee@gmail.com> and Daehwan Kim <infphilo@gmail.com>
---- hisat2.orig/hisat2-inspect
-+++ hisat2/hisat2-inspect
+Index: hisat2.git/hisat2-inspect
+===================================================================
+--- hisat2.git.orig/hisat2-inspect
++++ hisat2.git/hisat2-inspect
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
"""
Copyright 2015, Daehwan Kim <infphilo@gmail.com>
---- hisat2.orig/hisat2_extract_exons.py
-+++ hisat2/hisat2_extract_exons.py
-@@ -63,7 +63,7 @@
+Index: hisat2.git/hisat2_extract_exons.py
+===================================================================
+--- hisat2.git.orig/hisat2_extract_exons.py
++++ hisat2.git/hisat2_extract_exons.py
+@@ -63,7 +63,7 @@ def extract_exons(gtf_file, verbose = Fa
trans[transcript_id][2].append([left, right])
# Sort exons and merge where separating introns are <=5 bps
@@ -2657,7 +2695,7 @@ Description: 2to3
exons.sort()
tmp_exons = [exons[0]]
for i in range(1, len(exons)):
-@@ -75,7 +75,7 @@
+@@ -75,7 +75,7 @@ def extract_exons(gtf_file, verbose = Fa
# Calculate and print the unique junctions
tmp_exons = set()
@@ -2666,9 +2704,11 @@ Description: 2to3
for i in range(len(texons)):
tmp_exons.add((chrom, texons[i][0], texons[i][1], strand))
tmp_exons = sorted(tmp_exons)
---- hisat2.orig/hisat2_extract_snps_haplotypes_UCSC.py
-+++ hisat2/hisat2_extract_snps_haplotypes_UCSC.py
-@@ -155,7 +155,7 @@
+Index: hisat2.git/hisat2_extract_snps_haplotypes_UCSC.py
+===================================================================
+--- hisat2.git.orig/hisat2_extract_snps_haplotypes_UCSC.py
++++ hisat2.git/hisat2_extract_snps_haplotypes_UCSC.py
+@@ -155,7 +155,7 @@ def generate_haplotypes(snp_file,
vars_count[id] = 0
vars_count[id] += 1
vars_duplicate = set()
@@ -2677,7 +2717,7 @@ Description: 2to3
if count <= 1:
continue
vars_duplicate.add(id)
-@@ -302,7 +302,7 @@
+@@ -302,7 +302,7 @@ def generate_haplotypes(snp_file,
h_end += (int(h2_data) - 1)
assert h_begin <= h_end
h_new_begin = h_begin
@@ -2686,9 +2726,11 @@ Description: 2to3
hc = haplotypes[h_j].split('#')
_, hc_begin, hc_type, hc_data, _ = vars[int(hc[-1])]
hc_begin = int(hc_begin)
---- hisat2.orig/hisat2_extract_snps_haplotypes_VCF.py
-+++ hisat2/hisat2_extract_snps_haplotypes_VCF.py
-@@ -409,7 +409,7 @@
+Index: hisat2.git/hisat2_extract_snps_haplotypes_VCF.py
+===================================================================
+--- hisat2.git.orig/hisat2_extract_snps_haplotypes_VCF.py
++++ hisat2.git/hisat2_extract_snps_haplotypes_VCF.py
+@@ -409,7 +409,7 @@ def generate_haplotypes(snp_file,
h_end += (int(h2_data) - 1)
assert h_begin <= h_end
h_new_begin = h_begin
@@ -2697,7 +2739,7 @@ Description: 2to3
hc = haplotypes[h_j].split('#')
_, hc_begin, hc_type, hc_data, _ = vars[int(hc[-1])]
hc_begin = int(hc_begin)
-@@ -776,7 +776,7 @@
+@@ -776,7 +776,7 @@ def main(genome_file,
vars = []
else:
@@ -2706,9 +2748,11 @@ Description: 2to3
chr_seq = chr_dic[chr]
chr_genotype_vars = genotype_var_list[chr]
curr_right = -1
---- hisat2.orig/hisat2_extract_splice_sites.py
-+++ hisat2/hisat2_extract_splice_sites.py
-@@ -64,7 +64,7 @@
+Index: hisat2.git/hisat2_extract_splice_sites.py
+===================================================================
+--- hisat2.git.orig/hisat2_extract_splice_sites.py
++++ hisat2.git/hisat2_extract_splice_sites.py
+@@ -64,7 +64,7 @@ def extract_splice_sites(gtf_file, verbo
trans[transcript_id][2].append([left, right])
# Sort exons and merge where separating introns are <=5 bps
@@ -2717,7 +2761,7 @@ Description: 2to3
exons.sort()
tmp_exons = [exons[0]]
for i in range(1, len(exons)):
-@@ -76,7 +76,7 @@
+@@ -76,7 +76,7 @@ def extract_splice_sites(gtf_file, verbo
# Calculate and print the unique junctions
junctions = set()
@@ -2726,7 +2770,7 @@ Description: 2to3
for i in range(1, len(exons)):
junctions.add((chrom, exons[i-1][1], exons[i][0], strand))
junctions = sorted(junctions)
-@@ -88,7 +88,7 @@
+@@ -88,7 +88,7 @@ def extract_splice_sites(gtf_file, verbo
if verbose:
exon_lengths, intron_lengths, trans_lengths = \
Counter(), Counter(), Counter()
@@ -2735,7 +2779,7 @@ Description: 2to3
tran_len = 0
for i, exon in enumerate(exons):
exon_len = exon[1]-exon[0]+1
-@@ -100,7 +100,7 @@
+@@ -100,7 +100,7 @@ def extract_splice_sites(gtf_file, verbo
trans_lengths[tran_len] += 1
print('genes: {}, genes with multiple isoforms: {}'.format(
@@ -2744,9 +2788,11 @@ Description: 2to3
file=stderr)
print('transcripts: {}, transcript avg. length: {:.0f}'.format(
len(trans), sum(trans_lengths.elements())//len(trans)),
---- hisat2.orig/hisat2_simulate_reads.py
-+++ hisat2/hisat2_simulate_reads.py
-@@ -152,7 +152,7 @@
+Index: hisat2.git/hisat2_simulate_reads.py
+===================================================================
+--- hisat2.git.orig/hisat2_simulate_reads.py
++++ hisat2.git/hisat2_simulate_reads.py
+@@ -152,7 +152,7 @@ def read_transcript(genome_seq, gtf_file
transcripts[transcript_id][2].append([left, right])
# Sort exons and merge where separating introns are <=5 bps
@@ -2755,7 +2801,7 @@ Description: 2to3
exons.sort()
tmp_exons = [exons[0]]
for i in range(1, len(exons)):
-@@ -163,7 +163,7 @@
+@@ -163,7 +163,7 @@ def read_transcript(genome_seq, gtf_file
transcripts[tran] = [chr, strand, tmp_exons]
tmp_transcripts = {}
@@ -2764,7 +2810,7 @@ Description: 2to3
exon_lens = [e[1] - e[0] + 1 for e in exons]
transcript_len = sum(exon_lens)
if transcript_len >= frag_len:
-@@ -199,7 +199,7 @@
+@@ -199,7 +199,7 @@ def read_snp(snp_file):
"""
def sanity_check_input(genome_seq, genes, transcripts, snps, frag_len):
num_canon_ss, num_ss = 0, 0
@@ -2773,7 +2819,7 @@ Description: 2to3
assert transcript_len >= frag_len
if len(exons) <= 1:
continue
-@@ -274,7 +274,7 @@
+@@ -274,7 +274,7 @@ def generate_rna_expr_profile(expr_profi
"""
def generate_dna_expr_profile(genome_seq):
expr_profile = []
@@ -2782,15 +2828,17 @@ Description: 2to3
expr_profile.append(len(chr_seq))
expr_sum = float(sum(expr_profile))
expr_profile = [expr_profile[i] / expr_sum for i in range(len(expr_profile))]
---- hisat2.orig/hisat2lib/pymodule/ht2example.py
-+++ hisat2/hisat2lib/pymodule/ht2example.py
+Index: hisat2.git/hisat2lib/pymodule/ht2example.py
+===================================================================
+--- hisat2.git.orig/hisat2lib/pymodule/ht2example.py
++++ hisat2.git/hisat2lib/pymodule/ht2example.py
@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/python3
#
# Copyright 2018, Chanhee Park <parkchanhee@gmail.com> and Daehwan Kim <infphilo@gmail.com>
-@@ -30,7 +30,7 @@
+@@ -30,7 +30,7 @@ ht2_index = '../../evaluation/indexes/HI
# Get default options
ht2_options = ht2py.get_options()
@@ -2799,7 +2847,7 @@ Description: 2to3
ht2_options['gVerbose'] = 1
ht2_options['startVerbose'] = 1
# or
-@@ -38,7 +38,7 @@
+@@ -38,7 +38,7 @@ ht2_options = {}
handle = ht2py.init(ht2_index, ht2_options)
@@ -2808,7 +2856,7 @@ Description: 2to3
#print ht2py.index_getrefnamebyid(handle, 0, 1, 3, 5, 7, 9)
# outofindex
-@@ -62,7 +62,7 @@
+@@ -62,7 +62,7 @@ for pos in positions:
if direction == 1:
chr_dir = '-'
@@ -2817,23 +2865,27 @@ Description: 2to3
# close handle
ht2py.close(handle)
---- hisat2.orig/hisat2lib/pymodule/setup.py
-+++ hisat2/hisat2lib/pymodule/setup.py
+Index: hisat2.git/hisat2lib/pymodule/setup.py
+===================================================================
+--- hisat2.git.orig/hisat2lib/pymodule/setup.py
++++ hisat2.git/hisat2lib/pymodule/setup.py
@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/python3
#
# Copyright 2018, Chanhee Park <parkchanhee@gmail.com> and Daehwan Kim <infphilo@gmail.com>
---- hisat2.orig/scripts/sa.py
-+++ hisat2/scripts/sa.py
+Index: hisat2.git/scripts/sa.py
+===================================================================
+--- hisat2.git.orig/scripts/sa.py
++++ hisat2.git/scripts/sa.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
"""
sa.py
-@@ -15,7 +15,7 @@
+@@ -15,7 +15,7 @@ import struct
def loadBowtieSa(fh):
""" Load a .sa file from handle into an array of ints """
nsa = struct.unpack('I', fh.read(4))[0]
@@ -2842,7 +2894,7 @@ Description: 2to3
def loadBowtieSaFilename(fn):
""" Load a .sa file from filename into an array of ints """
-@@ -58,7 +58,7 @@
+@@ -58,7 +58,7 @@ if __name__ == "__main__":
# Suffix array is in sas; note that $ is considered greater than all
# other characters
if ref is not None:
@@ -2851,22 +2903,24 @@ Description: 2to3
sa1, sa2 = sas[i-1], sas[i]
assert sa1 != sa2
# Sanity check that suffixes are really in order
-@@ -76,4 +76,4 @@
+@@ -76,4 +76,4 @@ if __name__ == "__main__":
assert sas[-1] == len(ref)
go()
-
\ No newline at end of file
+
---- hisat2.orig/scripts/validate_repeat.py
-+++ hisat2/scripts/validate_repeat.py
+Index: hisat2.git/scripts/validate_repeat.py
+===================================================================
+--- hisat2.git.orig/scripts/validate_repeat.py
++++ hisat2.git/scripts/validate_repeat.py
@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/python3
import sys, subprocess
import re
from argparse import ArgumentParser, FileType
-@@ -203,7 +203,7 @@
+@@ -203,7 +203,7 @@ def main(genome_file, rpt_name):
seq = reverse_complement(seq)
if seq != repeat_sequence:
diff --git a/debian/patches/reproducibile b/debian/patches/reproducibile
index 6ab212d..f05713b 100644
--- a/debian/patches/reproducibile
+++ b/debian/patches/reproducibile
@@ -1,9 +1,11 @@
From: Michael R. Crusoe <crusoe@debian.org>
Subject: Make build reproducible
Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/250
---- hisat2.orig/Makefile
-+++ hisat2/Makefile
-@@ -26,7 +26,7 @@
+Index: hisat2.git/Makefile
+===================================================================
+--- hisat2.git.orig/Makefile
++++ hisat2.git/Makefile
+@@ -26,7 +26,7 @@ GCC_SUFFIX =
CC = $(GCC_PREFIX)/gcc$(GCC_SUFFIX)
CPP = $(GCC_PREFIX)/g++$(GCC_SUFFIX)
CXX = $(CPP)
@@ -12,7 +14,7 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/250
BOWTIE_MM = 1
BOWTIE_SHARED_MEM = 0
-@@ -149,7 +149,7 @@
+@@ -149,7 +149,7 @@ HISAT2_CPPS_MAIN = $(SEARCH_CPPS) hisat2
HISAT2_BUILD_CPPS_MAIN = $(BUILD_CPPS) hisat2_build_main.cpp
HISAT2_REPEAT_CPPS_MAIN = $(REPEAT_CPPS) $(BUILD_CPPS) hisat2_repeat_main.cpp
@@ -21,7 +23,7 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/250
VERSION = $(shell cat VERSION)
# Convert BITS=?? to a -m flag
-@@ -234,11 +234,11 @@
+@@ -234,11 +234,11 @@ HT2LIB_PKG_SRC = \
$(HT2LIB_DIR)/pymodule/ht2example.py
@@ -38,7 +40,7 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/250
example/reference/22_20-21M.fa \
example/reference/22_20-21M.snp \
$(PTHREAD_PKG) \
-@@ -261,10 +261,10 @@
+@@ -261,10 +261,10 @@ endif
# the Windows FIND tool instead.
FIND=$(shell which find)
@@ -53,7 +55,7 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/250
$(HT2LIB_PKG_SRC) \
Makefile \
CMakeLists.txt \
-@@ -286,15 +286,20 @@
+@@ -286,15 +286,20 @@ repeat: hisat2-repeat
repeat-debug: hisat2-repeat-debug
diff --git a/debian/patches/riscv5 b/debian/patches/riscv5
index 11c9a3e..2d19f71 100644
--- a/debian/patches/riscv5
+++ b/debian/patches/riscv5
@@ -1,8 +1,10 @@
From: Michael R. Crusoe <crusoe@debian.org>
Subject: Fix compilation on riscv5
---- hisat2.orig/alphabet.cpp
-+++ hisat2/alphabet.cpp
-@@ -400,7 +400,7 @@
+Index: hisat2.git/alphabet.cpp
+===================================================================
+--- hisat2.git.orig/alphabet.cpp
++++ hisat2.git/alphabet.cpp
+@@ -400,7 +400,7 @@ int dnacomp[5] = {
const char *iupacs = "!ACMGRSVTWYHKDBN!acmgrsvtwyhkdbn";
@@ -11,9 +13,11 @@ Subject: Fix compilation on riscv5
-1,
'A', // 0001
'C', // 0010
---- hisat2.orig/alphabet.h
-+++ hisat2/alphabet.h
-@@ -65,7 +65,7 @@
+Index: hisat2.git/alphabet.h
+===================================================================
+--- hisat2.git.orig/alphabet.h
++++ hisat2.git/alphabet.h
+@@ -65,7 +65,7 @@ extern uint8_t dinuc2color[5][5];
/// corresponding 2-bit nucleotide
extern uint8_t nuccol2nuc[5][5];
/// Convert a 4-bit mask into an IUPAC code
@@ -22,9 +26,11 @@ Subject: Fix compilation on riscv5
/// Convert an ascii color to an ascii dna char
extern char col2dna[];
---- hisat2.orig/sstring.h
-+++ hisat2/sstring.h
-@@ -3326,7 +3326,7 @@
+Index: hisat2.git/sstring.h
+===================================================================
+--- hisat2.git.orig/sstring.h
++++ hisat2.git/sstring.h
+@@ -3326,7 +3326,7 @@ public:
*/
char toChar(size_t idx) const {
assert_range((int)this->cs_[idx], 0, 15);
diff --git a/debian/patches/simde b/debian/patches/simde
index 61f774b..166b818 100644
--- a/debian/patches/simde
+++ b/debian/patches/simde
@@ -1,8 +1,10 @@
From: Michael R. Crusoe <crusoe@debian.org>
Subject: use the simde header library for greater compatibility
Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/251
---- hisat2.orig/aligner_sw.h
-+++ hisat2/aligner_sw.h
+Index: hisat2.git/aligner_sw.h
+===================================================================
+--- hisat2.git.orig/aligner_sw.h
++++ hisat2.git/aligner_sw.h
@@ -66,11 +66,12 @@
#define INLINE_CUPS
@@ -17,8 +19,10 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/251
#include "aligner_sw_common.h"
#include "aligner_sw_nuc.h"
#include "ds.h"
---- hisat2.orig/sse_util.h
-+++ hisat2/sse_util.h
+Index: hisat2.git/sse_util.h
+===================================================================
+--- hisat2.git.orig/sse_util.h
++++ hisat2.git/sse_util.h
@@ -20,11 +20,13 @@
#ifndef SSE_UTIL_H_
#define SSE_UTIL_H_
@@ -34,9 +38,11 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/251
class EList_m128i {
public:
---- hisat2.orig/Makefile
-+++ hisat2/Makefile
-@@ -55,7 +55,7 @@
+Index: hisat2.git/Makefile
+===================================================================
+--- hisat2.git.orig/Makefile
++++ hisat2.git/Makefile
+@@ -55,7 +55,7 @@ ifneq (,$(findstring Darwin,$(shell unam
MACOS = 1
endif
@@ -45,8 +51,10 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/251
INC += -I. -I third_party
MM_DEF =
---- hisat2.orig/processor_support.h
-+++ hisat2/processor_support.h
+Index: hisat2.git/processor_support.h
+===================================================================
+--- hisat2.git.orig/processor_support.h
++++ hisat2.git/processor_support.h
@@ -12,7 +12,7 @@
#if defined(__INTEL_COMPILER)
@@ -56,7 +64,7 @@ Forwarded: https://github.com/DaehwanKimLab/hisat2/pull/251
# define USING_GCC_COMPILER
# include <cpuid.h>
#elif defined(_MSC_VER)
-@@ -52,8 +52,7 @@
+@@ -52,8 +52,7 @@ public:
#elif defined(USING_GCC_COMPILER)
__get_cpuid(0x1, ®s.EAX, ®s.EBX, ®s.ECX, ®s.EDX);
#else
diff --git a/docs/_data/download-binary.yml b/docs/_data/download-binary.yml
index a779fa2..eaf5caf 100644
--- a/docs/_data/download-binary.yml
+++ b/docs/_data/download-binary.yml
@@ -5,8 +5,8 @@ release:
name: HISAT2
artifacts:
Source: https://cloud.biohpc.swmed.edu/index.php/s/fE9QCsX3NH4QwBi/download
- OSX_x86_64: https://cloud.biohpc.swmed.edu/index.php/s/an8KdGxxRdSRXjr/download
- Linux_x86_64: https://cloud.biohpc.swmed.edu/index.php/s/4pMgDq4oAF9QCfA/download
+ OSX_x86_64: https://cloud.biohpc.swmed.edu/index.php/s/zMgEtnF6LjnjFrr/download
+ Linux_x86_64: https://cloud.biohpc.swmed.edu/index.php/s/oTtGWbWjaxsQ2Ho/download
- version: 2.2.0
date: 2/6/2020
name: HISAT2
diff --git a/docs/_data/download-index.yml b/docs/_data/download-index.yml
index 08bb039..ad1f1dd 100644
--- a/docs/_data/download-index.yml
+++ b/docs/_data/download-index.yml
@@ -2,38 +2,38 @@
data:
GRCh38:
genome:
- url: https://cloud.biohpc.swmed.edu/index.php/s/grch38/download
+ url: https://genome-idx.s3.amazonaws.com/hisat/grch38_genome.tar.gz
genome_snp:
- url: https://cloud.biohpc.swmed.edu/index.php/s/grch38_snp/download
+ url: https://genome-idx.s3.amazonaws.com/hisat/grch38_snp.tar.gz
genome_tran:
- url: https://cloud.biohpc.swmed.edu/index.php/s/grch38_tran/download
+ url: https://genome-idx.s3.amazonaws.com/hisat/grch38_tran.tar.gz
genome_snp_tran:
- url: https://cloud.biohpc.swmed.edu/index.php/s/grch38_snp_tran/download
+ url: https://genome-idx.s3.amazonaws.com/hisat/grch38_snptran.tar.gz
genome_rep(above 2.2.0):
- url: https://cloud.biohpc.swmed.edu/index.php/s/grch38_rep/download
+ url: https://genome-idx.s3.amazonaws.com/hisat/grch38_rep.tar.gz
genome_snp_rep(above 2.2.0):
- url: https://cloud.biohpc.swmed.edu/index.php/s/grch38_snp_rep/download
+ url: https://genome-idx.s3.amazonaws.com/hisat/grch38_snprep.tar.gz
UCSC hg38:
genome:
- url: https://cloud.biohpc.swmed.edu/index.php/s/hg38/download
+ url: https://genome-idx.s3.amazonaws.com/hisat/hg38_genome.tar.gz
genome_tran:
- url: https://cloud.biohpc.swmed.edu/index.php/s/hg38_tran/download
+ url: https://genome-idx.s3.amazonaws.com/hisat/hg38_tran.tar.gz
GRCh37:
genome:
- url: https://cloud.biohpc.swmed.edu/index.php/s/grch37/download
+ url: https://genome-idx.s3.amazonaws.com/hisat/grch37_genome.tar.gz
genome_snp:
- url: https://cloud.biohpc.swmed.edu/index.php/s/grch37_snp/download
+ url: https://genome-idx.s3.amazonaws.com/hisat/grch37_snp.tar.gz
genome_tran:
- url: https://cloud.biohpc.swmed.edu/index.php/s/grch37_tran/download
+ url: https://genome-idx.s3.amazonaws.com/hisat/grch37_tran.tar.gz
genome_snp_tran:
- url: https://cloud.biohpc.swmed.edu/index.php/s/grch37_snp_tran/download
+ url: https://genome-idx.s3.amazonaws.com/hisat/grch37_snptran.tar.gz
UCSC hg19:
genome:
- url: https://cloud.biohpc.swmed.edu/index.php/s/hg19/download
+ url: https://genome-idx.s3.amazonaws.com/hisat/hg19_genome.tar.gz
- organism: M. musculus
data:
GRCm38:
- genome:
+ genome:
url: https://cloud.biohpc.swmed.edu/index.php/s/grcm38/download
genome_snp:
url: https://cloud.biohpc.swmed.edu/index.php/s/grcm38_snp/download
@@ -43,39 +43,39 @@
url: https://cloud.biohpc.swmed.edu/index.php/s/grcm38_snp_tran/download
UCSC mm10:
genome:
- url: https://cloud.biohpc.swmed.edu/index.php/s/mm10/download
+ url: https://genome-idx.s3.amazonaws.com/hisat/mm10_genome.tar.gz
- organism: R. norvegicus
data:
UCSC rn6:
genome:
- url: https://cloud.biohpc.swmed.edu/index.php/s/rn6/download
+ url: https://genome-idx.s3.amazonaws.com/hisat/rn6_genome.tar.gz
- organism: D. melanogaster
data:
BDGP6:
genome:
- url: ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/bdgp6.tar.gz
+ url: https://genome-idx.s3.amazonaws.com/hisat/bdgp6.tar.gz
genome_tran:
- url: ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/bdgp6_tran.tar.gz
+ url: https://genome-idx.s3.amazonaws.com/hisat/bdgp6_tran.tar.gz
UCSC dm6:
genome:
- url: ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/dm6.tar.gz
+ url: https://genome-idx.s3.amazonaws.com/hisat/dm6.tar.gz
- organism: C. elegans
data:
WBcel235:
genome:
- url: ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/wbcel235.tar.gz
+ url: https://genome-idx.s3.amazonaws.com/hisat/wbcel235.tar.gz
genome_tran:
- url: ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/wbcel235_tran.tar.gz
+ url: https://genome-idx.s3.amazonaws.com/hisat/wbcel235_tran.tar.gz
UCSC ce10:
genome:
- url: ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/ce10.tar.gz
+ url: https://cloud.biohpc.swmed.edu/index.php/s/bbynxoY2TPpRNQb/download
- organism: S. cerevisiae
data:
R64-1-1:
genome:
- url: ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/r64.tar.gz
- genome_tran:
- url: ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/r64_tran.tar.gz
+ url: https://cloud.biohpc.swmed.edu/index.php/s/JRSoKHD5cHfpCFE/download
+ genome_tran:
+ url: https://cloud.biohpc.swmed.edu/index.php/s/akeiMrGGtt5KoJY/download
UCSC sacCer3:
genome:
- url: ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/sc3.tar.gz
+ url: https://cloud.biohpc.swmed.edu/index.php/s/Gsq4goLW4TDAz4E/download
diff --git a/docs/_layouts/default.html b/docs/_layouts/default.html
index dbf5ab5..43975c4 100644
--- a/docs/_layouts/default.html
+++ b/docs/_layouts/default.html
@@ -113,6 +113,7 @@
<h2>Publications</h2>
<div style="font-size: 0.8em">
<ul>
+ <li>Zhang, Y., Park, C., Bennett, C., Thornton, M. and Kim, D. <a class="publication" href="https://doi.org/10.1101/gr.275193.120">Rapid and accurate alignment of nucleotide conversion sequencing reads with HISAT-3N.</a> <a class="publication" href="https://genome.cshlp.org"><i>Genome Research</i></a> <b>31</b>, 1290-1295 (2021).</li>
<li>Kim, D., Paggi, J.M., Park, C. <i>et al.</i> <a class="publication" href="https://doi.org/10.1038/s41587-019-0201-4">Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype.</a> <a class="publication" href="https://www.nature.com/nbt/"><i>Nat Biotechnol</i></a> <b>37</b>, 907–915 (2019).</li>
<li>Kim D, Langmead B and Salzberg SL. <a class="publication" href="https://doi.org/10.1038/nmeth.3317">HISAT: a fast spliced aligner with low memory requirements.</a> <a class="publication" href="https://www.nature.com/nmeth/"><i>Nature Methods</i></a> 2015</li>
<li>Pertea M, Kim D, Pertea G, Leek JT and Salzberg SL. <a class="publication" href="https://doi.org/10.1038/nprot.2016.095">Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown.</a> <a class="publication" href="https://www.nature.com/nprot/"><i>Nature Protocols</i></a> 2016</li>
diff --git a/docs/_pages/download.md b/docs/_pages/download.md
index 06b03a7..34dbba7 100644
--- a/docs/_pages/download.md
+++ b/docs/_pages/download.md
@@ -2,7 +2,7 @@
layout: page
title: Download
permalink: /download/
-order: 4
+order: 5
share: false
---
@@ -13,14 +13,21 @@ Please cite:
{:toc}
## Index
+HISAT2 indexes are hosted on AWS (Amazon Web Services), thanks to the AWS Public Datasets program. Click this [link](https://registry.opendata.aws/jhu-indexes/) for more details.
{% for item in site.data.download-index %}
### {{ item.organism }}
{% for data in item.data %}
<li>{{ data[0] }}</li>
-<table style="border-collapse: collapse; border: none;">
+<table class="table-noborder">
{% for genome in data[1] %}
-<tr style="border: none;"><td style="border: none;">{{ genome[0] }}</td><td style="border: none;"><a href="{{ genome[1].url }}">{{ genome[1].url }}</a></td></tr>
+<tr><td>{{ genome[0] }}</td>
+ <td>
+ {% for url in genome[1] %}
+ <a href="{{ url[1] }}">{{ url[1] }}</a><br/>
+ {% endfor %}
+ </td>
+</tr>
{% endfor %}
</table>
{% endfor %}
diff --git a/docs/_pages/hisat-3n.md b/docs/_pages/hisat-3n.md
new file mode 100644
index 0000000..692bdf1
--- /dev/null
+++ b/docs/_pages/hisat-3n.md
@@ -0,0 +1,261 @@
+---
+layout: page
+title: HISAT-3N
+permalink: /hisat-3n/
+order: 4
+share: false
+---
+
+HISAT-3N
+============
+
+Overview
+-----------------
+**HISAT-3N** (hierarchical indexing for spliced alignment of transcripts - 3 nucleotides)
+is designed for nucleotide conversion sequencing technologies and implemented based on HISAT2.
+There are two strategies for HISAT-3N to align nucleotide conversion sequencing reads: *standard mode* and *repeat mode*.
+The standard mode align reads with standard-3N index only, so it is fast and require small memory (~9GB for human genome alignment).
+The repeat mode align reads with both standard-3N index and repeat-3N index, then output 1,000 alignment result (the output number can be adjusted by `--repeat-limit`).
+The repeat mode can align nucleotide conversion reads more accurately,
+and it is only 10% slower than the standard mode with tiny more memory (repeat mode use about ~10.5GB) usage than the standard mode.
+
+HISAT-3N is developed based on [HISAT2], which is particularly optimized for RNA sequencing technology.
+HISAT-3N supports both strand-specific and non-strand reads.
+HISAT-3N can be used for any base-converted sequencing reads include [BS-seq], [SLAM-seq], [TAB-seq], [oxBS-seq], [TAPS], [scBS-seq], and [scSLAM-seq].
+
+[HISAT2]:https://github.com/DaehwanKimLab/hisat2
+[BS-seq]: https://en.wikipedia.org/wiki/Bisulfite_sequencing
+[SLAM-seq]: https://www.nature.com/articles/nmeth.4435
+[scBS-seq]: https://www.nature.com/articles/nmeth.3035
+[scSLAM-seq]: https://www.nature.com/articles/s41586-019-1369-y
+[TAPS]: https://www.nature.com/articles/s41587-019-0041-2
+[TAB-seq]: https://doi.org/10.1016/j.cell.2012.04.027
+[oxBS-seq]: https://science.sciencemag.org/content/336/6083/934
+
+
+Getting started
+============
+HISAT-3N requires a 64-bit computer running either Linux or Mac OS X and at least 16 GB of RAM.
+
+A few notes:
+
+1. The repeat 3N index building process requires 256 GB of RAM.
+2. The standard 3N index building requires no more than 16 GB of RAM.
+3. The alignment process with either standard or repeat index requires no more than 16 GB of RAM.
+4. [SAMtools] is required to sort SAM file for `hisat-3n-table`.
+
+Install
+------------
+
+ git clone https://github.com/DaehwanKimLab/hisat2.git hisat-3n
+ cd hisat-3n
+ git checkout -b hisat-3n origin/hisat-3n
+ make
+
+
+Make sure that you are in the `hisat-3n` branch
+
+
+Build a HISAT-3N index with `hisat-3n-build`
+-----------
+`hisat-3n-build` builds a 3N-index, which contains two hisat2 indexes, from a set of DNA sequences. For standard 3N-index,
+each index contains 16 files with suffix `.3n.*.*.ht2`.
+For repeat 3N-index, there are 16 more files in addition to the standard 3N-index, and they have the suffix
+`.3n.*.rep.*.ht2`.
+These files constitute the hisat-3n index and no other file is needed to alignment reads to the reference.
+
+* `--base-change <chr1,chr2>` argument is required for `hisat-3n-build` and `hisat-3n`.
+ Provide which base is converted in the sequencing process to another base. Please enter
+ 2 letters separated by ',' for this argument. The first letter(chr1) should be the converted base, the second letter(chr2) should be
+ the converted to base. For example, during slam-seq, some 'T' is converted to 'C',
+ please enter `--base-change T,C`. During bisulfite-seq, some 'C' is converted to 'T', please enter `--base-change C,T`.
+* Different conversion types may build the same hisat-3n index. Please check the table below for more detail.
+ Once you build the hisat-3n index with C to T conversion (for example BS-seq).
+ You can align the T to C conversion reads (for example SLAM-seq reads) with the same index.
+
+
+ | Conversion Types | HISAT-3N index suffix |
+ |:----------------------------------:|:-----------------------------:|
+ |C -> T<br>T -> C<br>A -> G<br>G -> A|.3n.CT.\*.ht2 <br>.3n.GA.\*.ht2|
+ |A -> C<br>C -> A<br>G -> T<br>T -> G|.3n.AC.\*.ht2 <br>.3n.TG.\*.ht2|
+ |A -> T<br>T -> A |.3n.AT.\*.ht2 <br>.3n.TA.\*.ht2|
+ |C -> G<br>G -> C |.3n.CG.\*.ht2 <br>.3n.GC.\*.ht2|
+
+#### Examples:
+ # Build the standard HISAT-3N index (with C to T conversion):
+ hisat-3n-build --base-change C,T genome.fa genome
+
+ # Build the repeat HISAT-3N index (with T to C conversion, require 256 GB memory for human genome index):
+ hisat-3n-build --base-change T,C --repeat-index genome.fa genome
+
+It is optional to make the graph index and add SNP or spice site information to the index, to increase the alignment accuracy.
+The graph index building may require more memory than the linear index building.
+For more detail, please check the [HISAT2 manual].
+
+[HISAT2 manual]:https://daehwankimlab.github.io/hisat2/manual/
+
+#### Examples:
+ # Build the standard HISAT-3N index integrated index with SNP information
+ hisat-3n-build --base-change C,T --snp genome.snp genome.fa genome
+
+ # Build the standard HISAT-3N integrated index with splice site information
+ hisat-3n-build --base-change C,T --ss genome.ss --exon genome.exon genome.fa genome
+
+ # Build the repeat HISAT-3N index integrated index with SNP information
+ hisat-3n-build --base-change C,T --repeat-index --snp genome.snp genome.fa genome
+
+ # Build the repeat HISAT-3N integrated index with splice site information
+ hisat-3n-build --base-change C,T --repeat-index --ss genome.ss --exon genome.exon genome.fa genome
+
+Alignment with `hisat-3n`
+------------
+After we build the HISAT-3N index, you are ready to use `hisat-3n` for alignment.
+HISAT-3N uses the HISAT2 argument but has some extra arguments. Please check [HISAT2 manual] for more detail.
+
+* `--base-change <chr1,chr2>`
+ Provide which base is converted in the sequencing process to another base. Please enter
+ 2 letters separated by ',' for this argument. The first letter(chr1) should be the converted base, the second letter(chr2) should be
+ the converted to base. For example, during slam-seq, some 'T' is converted to 'C',
+ please enter `--base-change T,C`. During bisulfite-seq, some 'C' is converted to 'T', please enter `--base-change C,T`.
+
+* `-x <hisat-3n-idx>`
+ The index for HISAT-3N. The basename is the name of the index files up to but not including the suffix `.3n.*.*.ht2` / etc.
+ For example, you build your index with basename 'genome' by HISAT-3N-build, please enter `-x genome`.
+
+* `--directional-mapping`
+ Make directional mapping. Please use this option only if your sequencing reads are generated from a strand-specific library.
+ The directional mapping mode is about 2x faster than the default (non-directional) mapping mode.
+
+* `--repeat-limit <int>`
+ You can set up the number of alignment will be checked for each repeat alignment. You may increase the number to let hisat-3n
+ output more, if a read has multiple mapping. We suggest the repeat limit number for paired-end reads alignment is no more
+ than 1,000,000. default: 1000.
+
+* `--unique-only`
+ Only output uniquely aligned reads. This option may conflict with the hisat2 option `--un`, `--al`,
+ `--un-conc`, and `--al-conc`. We will fix this problem in future version. Please do not use `--unique-only` if you want to use the hisat2 options above.
+
+#### Examples:
+ # Single-end slam-seq reads (T to C conversion, RNA) alignment with the standard 3N-index:
+ hisat-3n -x genome -f -U read.fa -S alignment_result.sam --base-change T,C --no-repeat-index
+
+ # Paired-end strand-specific bisulfite-seq read (C to T conversion) alignment with repeat 3N-index:
+ hisat-3n -x genome -f -1 read_1.fa -2 read_2.fa -S alignment_result.sam --base-change C,T --repeat --no-spliced-alignment --directional-mapping
+
+ # Single-end TAPS reads (have C to T conversion， RNA) alignment with the repeat 3N-index:
+ hisat-3n -x genome -q -U read.fq -S alignment_result.sam --base-change C,T --repeat
+
+
+#### Extra SAM tags generated by HISAT-3N:
+
+* `Yf:i:<N>`: Number of conversions are detected in the read.
+* `Zf:i:<N>`: Number of un-converted bases are detected in the read. `Yf` + `Zf` = total number of bases which can be converted in the read sequence.
+* `YZ:A:<A>`: The value `+` or `–` indicate the read is mapped to REF-3N (`+`) or REF-RC-3N (`-`).
+
+Generate a 3N-conversion-table with `hisat-3n-table`
+------------
+### Preparation
+
+To generate 3N-conversion-table, users need to sort the SAM file which generated by `hisat-3n`.
+[SAMtools] is required for this sorting process.
+
+Use `samtools sort` to convert the SAM file to a sorted SAM file.
+
+ samtools sort alignment_result.sam -o sorted_alignment_result.sam -O sam
+
+Generate 3N-conversion-table with `hisat-3n-table`:
+
+### Usage
+ hisat-3n-table [options]* --alignments <alignmentFile> --ref <refFile> --base-change <char1,char2>
+
+#### Main arguments
+* `--alignments <alignmentFile>`
+ SORTED SAM file. Please enter `-` for standard input.
+
+* `--ref <refFile>`
+ The reference genome file (FASTA format) for generating HISAT-3N index.
+
+* `--output-name <outputFile>`
+ Filename to write 3N-conversion-table (tsv format) to. By default, table is written to the “standard out” or “stdout” filehandle (i.e. the console).
+
+* `--base-change <char1,char2>`
+ The base-change rule. User should enter the exact same `--base-change` arguments in hisat-3n.
+ For example, please enter `--base-change C,T` for bisulfite sequencing reads.
+
+#### Input options
+* `-u/--unique-only`
+ Only count the unique aligned reads into 3N-conversion-table.
+
+* `-m/--multiple-only`
+ Only count the multiple aligned reads into 3N-conversion-table.
+
+* `-c/--CG-only`
+ Only count the CpG sites in reference genome. This option is designed for bisulfite sequencing reads.
+
+* `--added-chrname`
+ Please add this option if you use `--add-chrname` during `hisat-3n` alignment.
+ During `hisat-3n` alignment, the prefix "chr" is added in front of chromosome name and shows on SAM output, when user choose `--add-chrname`.
+ `hisat-3n-table` cannot find the chromosome name on reference because it has an additional "chr" prefix. This option is to help `hisat-3n-table`
+ find the matching chromosome name on reference file. The 3n-table provides the same chromosome name as SAM file.
+
+* `--removed-chrname`
+ Please add this option if you use `--remove-chrname` during `hisat-3n` alignment.
+ During `hisat-3n` alignment, the prefix "chr" is removed in front of chromosome name and shows on SAM output, when user choose `--remove-chrname`.
+ `hisat-3n-table` cannot find the chromosome name on reference because it has no "chr" prefix. This option is to help `hisat-3n-table`
+ find the matching chromosome name on reference file. The 3n-table provides the same chromosome name as SAM file.
+
+#### Other options:
+* `-p/--threads <int>`
+ Launch `int` parallel threads (default: 1) for table building.
+
+* `-h/--help`
+ Print usage information and quit.
+
+#### Examples:
+ # Generate the 3N-conversion-table for bisulfite sequencing data:
+ hisat-3n-table -p 16 --alignments sorted_alignment_result.sam --ref genome.fa --output-name output.tsv --base-change C,T
+
+ # Generate the 3N-conversion-table for TAPS data and only count base in CpG site and uniquely aligned:
+ hisat-3n-table -p 16 --alignments sorted_alignment_result.sam --ref genome.fa --output-name output.tsv --base-change C,T --CG-only --unique-only
+
+ # Generate the 3N-conversion-table for bisulfite sequencing data from sorted BAM file:
+ samtools view -h sorted_alignment_result.bam | hisat-3n-table --ref genome.fa --alignments - --output-name output.tsv --base-change C,T
+
+ # Generate the 3N-conversion-table for bisulfite sequencing data from unsorted BAM file:
+ samtools sort alignment_result.bam -O sam | hisat-3n-table --ref genome.fa --alignments - --output-name output.tsv --base-change C,T
+
+
+#### Note:
+There are 7 columns in the 3N-conversion-table:
+
+1. `ref`: the chromosome name.
+2. `pos`: 1-based position in ref.
+3. `strand`: '+' for forward strand. '-' for reverse strand.
+4. `convertedBaseQualities`: the qualities for converted base in read-level measurement. Length of this string is equal to
+ the number of converted Base in read-level measurement.
+5. `convertedBaseCount`: number of distinct read positions where converted base in read-level measurements were found.
+ this number should equal to the length of convertedBaseQualities.
+6. `unconvertedBaseQualities`: the qualities for unconverted base in read-level measurement. Length of this string is equal to
+ the number of unconverted Base in read-level measurement.
+7. `unconvertedBaseCount`: number of distinct read positions where unconverted base in read-level measurements were found.
+ this number should equal to the length of unconvertedBaseQualities.
+
+##### Sample 3N-conversion-table:
+ ref pos strand convertedBaseQualities convertedBaseCount unconvertedBaseQualities unconvertedBaseCount
+ 1 11874 + FFFFFB<BF<F 11 0
+ 1 11877 - FFFFFF< 7 0
+ 1 11878 + FFFBB//F/BB 11 0
+ 1 11879 + 0 FFFBB//FB/ 10
+ 1 11880 - F 1 FFFF/ 5
+
+[SAMtools]: http://samtools.sourceforge.net
+
+Publication
+============
+
+* HISAT-3N
+ Zhang, Y., Park, C., Bennett, C., Thornton, M. and Kim, D. [Rapid and accurate alignment of nucleotide conversion sequencing reads with HISAT-3N](https://doi.org/10.1101/gr.275193.120). _Genome Research_ **31(7)**: 1290-1295 (2021)
+
+
+* HIAST2
+ Kim, D., Paggi, J.M., Park, C. _et al._ [Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype](https://doi.org/10.1038/s41587-019-0201-4). _Nat Biotechnol_ **37**, 907–915 (2019)
diff --git a/docs/_pages/hisat2.md b/docs/_pages/hisat2.md
index 41b0458..2c7b47a 100644
--- a/docs/_pages/hisat2.md
+++ b/docs/_pages/hisat2.md
@@ -8,6 +8,18 @@ share: false
**HISAT2** is a fast and sensitive alignment program for mapping next-generation sequencing reads (both DNA and RNA) to a population of human genomes as well as to a single reference genome. Based on an extension of BWT for graphs ([Sirén et al. 2014](http://dl.acm.org/citation.cfm?id=2674828)), we designed and implemented a graph FM index (GFM), an original approach and its first implementation. In addition to using one global GFM index that represents a population of human genomes, **HISAT2** uses a large set of small GFM indexes that collectively cover the whole genome. These small indexes (called local indexes), combined with several alignment strategies, enable rapid and accurate alignment of sequencing reads. This new indexing scheme is called a Hierarchical Graph FM index (HGFM).
+### [The HISAT-3N paper](https://genome.cshlp.org/content/31/7/1290.abstract) published at *Genome Research*. 7/1/2021
+
+### HISAT-3N beta release 12/14/2020
+
+HISAT-3N is a software system for analyzing nucleotide conversion sequencing reads. See the [HISAT-3N] for more details.
+
+[HISAT-3N]: {{ site.baseurl }}{% link _pages/hisat-3n.md %}
+
+### Index files are moved to the AWS Public Dataset Program. 9/3/2020
+
+We have moved HISAT2 index files to the AWS Public Dataset Program. See the [link](https://registry.opendata.aws/jhu-indexes/) for more details.
+
### HISAT 2.2.1 release 7/24/2020
diff --git a/docs/_pages/howto.md b/docs/_pages/howto.md
index 3db1f75..3fc968d 100644
--- a/docs/_pages/howto.md
+++ b/docs/_pages/howto.md
@@ -2,7 +2,7 @@
layout: page
title: HowTo
permalink: /howto/
-order: 5
+order: 6
share: false
---
diff --git a/docs/_pages/links.md b/docs/_pages/links.md
index 861b971..db71e34 100644
--- a/docs/_pages/links.md
+++ b/docs/_pages/links.md
@@ -2,12 +2,13 @@
layout: page
title: Links
permalink: /links/
-order: 6
+order: 7
share: false
---
* KimLab - <https://kim-lab.org>
- * github - <https://github.com/DaehwanKimLab>
+ * lab github - <https://github.com/DaehwanKimLab>
+ * github for HISAT2 - <https://github.com/DaehwanKimLab/hisat2>
* hisat-genotype - <https://daehwankimlab.github.io/hisat-genotype>
* github for hisat-genotype - <https://github.com/DaehwanKimLab/hisat-genotype>
diff --git a/docs/_sass/_content.scss b/docs/_sass/_content.scss
index b028d5d..4f7503f 100644
--- a/docs/_sass/_content.scss
+++ b/docs/_sass/_content.scss
@@ -200,4 +200,11 @@
}
}
}
+
+ .table-noborder {
+ table, tbody, td, tr, th {
+ border-collapse: collapse;
+ border: none;
+ }
+ }
}
Debdiff
[The following lists of changes regard files as different if they have different names, permissions or owners.]
Files in second set of .debs but not in first
-rw-r--r-- root/root /usr/lib/debug/.build-id/06/c7cf0dfd9d395461b35679abae11f0fb2b6676.debug -rw-r--r-- root/root /usr/lib/debug/.build-id/07/2ffba8fa91651e0ba1e63904905e05d47caf94.debug -rw-r--r-- root/root /usr/lib/debug/.build-id/b0/d66980da44686dfd7eba3af0c3ccac53df1da1.debug -rw-r--r-- root/root /usr/lib/debug/.build-id/c3/bd18a57230bb22f143e83e90182da8bd3820e8.debug -rw-r--r-- root/root /usr/lib/debug/.build-id/de/51f7706198b6e90c4ce42cde81f3dd490a8675.debug -rw-r--r-- root/root /usr/lib/debug/.build-id/ef/95a9cc5cf42b259a7df952e21f386b17d09722.debug
Files in first set of .debs but not in second
-rw-r--r-- root/root /usr/lib/debug/.build-id/05/8c8f10abcd9c02b3beb6d167b7f3d7f3d864ba.debug -rw-r--r-- root/root /usr/lib/debug/.build-id/07/5bea55b34f8ce44b2cfef88ab3a196658b3b34.debug -rw-r--r-- root/root /usr/lib/debug/.build-id/28/28939e01f7cd252979c0d12e72a7de13adb570.debug -rw-r--r-- root/root /usr/lib/debug/.build-id/53/8de299f9d254bd9fadeeffed668f9e332135d7.debug -rw-r--r-- root/root /usr/lib/debug/.build-id/a8/db515488eee9d3e3ca7f4ce9f1de4dcfb90a0a.debug -rw-r--r-- root/root /usr/lib/debug/.build-id/e7/40551e01e92616ddf1df97574d22bc7c1bb507.debug
Control files of package hisat2: lines which differ (wdiff format)
Built-Using: simde (= 0.7.4~rc1-1) 0.7.2-6)
Control files of package hisat2-dbgsym: lines which differ (wdiff format)
Build-Ids: 058c8f10abcd9c02b3beb6d167b7f3d7f3d864ba 075bea55b34f8ce44b2cfef88ab3a196658b3b34 2828939e01f7cd252979c0d12e72a7de13adb570 538de299f9d254bd9fadeeffed668f9e332135d7 a8db515488eee9d3e3ca7f4ce9f1de4dcfb90a0a e740551e01e92616ddf1df97574d22bc7c1bb507 06c7cf0dfd9d395461b35679abae11f0fb2b6676 072ffba8fa91651e0ba1e63904905e05d47caf94 b0d66980da44686dfd7eba3af0c3ccac53df1da1 c3bd18a57230bb22f143e83e90182da8bd3820e8 de51f7706198b6e90c4ce42cde81f3dd490a8675 ef95a9cc5cf42b259a7df952e21f386b17d09722
Control files of package python3-hisat2: lines which differ (wdiff format)
Built-Using: simde (= 0.7.4~rc1-1) 0.7.2-6)