The Rhisat2 package provides an interface to the HISAT2 short-read aligner software (Kim, Langmead, and Salzberg 2015).
Use the BiocManager package to download and install the
package from Bioconductor as follows:
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("Rhisat2")If required, the latest development version of the package can also be installed from GitHub.
Once the package is installed, load it into your R session:
To build an index for the alignment, use the
hisat2_build function. You need to provide one or more
fasta file with reference sequences, as well as an output directory
where the index will be stored, and a “prefix” (that will determine the
name of the index files in the output directory). Any additional
arguments to the hisat2-build binary can also be supplied
to the function (see hisat2_build_usage() or https://ccb.jhu.edu/software/hisat2/manual.shtml#the-hisat2-build-indexer
for an overview of the available options).
The package contains three example reference sequences, corresponding
to short pieces of three different chromosomes. We show how to create an
index (named myindex) based on these sequences.
list.files(system.file("extdata/refs", package="Rhisat2"), pattern="\\.fa$")
#> [1] "chr1.fa" "chr2.fa" "chr3.fa"
refs <- list.files(system.file("extdata/refs", package="Rhisat2"),
full.names=TRUE, pattern="\\.fa$")
td <- tempdir()
hisat2_build(references=refs, outdir=td, prefix="myindex",
force=TRUE, strict=TRUE, execute=TRUE)By instead setting execute=FALSE in the command above,
hisat2_build() will return the index building shell command
for inspection, without executing it.
print(hisat2_build(references=refs, outdir=td, prefix="myindex",
force=TRUE, strict=TRUE, execute=FALSE))
#> [1] "'/tmp/RtmpiCXGfL/Rinst181e6e7307a9/Rhisat2/hisat2-build' '/tmp/RtmpiCXGfL/Rinst181e6e7307a9/Rhisat2/extdata/refs/chr1.fa','/tmp/RtmpiCXGfL/Rinst181e6e7307a9/Rhisat2/extdata/refs/chr2.fa','/tmp/RtmpiCXGfL/Rinst181e6e7307a9/Rhisat2/extdata/refs/chr3.fa' '/tmp/RtmpMUI5gC/myindex'"After creating the index, reads can be aligned using the
hisat2 wrapper function. Most commonly, the reads will be
provided in fastq files (one file for single-end reads, two
files for paired-end reads). The names of these files can be provided
directly to the hisat2 function, either as a vector (for
single-end reads) or as a list of length 2 (for paired-end reads, each
list element corresponds to one mate). You also need to provide the path
to the index generated by hisat2_build (the output
directory combined with the prefix), and the output file name where the
alignments should be written.
list.files(system.file("extdata/reads", package="Rhisat2"),
pattern="\\.fastq$")
#> [1] "reads1.fastq" "reads2.fastq"
reads <- list.files(system.file("extdata/reads", package="Rhisat2"),
pattern="\\.fastq$", full.names=TRUE)
hisat2(sequences=as.list(reads), index=file.path(td, "myindex"),
type="paired", outfile=file.path(td, "out.sam"),
force=TRUE, strict=TRUE, execute=TRUE)As for hisat2_build(), any additional arguments to the
hisat2 binary can be provided to the hisat2()
function (see hisat2_usage() or https://ccb.jhu.edu/software/hisat2/manual.shtml#running-hisat2
for an overview of the available options).
With HISAT2, you can provide a file with known splice sites at the
alignment step, which can help in finding the correct alignments across
known splice junctions. A text file in the required format can be
generated using the extract_splice_sites() function,
starting from an annotation file in gtf or gff3 format, or from a
GRanges or TxDb object.
spsfile <- tempfile()
gtf <- system.file("extdata/refs/genes.gtf", package="Rhisat2")
extract_splice_sites(features=gtf, outfile=spsfile)
#> Import genomic features from the file as a GRanges object ... OK
#> Prepare the 'metadata' data frame ... OK
#> Make the TxDb object ...
#> Warning in .makeTxDb_normarg_chrominfo(chrominfo): genome version information is not available for
#> this TxDb object
#> OK
hisat2(sequences=as.list(reads), index=file.path(td, "myindex"),
type="paired", outfile=file.path(td, "out_sps.sam"),
`known-splicesite-infile`=spsfile,
force=TRUE, strict=TRUE, execute=TRUE)To get the version of HISAT2:
hisat2_version()
#> [1] "/tmp/RtmpiCXGfL/Rinst181e6e7307a9/Rhisat2/hisat2 version 2.2.1"
#> [2] "64-bit"
#> [3] "Built on 11da3bb9ce2c"
#> [4] "Sun Dec 22 03:38:51 AM UTC 2024"
#> [5] "Compiler: gcc version 13.3.0 (Ubuntu 13.3.0-6ubuntu2~24.04) "
#> [6] "Options: -O3 -m64 -msse2 -funroll-loops -g3 -std=c++11"
#> [7] "Sizeof {int, long, long long, void*, size_t, off_t}: {4, 8, 8, 8, 8, 8}"To see the usage of hisat2_build():
hisat2_build_usage()
#> [1] "HISAT2 version 2.2.1 by Daehwan Kim ([email protected], http://www.ccb.jhu.edu/people/infphilo)"
#> [2] "Usage: hisat2-build [options]* <reference_in> <ht2_index_base>"
#> [3] " reference_in comma-separated list of files with ref sequences"
#> [4] " hisat2_index_base write ht2 data to files with this dir/basename"
#> [5] "Options:"
#> [6] " -c reference sequences given on cmd line (as"
#> [7] " <reference_in>)"
#> [8] " -a/--noauto disable automatic -p/--bmax/--dcv memory-fitting"
#> [9] " -p <int> number of threads"
#> [10] " --bmax <int> max bucket sz for blockwise suffix-array builder"
#> [11] " --bmaxdivn <int> max bucket sz as divisor of ref len (default: 4)"
#> [12] " --dcv <int> diff-cover period for blockwise (default: 1024)"
#> [13] " --nodc disable diff-cover (algorithm becomes quadratic)"
#> [14] " -r/--noref don't build .3/.4.ht2 (packed reference) portion"
#> [15] " -3/--justref just build .3/.4.ht2 (packed reference) portion"
#> [16] " -o/--offrate <int> SA is sampled every 2^offRate BWT chars (default: 5)"
#> [17] " -t/--ftabchars <int> # of chars consumed in initial lookup (default: 10)"
#> [18] " --localoffrate <int> SA (local) is sampled every 2^offRate BWT chars (default: 3)"
#> [19] " --localftabchars <int> # of chars consumed in initial lookup in a local index (default: 6)"
#> [20] " --snp <path> SNP file name"
#> [21] " --haplotype <path> haplotype file name"
#> [22] " --ss <path> Splice site file name"
#> [23] " --exon <path> Exon file name"
#> [24] " --repeat-ref <path> Repeat reference file name"
#> [25] " --repeat-info <path> Repeat information file name"
#> [26] " --repeat-snp <path> Repeat snp file name"
#> [27] " --repeat-haplotype <path> Repeat haplotype file name"
#> [28] " --seed <int> seed for random number generator"
#> [29] " -q/--quiet disable verbose output (for debugging)"
#> [30] " -h/--help print detailed description of tool and its options"
#> [31] " --usage print this usage message"
#> [32] " --version print version information and quit"And similarly to see the usage of hisat2():
hisat2_usage()
#> [1] "HISAT2 version 2.2.1 by Daehwan Kim ([email protected], www.ccb.jhu.edu/people/infphilo)"
#> [2] "Usage: "
#> [3] " hisat2-align [options]* -x <ht2-idx> {-1 <m1> -2 <m2> | -U <r>} [-S <sam>]"
#> [4] ""
#> [5] " <ht2-idx> Index filename prefix (minus trailing .X.ht2)."
#> [6] " <m1> Files with #1 mates, paired with files in <m2>."
#> [7] " <m2> Files with #2 mates, paired with files in <m1>."
#> [8] " <r> Files with unpaired reads."
#> [9] " <sam> File for SAM output (default: stdout)"
#> [10] ""
#> [11] " <m1>, <m2>, <r> can be comma-separated lists (no whitespace) and can be"
#> [12] " specified many times. E.g. '-U file1.fq,file2.fq -U file3.fq'."
#> [13] ""
#> [14] "Options (defaults in parentheses):"
#> [15] ""
#> [16] " Input:"
#> [17] " -q query input files are FASTQ .fq/.fastq (default)"
#> [18] " --qseq query input files are in Illumina's qseq format"
#> [19] " -f query input files are (multi-)FASTA .fa/.mfa"
#> [20] " -r query input files are raw one-sequence-per-line"
#> [21] " -c <m1>, <m2>, <r> are sequences themselves, not files"
#> [22] " -s/--skip <int> skip the first <int> reads/pairs in the input (none)"
#> [23] " -u/--upto <int> stop after first <int> reads/pairs (no limit)"
#> [24] " -5/--trim5 <int> trim <int> bases from 5'/left end of reads (0)"
#> [25] " -3/--trim3 <int> trim <int> bases from 3'/right end of reads (0)"
#> [26] " --phred33 qualities are Phred+33 (default)"
#> [27] " --phred64 qualities are Phred+64"
#> [28] " --int-quals qualities encoded as space-delimited integers"
#> [29] ""
#> [30] " Presets: Same as:"
#> [31] " --fast --no-repeat-index"
#> [32] " --sensitive --bowtie2-dp 1 -k 30 --score-min L,0,-0.5"
#> [33] " --very-sensitive --bowtie2-dp 2 -k 50 --score-min L,0,-1"
#> [34] ""
#> [35] " Alignment:"
#> [36] " --bowtie2-dp <int> use Bowtie2's dynamic programming alignment algorithm (0) - 0: no dynamic programming, 1: conditional dynamic programming, and 2: unconditional dynamic programming (slowest)"
#> [37] " --n-ceil <func> func for max # non-A/C/G/Ts permitted in aln (L,0,0.15)"
#> [38] " --ignore-quals treat all quality values as 30 on Phred scale (off)"
#> [39] " --nofw do not align forward (original) version of read (off)"
#> [40] " --norc do not align reverse-complement version of read (off)"
#> [41] " --no-repeat-index do not use repeat index"
#> [42] ""
#> [43] " Spliced Alignment:"
#> [44] " --pen-cansplice <int> penalty for a canonical splice site (0)"
#> [45] " --pen-noncansplice <int> penalty for a non-canonical splice site (12)"
#> [46] " --pen-canintronlen <func> penalty for long introns (G,-8,1) with canonical splice sites"
#> [47] " --pen-noncanintronlen <func> penalty for long introns (G,-8,1) with noncanonical splice sites"
#> [48] " --min-intronlen <int> minimum intron length (20)"
#> [49] " --max-intronlen <int> maximum intron length (500000)"
#> [50] " --known-splicesite-infile <path> provide a list of known splice sites"
#> [51] " --novel-splicesite-outfile <path> report a list of splice sites"
#> [52] " --novel-splicesite-infile <path> provide a list of novel splice sites"
#> [53] " --no-temp-splicesite disable the use of splice sites found"
#> [54] " --no-spliced-alignment disable spliced alignment"
#> [55] " --rna-strandness <string> specify strand-specific information (unstranded)"
#> [56] " --tmo reports only those alignments within known transcriptome"
#> [57] " --dta reports alignments tailored for transcript assemblers"
#> [58] " --dta-cufflinks reports alignments tailored specifically for cufflinks"
#> [59] " --avoid-pseudogene tries to avoid aligning reads to pseudogenes (experimental option)"
#> [60] " --no-templatelen-adjustment disables template length adjustment for RNA-seq reads"
#> [61] ""
#> [62] " Scoring:"
#> [63] " --mp <int>,<int> max and min penalties for mismatch; lower qual = lower penalty <6,2>"
#> [64] " --sp <int>,<int> max and min penalties for soft-clipping; lower qual = lower penalty <2,1>"
#> [65] " --no-softclip no soft-clipping"
#> [66] " --np <int> penalty for non-A/C/G/Ts in read/ref (1)"
#> [67] " --rdg <int>,<int> read gap open, extend penalties (5,3)"
#> [68] " --rfg <int>,<int> reference gap open, extend penalties (5,3)"
#> [69] " --score-min <func> min acceptable alignment score w/r/t read length"
#> [70] " (L,0.0,-0.2)"
#> [71] ""
#> [72] " Reporting:"
#> [73] " -k <int> It searches for at most <int> distinct, primary alignments for each read. Primary alignments mean "
#> [74] " alignments whose alignment score is equal to or higher than any other alignments. The search terminates "
#> [75] " when it cannot find more distinct valid alignments, or when it finds <int>, whichever happens first. "
#> [76] " The alignment score for a paired-end alignment equals the sum of the alignment scores of "
#> [77] " the individual mates. Each reported read or pair alignment beyond the first has the SAM ‘secondary’ bit "
#> [78] " (which equals 256) set in its FLAGS field. For reads that have more than <int> distinct, "
#> [79] " valid alignments, hisat2 does not guarantee that the <int> alignments reported are the best possible "
#> [80] " in terms of alignment score. Default: 5 (linear index) or 10 (graph index)."
#> [81] " Note: HISAT2 is not designed with large values for -k in mind, and when aligning reads to long, "
#> [82] " repetitive genomes, large -k could make alignment much slower."
#> [83] " --max-seeds <int> HISAT2, like other aligners, uses seed-and-extend approaches. HISAT2 tries to extend seeds to "
#> [84] " full-length alignments. In HISAT2, --max-seeds is used to control the maximum number of seeds that "
#> [85] " will be extended. For DNA-read alignment (--no-spliced-alignment), HISAT2 extends up to these many seeds"
#> [86] " and skips the rest of the seeds. For RNA-read alignment, HISAT2 skips extending seeds and reports "
#> [87] " no alignments if the number of seeds is larger than the number specified with the option, "
#> [88] " to be compatible with previous versions of HISAT2. Large values for --max-seeds may improve alignment "
#> [89] " sensitivity, but HISAT2 is not designed with large values for --max-seeds in mind, and when aligning "
#> [90] " reads to long, repetitive genomes, large --max-seeds could make alignment much slower. "
#> [91] " The default value is the maximum of 5 and the value that comes with -k times 2."
#> [92] " -a/--all HISAT2 reports all alignments it can find. Using the option is equivalent to using both --max-seeds "
#> [93] " and -k with the maximum value that a 64-bit signed integer can represent (9,223,372,036,854,775,807)."
#> [94] " --repeat report alignments to repeat sequences directly"
#> [95] ""
#> [96] " Paired-end:"
#> [97] " -I/--minins <int> minimum fragment length (0), only valid with --no-spliced-alignment"
#> [98] " -X/--maxins <int> maximum fragment length (500), only valid with --no-spliced-alignment"
#> [99] " --fr/--rf/--ff -1, -2 mates align fw/rev, rev/fw, fw/fw (--fr)"
#> [100] " --no-mixed suppress unpaired alignments for paired reads"
#> [101] " --no-discordant suppress discordant alignments for paired reads"
#> [102] ""
#> [103] " Output:"
#> [104] " -t/--time print wall-clock time taken by search phases"
#> [105] " --summary-file <path> print alignment summary to this file."
#> [106] " --new-summary print alignment summary in a new style, which is more machine-friendly."
#> [107] " --quiet print nothing to stderr except serious errors"
#> [108] " --met-file <path> send metrics to file at <path> (off)"
#> [109] " --met-stderr send metrics to stderr (off)"
#> [110] " --met <int> report internal counters & metrics every <int> secs (1)"
#> [111] " --no-head suppress header lines, i.e. lines starting with @"
#> [112] " --no-sq suppress @SQ header lines"
#> [113] " --rg-id <text> set read group id, reflected in @RG line and RG:Z: opt field"
#> [114] " --rg <text> add <text> (\"lab:value\") to @RG line of SAM header."
#> [115] " Note: @RG line only printed when --rg-id is set."
#> [116] " --omit-sec-seq put '*' in SEQ and QUAL fields for secondary alignments."
#> [117] ""
#> [118] " Performance:"
#> [119] " -o/--offrate <int> override offrate of index; must be >= index's offrate"
#> [120] " -p/--threads <int> number of alignment threads to launch (1)"
#> [121] " --reorder force SAM output order to match order of input reads"
#> [122] " --mm use memory-mapped I/O for index; many 'hisat2's can share"
#> [123] ""
#> [124] " Other:"
#> [125] " --qc-filter filter out reads that are bad according to QSEQ filter"
#> [126] " --seed <int> seed for random number generator (0)"
#> [127] " --non-deterministic seed rand. gen. arbitrarily instead of using read attributes"
#> [128] " --remove-chrname remove 'chr' from reference names in alignment"
#> [129] " --add-chrname add 'chr' to reference names in alignment "
#> [130] " --version print version information and quit"
#> [131] " -h/--help print this usage message"sessionInfo()
#> R version 4.4.2 (2024-10-31)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8
#> [4] LC_COLLATE=C LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C
#> [10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] Rhisat2_1.23.0 BiocStyle_2.35.0
#>
#> loaded via a namespace (and not attached):
#> [1] tidyselect_1.2.1 dplyr_1.1.4 blob_1.2.4
#> [4] filelock_1.0.3 Biostrings_2.75.3 bitops_1.0-9
#> [7] fastmap_1.2.0 RCurl_1.98-1.16 BiocFileCache_2.15.0
#> [10] GenomicAlignments_1.43.0 XML_3.99-0.17 digest_0.6.37
#> [13] lifecycle_1.0.4 KEGGREST_1.47.0 RSQLite_2.3.9
#> [16] magrittr_2.0.3 compiler_4.4.2 rlang_1.1.4
#> [19] sass_0.4.9 progress_1.2.3 tools_4.4.2
#> [22] igraph_2.1.2 yaml_2.3.10 rtracklayer_1.67.0
#> [25] knitr_1.49 prettyunits_1.2.0 S4Arrays_1.7.1
#> [28] bit_4.5.0.1 curl_6.0.1 DelayedArray_0.33.3
#> [31] xml2_1.3.6 abind_1.4-8 BiocParallel_1.41.0
#> [34] txdbmaker_1.3.1 BiocGenerics_0.53.3 sys_3.4.3
#> [37] grid_4.4.2 stats4_4.4.2 biomaRt_2.63.0
#> [40] SummarizedExperiment_1.37.0 cli_3.6.3 rmarkdown_2.29
#> [43] crayon_1.5.3 generics_0.1.3 httr_1.4.7
#> [46] rjson_0.2.23 RUnit_0.4.33 DBI_1.2.3
#> [49] cachem_1.1.0 stringr_1.5.1 zlibbioc_1.52.0
#> [52] parallel_4.4.2 AnnotationDbi_1.69.0 BiocManager_1.30.25
#> [55] XVector_0.47.1 restfulr_0.0.15 matrixStats_1.4.1
#> [58] vctrs_0.6.5 Matrix_1.7-1 jsonlite_1.8.9
#> [61] IRanges_2.41.2 hms_1.1.3 S4Vectors_0.45.2
#> [64] bit64_4.5.2 GenomicFeatures_1.59.1 maketools_1.3.1
#> [67] SGSeq_1.41.0 jquerylib_0.1.4 glue_1.8.0
#> [70] codetools_0.2-20 stringi_1.8.4 GenomeInfoDb_1.43.2
#> [73] GenomicRanges_1.59.1 BiocIO_1.17.1 UCSC.utils_1.3.0
#> [76] tibble_3.2.1 pillar_1.10.0 rappdirs_0.3.3
#> [79] htmltools_0.5.8.1 GenomeInfoDbData_1.2.13 R6_2.5.1
#> [82] dbplyr_2.5.0 httr2_1.0.7 evaluate_1.0.1
#> [85] lattice_0.22-6 Biobase_2.67.0 png_0.1-8
#> [88] Rsamtools_2.23.1 memoise_2.0.1 bslib_0.8.0
#> [91] SparseArray_1.7.2 xfun_0.49 MatrixGenerics_1.19.0
#> [94] buildtools_1.0.0 pkgconfig_2.0.3