diff --git a/lib/galaxy/datatypes/binary.py b/lib/galaxy/datatypes/binary.py
index 8148133decf..daefe7a7e97 100644
--- a/lib/galaxy/datatypes/binary.py
+++ b/lib/galaxy/datatypes/binary.py
@@ -16,12 +16,13 @@ from json import dumps
 
 import h5py
 import pysam
+import pysam.bcftools
 from bx.seq.twobit import TWOBIT_MAGIC_NUMBER, TWOBIT_MAGIC_NUMBER_SWAP, TWOBIT_MAGIC_SIZE
 
 from galaxy import util
 from galaxy.datatypes import metadata
 from galaxy.datatypes.metadata import DictParameter, ListParameter, MetadataElement, MetadataParameter
-from galaxy.util import FILENAME_VALID_CHARS, nice_size, sqlite, which
+from galaxy.util import FILENAME_VALID_CHARS, nice_size, sqlite
 from galaxy.util.checkers import is_bz2, is_gzip
 from . import data, dataproviders
 
@@ -206,186 +207,85 @@ class Bam(Binary):
     MetadataElement(name="column_types", default=['str', 'int', 'str', 'int', 'int', 'str', 'str', 'int', 'int', 'str', 'str', 'str'], desc="Column types", param=metadata.ColumnTypesParameter, readonly=True, visible=False, no_value=[])
     MetadataElement(name="column_names", default=['QNAME', 'FLAG', 'RNAME', 'POS', 'MAPQ', 'CIGAR', 'MRNM', 'MPOS', 'ISIZE', 'SEQ', 'QUAL', 'OPT'], desc="Column names", readonly=True, visible=False, optional=True, no_value=[])
 
-    def _get_samtools_version(self):
-        version = '0.0.0'
-        samtools_exec = which('samtools')
-        if not samtools_exec:
-            message = 'Attempting to use functionality requiring samtools, but it cannot be located on Galaxy\'s PATH.'
-            raise Exception(message)
-
-        p = subprocess.Popen(['samtools', '--version-only'], stdout=subprocess.PIPE, stderr=subprocess.PIPE)
-        output, error = p.communicate()
-        # --version-only is available
-        # Format is <version x.y.z>+htslib-<a.b.c>
-        if p.returncode == 0:
-            version = output.split('+')[0]
-            return version
-
-        output = subprocess.Popen(['samtools'], stderr=subprocess.PIPE, stdout=subprocess.PIPE).communicate()[1]
-        lines = output.split('\n')
-        for line in lines:
-            if line.lower().startswith('version'):
-                # Assuming line looks something like: version: 0.1.12a (r862)
-                version = line.split()[1]
-                break
-        return version
-
     @staticmethod
     def merge(split_files, output_file):
+        """
+        Merges BAM files
 
-        tmp_dir = tempfile.mkdtemp()
-        stderr_name = tempfile.NamedTemporaryFile(dir=tmp_dir, prefix="bam_merge_stderr").name
-        command = ["samtools", "merge", "-f", output_file] + split_files
-        proc = subprocess.Popen(args=command, stderr=open(stderr_name, 'wb'))
-        exit_code = proc.wait()
-        # Did merge succeed?
-        stderr = open(stderr_name).read().strip()
-        if stderr:
-            if exit_code != 0:
-                shutil.rmtree(tmp_dir)  # clean up
-                raise Exception("Error merging BAM files: %s" % stderr)
-            else:
-                print(stderr)
-        os.unlink(stderr_name)
-        os.rmdir(tmp_dir)
-
-    def _is_coordinate_sorted(self, file_name):
-        """See if the input BAM file is sorted from the header information."""
-        output = subprocess.check_output(["samtools", "view", "-H", file_name])
-        return 'SO:coordinate' in output or 'SO:sorted' in output
+        :param split_files: List of bam file paths to merge
+        :param output_file: Write merged bam file to this location
+        """
+        pysam.merge('-O', 'BAM', output_file, *split_files)
 
     def dataset_content_needs_grooming(self, file_name):
-        """See if file_name is a sorted BAM file"""
-        version = self._get_samtools_version()
-        if version < '0.1.13':
-            return not self._is_coordinate_sorted(file_name)
-        else:
-            # Samtools version 0.1.13 or newer produces an error condition when attempting to index an
-            # unsorted bam file - see http://biostar.stackexchange.com/questions/5273/is-my-bam-file-sorted.
-            # So when using a newer version of samtools, we'll first check if the input BAM file is sorted
-            # from the header information.  If the header is present and sorted, we do nothing by returning False.
-            # If it's present and unsorted or if it's missing, we'll index the bam file to see if it produces the
-            # error.  If it does, sorting is needed so we return True (otherwise False).
-            #
-            # TODO: we're creating an index file here and throwing it away.  We then create it again when
-            # the set_meta() method below is called later in the job process.  We need to enhance this overall
-            # process so we don't create an index twice.  In order to make it worth the time to implement the
-            # upload tool / framework to allow setting metadata from directly within the tool itself, it should be
-            # done generically so that all tools will have the ability.  In testing, a 6.6 gb BAM file took 128
-            # seconds to index with samtools, and 45 minutes to sort, so indexing is relatively inexpensive.
-            if self._is_coordinate_sorted(file_name):
-                return False
-            index_name = tempfile.NamedTemporaryFile(prefix="bam_index").name
-            stderr_name = tempfile.NamedTemporaryFile(prefix="bam_index_stderr").name
-            proc = subprocess.Popen(['samtools', 'index', file_name, index_name], stderr=open(stderr_name, 'wb'))
-            proc.wait()
-            stderr = open(stderr_name).read().strip()
-            if stderr:
-                try:
-                    os.unlink(index_name)
-                except OSError:
-                    pass
-                try:
-                    os.unlink(stderr_name)
-                except OSError:
-                    pass
-                # Return True if unsorted error condition is found (find returns -1 if string is not found).
-                return stderr.find("[bam_index_core] the alignment is not sorted") != -1
-            try:
-                os.unlink(index_name)
-            except OSError:
-                pass
-            try:
-                os.unlink(stderr_name)
-            except OSError:
-                pass
-            return False
+        """
+        Check if file_name is a coordinate-sorted BAM file
+        """
+        # The best way to ensure that BAM files are coordinate-sorted and indexable
+        # is to actually index them.
+        index_name = tempfile.NamedTemporaryFile(prefix="bam_index").name
+        try:
+            # If pysam fails to index a file it will write to stderr,
+            # and this causes the set_meta script to fail. So instead
+            # we start another process and discard stderr.
+            cmd = ['python', '-c', "import pysam; pysam.index('%s', '%s')" % (file_name, index_name)]
+            with open(os.devnull, 'w') as devnull:
+                subprocess.check_call(cmd, stderr=devnull, shell=False)
+            needs_sorting = False
+        except subprocess.CalledProcessError:
+            needs_sorting = True
+        try:
+            os.unlink(index_name)
+        except Exception:
+            pass
+        return needs_sorting
 
     def groom_dataset_content(self, file_name):
         """
-        Ensures that the Bam file contents are sorted.  This function is called
+        Ensures that the BAM file contents are sorted.  This function is called
         on an output dataset after the content is initially generated.
         """
-        # Use samtools to sort the Bam file
-        # $ samtools sort
-        # Usage: samtools sort [-on] [-m <maxMem>] <in.bam> <out.prefix>
-        # Sort alignments by leftmost coordinates. File <out.prefix>.bam will be created.
-        # This command may also create temporary files <out.prefix>.%d.bam when the
-        # whole alignment cannot be fitted into memory ( controlled by option -m ).
+        # Use pysam to sort the BAM file
+        # This command may also creates temporary files <out.prefix>.%d.bam when the
+        # whole alignment cannot fit into memory.
         # do this in a unique temp directory, because of possible <out.prefix>.%d.bam temp files
         if not self.dataset_content_needs_grooming(file_name):
             # Don't re-sort if already sorted
             return
         tmp_dir = tempfile.mkdtemp()
         tmp_sorted_dataset_file_name_prefix = os.path.join(tmp_dir, 'sorted')
-        stderr_name = tempfile.NamedTemporaryFile(dir=tmp_dir, prefix="bam_sort_stderr").name
-        samtools_created_sorted_file_name = "%s.bam" % tmp_sorted_dataset_file_name_prefix  # samtools accepts a prefix, not a filename, it always adds .bam to the prefix
-        proc = subprocess.Popen(['samtools', 'sort', file_name, tmp_sorted_dataset_file_name_prefix],
-                                cwd=tmp_dir, stderr=open(stderr_name, 'wb'))
-        exit_code = proc.wait()
-        # Did sort succeed?
-        stderr = open(stderr_name).read().strip()
-        if stderr:
-            if exit_code != 0:
-                shutil.rmtree(tmp_dir)  # clean up
-                raise Exception("Error Grooming BAM file contents: %s" % stderr)
-            else:
-                print(stderr)
+        sorted_file_name = "%s.bam" % tmp_sorted_dataset_file_name_prefix
+        slots = os.environ.get('GALAXY_SLOTS', 1)
+        try:
+            pysam.sort("-@%s" % slots, file_name, '-T', tmp_sorted_dataset_file_name_prefix, '-O', 'BAM', '-o', sorted_file_name)
+        except Exception:
+            shutil.rmtree(tmp_dir, ignore_errors=True)
+            raise
         # Move samtools_created_sorted_file_name to our output dataset location
-        shutil.move(samtools_created_sorted_file_name, file_name)
+        shutil.move(sorted_file_name, file_name)
         # Remove temp file and empty temporary directory
-        os.unlink(stderr_name)
         os.rmdir(tmp_dir)
 
     def init_meta(self, dataset, copy_from=None):
         Binary.init_meta(self, dataset, copy_from=copy_from)
 
     def set_meta(self, dataset, overwrite=True, **kwd):
-        """ Creates the index for the BAM file. """
         # These metadata values are not accessible by users, always overwrite
         index_file = dataset.metadata.bam_index
         if not index_file:
             index_file = dataset.metadata.spec['bam_index'].param.new_file(dataset=dataset)
-        # Create the Bam index
-        # $ samtools index
-        # Usage: samtools index <in.bam> [<out.index>]
-        stderr_name = tempfile.NamedTemporaryFile(prefix="bam_index_stderr").name
-        command = ['samtools', 'index', dataset.file_name, index_file.file_name]
-        exit_code = subprocess.call(args=command, stderr=open(stderr_name, 'wb'))
-        # Did index succeed?
-        if exit_code == -6:
-            # SIGABRT, most likely samtools 1.0+ which does not accept the index name parameter.
-            dataset_symlink = os.path.join(os.path.dirname(index_file.file_name),
-                                           '__dataset_%d_%s' % (dataset.id, os.path.basename(index_file.file_name)))
-            os.symlink(dataset.file_name, dataset_symlink)
-            try:
-                command = ['samtools', 'index', dataset_symlink]
-                exit_code = subprocess.call(args=command, stderr=open(stderr_name, 'wb'))
-                shutil.move(dataset_symlink + '.bai', index_file.file_name)
-            except Exception as e:
-                open(stderr_name, 'ab+').write('Galaxy attempted to build the BAM index with samtools 1.0+ but failed: %s\n' % e)
-                exit_code = 1  # Make sure an exception raised by shutil.move() is re-raised below
-            finally:
-                os.unlink(dataset_symlink)
-        stderr = open(stderr_name).read().strip()
-        if stderr:
-            if exit_code != 0:
-                os.unlink(stderr_name)  # clean up
-                raise Exception("Error Setting BAM Metadata: %s" % stderr)
-            else:
-                print(stderr)
+        pysam.index(dataset.file_name, index_file.file_name)
         dataset.metadata.bam_index = index_file
-        # Remove temp file
-        os.unlink(stderr_name)
         # Now use pysam with BAI index to determine additional metadata
         try:
             bam_file = pysam.AlignmentFile(dataset.file_name, mode='rb', index_filename=index_file.file_name)
+            # TODO: Reference names, lengths, read_groups and headers can become very large, truncate when necessary
             dataset.metadata.reference_names = list(bam_file.references)
             dataset.metadata.reference_lengths = list(bam_file.lengths)
             dataset.metadata.bam_header = bam_file.header
             dataset.metadata.read_groups = [read_group['ID'] for read_group in dataset.metadata.bam_header.get('RG', []) if 'ID' in read_group]
-            dataset.metadata.sort_order = dataset.metadata.bam_header.get('HD', {}).get('SO', None)
-            dataset.metadata.bam_version = dataset.metadata.bam_header.get('HD', {}).get('VN', None)
+            dataset.metadata.sort_order = bam_file.header.get('HD', {}).get('SO', None)
+            dataset.metadata.bam_version = bam_file.header.get('HD', {}).get('VN', None)
         except Exception:
             # Per Dan, don't log here because doing so will cause datasets that
             # fail metadata to end in the error state
@@ -586,24 +486,7 @@ class CRAM(Binary):
 
     def set_index_file(self, dataset, index_file):
         try:
-            # @todo when pysam 1.2.1 or pysam 1.3.0 gets released and becomes
-            # a dependency of galaxy, use pysam.index(alignment, target_idx)
-            # This currently gives coredump in the current release but is
-            # fixed in the dev branch:
-            # xref: https://github.com/samtools/samtools/issues/199
-
-            dataset_symlink = os.path.join(os.path.dirname(index_file.file_name), '__dataset_%d_%s' % (dataset.id, os.path.basename(index_file.file_name)))
-            os.symlink(dataset.file_name, dataset_symlink)
-            pysam.index(dataset_symlink)
-
-            tmp_index = dataset_symlink + ".crai"
-            if os.path.isfile(tmp_index):
-                shutil.move(tmp_index, index_file.file_name)
-                return index_file.file_name
-            else:
-                os.unlink(dataset_symlink)
-                log.warning('%s, expected crai index not created for: %s', self, dataset.file_name)
-                return False
+            pysam.index(dataset.file_name, index_file.file_name)
         except Exception as exc:
             log.warning('%s, set_index_file Exception: %s', self, exc)
             return False
@@ -635,13 +518,6 @@ class Bcf(BaseBcf):
     """
     Class describing a (BGZF-compressed) BCF file
 
-    >>> from galaxy.datatypes.sniff import get_test_fname
-    >>> fname = get_test_fname('1.bcf')
-    >>> Bcf().sniff(fname)
-    True
-    >>> fname = get_test_fname('1.bcf_uncompressed')
-    >>> Bcf().sniff(fname)
-    False
     """
     file_ext = "bcf"
 
@@ -665,24 +541,16 @@ class Bcf(BaseBcf):
         if not index_file:
             index_file = dataset.metadata.spec['bcf_index'].param.new_file(dataset=dataset)
         # Create the bcf index
-        # $ bcftools index
-        # Usage: bcftools index <in.bcf>
-
         dataset_symlink = os.path.join(os.path.dirname(index_file.file_name),
                                        '__dataset_%d_%s' % (dataset.id, os.path.basename(index_file.file_name)))
         os.symlink(dataset.file_name, dataset_symlink)
-
-        stderr_name = tempfile.NamedTemporaryFile(prefix="bcf_index_stderr").name
-        command = ['bcftools', 'index', dataset_symlink]
         try:
-            subprocess.check_call(args=command, stderr=open(stderr_name, 'wb'))
-            shutil.move(dataset_symlink + '.csi', index_file.file_name)  # this will fail if bcftools < 1.0 is used, because it creates a .bci index file instead of .csi
+            pysam.bcftools.index(dataset_symlink)
+            shutil.move(dataset_symlink + '.csi', index_file.file_name)
         except Exception as e:
-            stderr = open(stderr_name).read().strip()
-            raise Exception('Error setting BCF metadata: %s' % (stderr or str(e)))
+            raise Exception('Error setting BCF metadata: %s' % (str(e)))
         finally:
             # Remove temp file and symlink
-            os.remove(stderr_name)
             os.remove(dataset_symlink)
         dataset.metadata.bcf_index = index_file
 
diff --git a/lib/galaxy/datatypes/converters/cram_to_bam.py b/lib/galaxy/datatypes/converters/cram_to_bam.py
index 6f2c6142601..c0634d77e3b 100644
--- a/lib/galaxy/datatypes/converters/cram_to_bam.py
+++ b/lib/galaxy/datatypes/converters/cram_to_bam.py
@@ -5,6 +5,7 @@ Uses pysam to convert a CRAM file to a sorted bam file.
 usage: %prog in_file out_file
 """
 import optparse
+import os
 
 import pysam
 
@@ -14,7 +15,8 @@ def main():
     parser = optparse.OptionParser()
     (options, args) = parser.parse_args()
     input_fname, output_fname = args
-    pysam.sort('-o', output_fname, '-O', 'bam', '-T', '.', input_fname)
+    slots = os.getenv('GALAXY_SLOTS', 1)
+    pysam.sort("-@%s" % slots, '-o', output_fname, '-O', 'bam', '-T', '.', input_fname)
 
 
 if __name__ == "__main__":
diff --git a/lib/galaxy/datatypes/converters/interval_to_tabix_converter.py b/lib/galaxy/datatypes/converters/interval_to_tabix_converter.py
index a8bee5eafc4..faf988da48e 100644
--- a/lib/galaxy/datatypes/converters/interval_to_tabix_converter.py
+++ b/lib/galaxy/datatypes/converters/interval_to_tabix_converter.py
@@ -21,20 +21,29 @@ def main():
     parser.add_option('-e', '--end-col', type='int', dest='end_col')
     parser.add_option('-P', '--preset', dest='preset')
     (options, args) = parser.parse_args()
-    input_fname, index_fname, out_fname = args
+    _, bgzip_fname, out_fname = args
+    to_tabix(bgzip_fname=bgzip_fname,
+             out_fname=out_fname,
+             preset=options.preset,
+             chrom_col=options.chrom_col,
+             start_col=options.start_col,
+             end_col=options.end_col)
 
+
+def to_tabix(bgzip_fname, out_fname, preset=None, chrom_col=None, start_col=None, end_col=None):
     # Create index.
-    if options.preset:
+    if preset:
         # Preset type.
-        pysam.tabix_index(filename=index_fname, preset=options.preset, keep_original=True,
-                          index_filename=out_fname)
+        bgzip_fname = pysam.tabix_index(filename=bgzip_fname, preset=preset, keep_original=True,
+                                        index=out_fname, force=True)
     else:
         # For interval files; column indices are 0-based.
-        pysam.tabix_index(filename=index_fname, seq_col=(options.chrom_col - 1),
-                          start_col=(options.start_col - 1), end_col=(options.end_col - 1),
-                          keep_original=True, index_filename=out_fname)
-    if os.path.getsize(index_fname) == 0:
+        bgzip_fname = pysam.tabix_index(filename=bgzip_fname, seq_col=(chrom_col - 1),
+                                        start_col=(start_col - 1), end_col=(end_col - 1),
+                                        keep_original=True, index=out_fname, force=True)
+    if os.path.getsize(out_fname) == 0:
         sys.stderr.write("The converted tabix index file is empty, meaning the input data is invalid.")
+    return bgzip_fname
 
 
 if __name__ == "__main__":
diff --git a/lib/galaxy/datatypes/set_metadata_tool.xml b/lib/galaxy/datatypes/set_metadata_tool.xml
index 4daaed1d23f..69c84dd9c5c 100644
--- a/lib/galaxy/datatypes/set_metadata_tool.xml
+++ b/lib/galaxy/datatypes/set_metadata_tool.xml
@@ -1,7 +1,6 @@
 <tool id="__SET_METADATA__" name="Set External Metadata" version="1.0.1" tool_type="set_metadata">
     <type class="SetMetadataTool" module="galaxy.tools"/>
     <requirements>
-        <requirement type="package">samtools</requirement>
         <requirement type="package" version="1.5">bcftools</requirement>
     </requirements>
     <action module="galaxy.tools.actions.metadata" class="SetMetadataToolAction"/>
diff --git a/lib/galaxy/datatypes/tabular.py b/lib/galaxy/datatypes/tabular.py
index 7f94565fb89..51886d7b1e3 100644
--- a/lib/galaxy/datatypes/tabular.py
+++ b/lib/galaxy/datatypes/tabular.py
@@ -15,6 +15,8 @@ import tempfile
 from cgi import escape
 from json import dumps
 
+import pysam
+
 from galaxy import util
 from galaxy.datatypes import binary, data, metadata
 from galaxy.datatypes.metadata import MetadataElement
@@ -23,6 +25,7 @@ from galaxy.datatypes.sniff import (
     iter_headers
 )
 from galaxy.util import compression_utils
+from galaxy.util.checkers import is_gzip
 from . import dataproviders
 
 if sys.version_info > (3,):
@@ -733,6 +736,11 @@ class BaseVcf(Tabular):
 class Vcf(BaseVcf):
     file_ext = 'vcf'
 
+    def sniff(self, filename):
+        if is_gzip(filename):
+            return False
+        return super(Vcf, self).sniff(filename)
+
 
 class VcfGz(BaseVcf, binary.Binary):
     file_ext = 'vcf_bgzip'
@@ -740,33 +748,23 @@ class VcfGz(BaseVcf, binary.Binary):
 
     MetadataElement(name="tabix_index", desc="Vcf Index File", param=metadata.FileParameter, file_ext="tbi", readonly=True, no_value=None, visible=False, optional=True)
 
+    def sniff(self, filename):
+        if not is_gzip(filename):
+            return False
+        return super(VcfGz, self).sniff(filename)
+
     def set_meta(self, dataset, **kwd):
         super(BaseVcf, self).set_meta(dataset, **kwd)
         """ Creates the index for the VCF file. """
         # These metadata values are not accessible by users, always overwrite
-        index_file = dataset.metadata.bcf_index
+        index_file = dataset.metadata.tabix_index
         if not index_file:
             index_file = dataset.metadata.spec['tabix_index'].param.new_file(dataset=dataset)
-        # Create the bcf index
-        # $ bcftools index
-        # Usage: bcftools index <in.bcf>
 
-        dataset_symlink = os.path.join(os.path.dirname(index_file.file_name),
-                                       '__dataset_%d_%s' % (dataset.id, os.path.basename(index_file.file_name))) + ".vcf.gz"
-        os.symlink(dataset.file_name, dataset_symlink)
-
-        stderr_name = tempfile.NamedTemporaryFile(prefix="bcf_index_stderr").name
-        command = ['bcftools', 'index', '-t', dataset_symlink]
         try:
-            subprocess.check_call(args=command, stderr=open(stderr_name, 'wb'))
-            shutil.move(dataset_symlink + '.tbi', index_file.file_name)  # this will fail if bcftools < 1.0 is used, because it creates a .bci index file instead of .csi
+            pysam.tabix_index(dataset.file_name, index=index_file.file_name, preset='vcf', force=True)
         except Exception as e:
-            stderr = open(stderr_name).read().strip()
-            raise Exception('Error setting BCF metadata: %s' % (stderr or str(e)))
-        finally:
-            # Remove temp file and symlink
-            os.remove(stderr_name)
-            os.remove(dataset_symlink)
+            raise Exception('Error setting VCF.gz metadata: %s' % (str(e)))
         dataset.metadata.tabix_index = index_file
 
 
diff --git a/lib/galaxy/datatypes/test/1.bam b/lib/galaxy/datatypes/test/1.bam
index 95c65de8572..00ad79a950e 100644
Binary files a/lib/galaxy/datatypes/test/1.bam and b/lib/galaxy/datatypes/test/1.bam differ
diff --git a/lib/galaxy/datatypes/test/1.vcf b/lib/galaxy/datatypes/test/1.vcf
new file mode 100644
index 00000000000..f75da66f171
--- /dev/null
+++ b/lib/galaxy/datatypes/test/1.vcf
@@ -0,0 +1,34 @@
+##fileformat=VCFv4.1
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth; some reads may have been filtered">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	wgs2
+chr1	500	.	C	T	.	.	DP=10	GT	1/1
+chr1	1000	.	A	G	.	.	DP=10	GT	0/1
+chr1	1500	.	G	A	.	.	DP=10	GT	1/0
+chr1	2000	.	G	A	.	.	DP=10	GT	1/1
+chr1	2500	.	A	G	.	.	DP=10	GT	0/1
+chr1	3000	.	G	C	.	.	DP=10	GT	1/0
+chr1	3500	.	G	T	.	.	DP=10	GT	1/1
+chr1	4000	.	A	G	.	.	DP=10	GT	0/1
+chr1	4500	.	A	C	.	.	DP=10	GT	1/0
+chr1	5000	.	A	G	.	.	DP=10	GT	1/1
+chr1	5500	.	C	T	.	.	DP=10	GT	0/1
+chr1	6000	.	C	A	.	.	DP=10	GT	1/0
+chr1	6500	.	C	T	.	.	DP=10	GT	1/1
+chr1	7000	.	G	T	.	.	DP=10	GT	0/1
+chr1	7500	.	G	C	.	.	DP=10	GT	1/0
+chr1	8000	.	T	G	.	.	DP=10	GT	1/1
+chr1	8500	.	T	A	.	.	DP=10	GT	0/1
+chr1	9000	.	C	T	.	.	DP=10	GT	1/0
+chr1	9500	.	G	A	.	.	DP=10	GT	1/1
+chr1	10000	.	A	G	.	.	DP=10	GT	0/1
+chr1	10500	.	T	C	.	.	DP=10	GT	1/0
+chr1	11000	.	A	G	.	.	DP=10	GT	1/1
+chr1	11500	.	C	T	.	.	DP=10	GT	0/1
+chr1	12000	.	A	G	.	.	DP=10	GT	1/0
+chr1	12500	.	C	G	.	.	DP=10	GT	1/1
+chr1	13000	.	T	C	.	.	DP=10	GT	0/1
+chr1	13500	.	G	T	.	.	DP=10	GT	1/0
+chr1	14000	.	A	C	.	.	DP=10	GT	1/1
+chr1	14500	.	A	G	.	.	DP=10	GT	0/1
+chr1	15000	.	T	G	.	.	DP=10	GT	1/0
diff --git a/lib/galaxy/datatypes/test/1.vcf.gz b/lib/galaxy/datatypes/test/1.vcf.gz
new file mode 100644
index 00000000000..1b5b3146fd8
Binary files /dev/null and b/lib/galaxy/datatypes/test/1.vcf.gz differ
diff --git a/lib/galaxy/datatypes/test/2.cram b/lib/galaxy/datatypes/test/2.cram
new file mode 100644
index 00000000000..dacd2d459e3
Binary files /dev/null and b/lib/galaxy/datatypes/test/2.cram differ
diff --git a/lib/galaxy/datatypes/test/2.shuffled.bam b/lib/galaxy/datatypes/test/2.shuffled.bam
new file mode 100644
index 00000000000..72aafe572c1
Binary files /dev/null and b/lib/galaxy/datatypes/test/2.shuffled.bam differ
diff --git a/lib/galaxy/visualization/data_providers/genome.py b/lib/galaxy/visualization/data_providers/genome.py
index 1f7445104ea..d78fbccadc5 100644
--- a/lib/galaxy/visualization/data_providers/genome.py
+++ b/lib/galaxy/visualization/data_providers/genome.py
@@ -8,6 +8,9 @@ import os
 import random
 import re
 import sys
+import tempfile
+from contextlib import contextmanager
+from distutils.version import LooseVersion
 from json import loads
 
 import pysam
@@ -24,6 +27,8 @@ from galaxy.visualization.data_providers.cigar import get_ref_based_read_seq_and
 # Utility functions.
 #
 
+PYSAM_INDEX_SYMLINK_NECESSARY = LooseVersion(pysam.__version__) <= LooseVersion('0.13.0')
+
 
 def float_nan(n):
     '''
@@ -158,6 +163,7 @@ class GenomeDataProvider(BaseDataProvider):
         """
         raise Exception("Unimplemented Function")
 
+    @contextmanager
     def open_data_file(self):
         """
         Open data file for reading data.
@@ -186,17 +192,9 @@ class GenomeDataProvider(BaseDataProvider):
             dataset_type, data
         """
         start, end = int(low), int(high)
-        data_file = self.open_data_file()
-        iterator = self.get_iterator(data_file, chrom, start, end, **kwargs)
-        data = self.process_data(iterator, start_val, max_vals, start=start, end=end, **kwargs)
-        try:
-            data_file.close()
-        except AttributeError:
-            # FIXME: some data providers do not have a close function implemented.
-            # Providers without a close function include:
-            #  bx IntervalIndex
-            pass
-
+        with self.open_data_file() as data_file:
+            iterator = self.get_iterator(data_file, chrom, start, end, **kwargs)
+            data = self.process_data(iterator, start_val, max_vals, start=start, end=end, **kwargs)
         return data
 
     def get_genome_data(self, chroms_info, **kwargs):
@@ -322,9 +320,21 @@ class TabixDataProvider(FilterableMixin, GenomeDataProvider):
 
     col_name_data_attr_mapping = {4: {'index': 4, 'name': 'Score'}}
 
+    @contextmanager
     def open_data_file(self):
-        return pysam.Tabixfile(self.dependencies['bgzip'].file_name,
-                               index=self.converted_dataset.file_name)
+        # We create a symlink to the index file. This is
+        # required until https://github.com/pysam-developers/pysam/pull/586 is merged.
+        if PYSAM_INDEX_SYMLINK_NECESSARY:
+            fd, index_path = tempfile.mkstemp(suffix='.tbi')
+            os.close(fd)
+            os.unlink(index_path)
+            os.symlink(self.converted_dataset.file_name, index_path)
+        else:
+            index_path = self.converted_dataset.file_name
+        with pysam.TabixFile(self.dependencies['bgzip'].file_name, index=index_path) as f:
+            yield f
+        if PYSAM_INDEX_SYMLINK_NECESSARY:
+            os.unlink(index_path)
 
     def get_iterator(self, data_file, chrom, start, end, **kwargs):
         # chrom must be a string, start/end integers.
@@ -347,19 +357,12 @@ class TabixDataProvider(FilterableMixin, GenomeDataProvider):
         return iterator
 
     def write_data_to_file(self, regions, filename):
-        out = open(filename, "w")
-
-        data_file = self.open_data_file()
-        for region in regions:
-            # Write data in region.
-            iterator = self.get_iterator(data_file, region.chrom, region.start, region.end)
-            for line in iterator:
-                out.write("%s\n" % line)
-
-        # TODO: once Pysam is updated and Tabixfile has a close() method,
-        # data_file.close()
-
-        out.close()
+        with self.open_data_file() as data_file, open(filename, 'w') as out:
+            for region in regions:
+                # Write data in region.
+                iterator = self.get_iterator(data_file, region.chrom, region.start, region.end)
+                for line in iterator:
+                    out.write("%s\n" % line)
 
 #
 # -- Interval data providers --
@@ -528,20 +531,16 @@ class BedDataProvider(GenomeDataProvider):
         return {'data': rval, 'dataset_type': self.dataset_type, 'message': message}
 
     def write_data_to_file(self, regions, filename):
-        out = open(filename, "w")
-
-        for region in regions:
-            # Write data in region.
-            chrom = region.chrom
-            start = region.start
-            end = region.end
-            data_file = self.open_data_file()
-            iterator = self.get_iterator(data_file, chrom, start, end)
-            for line in iterator:
-                out.write("%s\n" % line)
-            data_file.close()
-
-        out.close()
+        with open(filename, "w") as out:
+            for region in regions:
+                # Write data in region.
+                chrom = region.chrom
+                start = region.start
+                end = region.end
+                with self.open_data_file() as data_file:
+                    iterator = self.get_iterator(data_file, chrom, start, end)
+                    for line in iterator:
+                        out.write("%s\n" % line)
 
 
 class BedTabixDataProvider(TabixDataProvider, BedDataProvider):
@@ -569,18 +568,19 @@ class RawBedDataProvider(BedDataProvider):
         data_file.seek(0)
 
         def line_filter_iter():
-            for line in open(self.original_dataset.file_name):
-                if line.startswith("track") or line.startswith("browser"):
-                    continue
-                feature = line.split()
-                feature_chrom = feature[0]
-                feature_start = int(feature[1])
-                feature_end = int(feature[2])
-                if (chrom is not None and feature_chrom != chrom) \
-                        or (start is not None and feature_start > end) \
-                        or (end is not None and feature_end < start):
-                    continue
-                yield line
+            with open(self.original_dataset.file_name) as data_file:
+                for line in data_file:
+                    if line.startswith("track") or line.startswith("browser"):
+                        continue
+                    feature = line.split()
+                    feature_chrom = feature[0]
+                    feature_start = int(feature[1])
+                    feature_end = int(feature[2])
+                    if (chrom is not None and feature_chrom != chrom) \
+                            or (start is not None and feature_start > end) \
+                            or (end is not None and feature_end < start):
+                        continue
+                    yield line
 
         return line_filter_iter()
 
@@ -733,14 +733,12 @@ class VcfDataProvider(GenomeDataProvider):
 
     def write_data_to_file(self, regions, filename):
         out = open(filename, "w")
-        data_file = self.open_data_file()
-
-        for region in regions:
-            # Write data in region.
-            iterator = self.get_iterator(data_file, region.chrom, region.start, region.end)
-            for line in iterator:
-                out.write("%s\n" % line)
-        out.close()
+        with self.open_data_file() as data_file:
+            for region in regions:
+                # Write data in region.
+                iterator = self.get_iterator(data_file, region.chrom, region.start, region.end)
+                for line in iterator:
+                    out.write("%s\n" % line)
 
 
 class VcfTabixDataProvider(TabixDataProvider, VcfDataProvider):
@@ -759,8 +757,10 @@ class RawVcfDataProvider(VcfDataProvider):
     for large datasets.
     """
 
+    @contextmanager
     def open_data_file(self):
-        return open(self.original_dataset.file_name)
+        with open(self.original_dataset.file_name) as f:
+            yield f
 
     def get_iterator(self, data_file, chrom, start, end, **kwargs):
         # Skip comments.
@@ -856,10 +856,12 @@ class BamDataProvider(GenomeDataProvider, FilterableMixin):
         new_bamfile.close()
         bamfile.close()
 
+    @contextmanager
     def open_data_file(self):
         # Attempt to open the BAM file with index
-        return pysam.AlignmentFile(self.original_dataset.file_name, mode='rb',
-                                   index_filename=self.converted_dataset.file_name)
+        with pysam.AlignmentFile(self.original_dataset.file_name, mode='rb',
+                                 index_filename=self.converted_dataset.file_name) as f:
+            yield f
 
     def get_iterator(self, data_file, chrom, start, end, **kwargs):
         """
@@ -1268,33 +1270,30 @@ class IntervalIndexDataProvider(FilterableMixin, GenomeDataProvider):
     dataset_type = 'interval_index'
 
     def write_data_to_file(self, regions, filename):
-        source = open(self.original_dataset.file_name)
         index = Indexes(self.converted_dataset.file_name)
-        out = open(filename, 'w')
+        with open(self.original_dataset.file_name) as source, open(filename, 'w') as out:
+            for region in regions:
+                # Write data from region.
+                chrom = region.chrom
+                start = region.start
+                end = region.end
+                for start, end, offset in index.find(chrom, start, end):
+                    source.seek(offset)
 
-        for region in regions:
-            # Write data from region.
-            chrom = region.chrom
-            start = region.start
-            end = region.end
-            for start, end, offset in index.find(chrom, start, end):
-                source.seek(offset)
-
-                # HACK: write differently depending on original dataset format.
-                if self.original_dataset.ext not in ['gff', 'gff3', 'gtf']:
-                    line = source.readline()
-                    out.write(line)
-                else:
-                    reader = GFFReaderWrapper(source, fix_strand=True)
-                    feature = reader.next()
-                    for interval in feature.intervals:
-                        out.write('\t'.join(interval.fields) + '\n')
-
-        source.close()
-        out.close()
+                    # HACK: write differently depending on original dataset format.
+                    if self.original_dataset.ext not in ['gff', 'gff3', 'gtf']:
+                        line = source.readline()
+                        out.write(line)
+                    else:
+                        reader = GFFReaderWrapper(source, fix_strand=True)
+                        feature = reader.next()
+                        for interval in feature.intervals:
+                            out.write('\t'.join(interval.fields) + '\n')
 
+    @contextmanager
     def open_data_file(self):
-        return Indexes(self.converted_dataset.file_name)
+        i = Indexes(self.converted_dataset.file_name)
+        yield i
 
     def get_iterator(self, data_file, chrom, start, end, **kwargs):
         """
@@ -1309,34 +1308,31 @@ class IntervalIndexDataProvider(FilterableMixin, GenomeDataProvider):
     def process_data(self, iterator, start_val=0, max_vals=None, **kwargs):
         results = []
         message = None
-        source = open(self.original_dataset.file_name)
+        with open(self.original_dataset.file_name) as source:
+            # Build data to return. Payload format is:
+            # [ <guid/offset>, <start>, <end>, <name>, <score>, <strand>, <thick_start>,
+            #   <thick_end>, <blocks> ]
+            #
+            # First three entries are mandatory, others are optional.
+            filter_cols = loads(kwargs.get("filter_cols", "[]"))
+            no_detail = ("no_detail" in kwargs)
+            for count, val in enumerate(iterator):
+                offset = val[2]
+                if count < start_val:
+                    continue
+                if count - start_val >= max_vals:
+                    message = self.error_max_vals % (max_vals, "features")
+                    break
+                source.seek(offset)
+                # TODO: can we use column metadata to fill out payload?
 
-        #
-        # Build data to return. Payload format is:
-        # [ <guid/offset>, <start>, <end>, <name>, <score>, <strand>, <thick_start>,
-        #   <thick_end>, <blocks> ]
-        #
-        # First three entries are mandatory, others are optional.
-        #
-        filter_cols = loads(kwargs.get("filter_cols", "[]"))
-        no_detail = ("no_detail" in kwargs)
-        for count, val in enumerate(iterator):
-            offset = val[2]
-            if count < start_val:
-                continue
-            if count - start_val >= max_vals:
-                message = self.error_max_vals % (max_vals, "features")
-                break
-            source.seek(offset)
-            # TODO: can we use column metadata to fill out payload?
+                # GFF dataset.
+                reader = GFFReaderWrapper(source, fix_strand=True)
+                feature = reader.next()
+                payload = package_gff_feature(feature, no_detail, filter_cols)
+                payload.insert(0, offset)
 
-            # GFF dataset.
-            reader = GFFReaderWrapper(source, fix_strand=True)
-            feature = reader.next()
-            payload = package_gff_feature(feature, no_detail, filter_cols)
-            payload.insert(0, offset)
-
-            results.append(payload)
+                results.append(payload)
 
         return {'data': results, 'message': message}
 
diff --git a/test/unit/datatypes/converters/__init__.py b/test/unit/datatypes/converters/__init__.py
new file mode 100644
index 00000000000..e69de29bb2d
diff --git a/test/unit/datatypes/converters/test_interval_to_tabix.py b/test/unit/datatypes/converters/test_interval_to_tabix.py
new file mode 100644
index 00000000000..574de716ed5
--- /dev/null
+++ b/test/unit/datatypes/converters/test_interval_to_tabix.py
@@ -0,0 +1,17 @@
+import pysam
+
+from galaxy.datatypes.converters.interval_to_tabix_converter import to_tabix
+from ..util import (
+    get_input_files,
+    get_tmp_path
+)
+
+
+def test_to_tabix():
+    with get_input_files('1.vcf') as input_files:
+        with get_tmp_path(suffix='.tbi') as index:
+            bgzip_fname = to_tabix(input_files[0],
+                                   index,
+                                   preset='vcf')
+            f = pysam.TabixFile(bgzip_fname, index=index)
+            f.close()
diff --git a/test/unit/datatypes/test_bam.py b/test/unit/datatypes/test_bam.py
new file mode 100644
index 00000000000..c06efc18962
--- /dev/null
+++ b/test/unit/datatypes/test_bam.py
@@ -0,0 +1,48 @@
+import pysam
+
+from galaxy.datatypes.binary import Bam
+from .util import (
+    get_dataset,
+    get_input_files,
+    get_tmp_path
+)
+
+
+def test_merge_bam():
+    with get_input_files('1.bam', '1.bam') as input_files, get_tmp_path() as outpath:
+        Bam.merge(input_files, outpath)
+        alignment_count_output = int(pysam.view('-c', outpath).strip())
+        alignment_count_input = int(pysam.view('-c', input_files[0]).strip()) * 2
+        assert alignment_count_input == alignment_count_output
+
+
+def test_dataset_content_needs_grooming():
+    b = Bam()
+    with get_input_files('1.bam', '2.shuffled.bam') as input_files:
+        sorted_bam, shuffled_bam = input_files
+        assert b.dataset_content_needs_grooming(sorted_bam) is False
+        assert b.dataset_content_needs_grooming(shuffled_bam) is True
+
+
+def test_groom_dataset_content():
+    b = Bam()
+    try:
+        with get_input_files('2.shuffled.bam') as input_files:
+            b.groom_dataset_content(input_files[0])
+            assert b.dataset_content_needs_grooming(input_files[0]) is False
+    except AssertionError as e:
+        # Grooming modifies files in-place, so the md5 hash comparison has to fail
+        assert 'Unexpected change' in e.message
+        return
+    # should not reach this part of the test
+    raise Exception('Bam grooming did not occur in-place')
+
+
+def test_set_meta_presorted():
+    b = Bam()
+    with get_dataset('1.bam') as dataset:
+        b.set_meta(dataset=dataset)
+        assert dataset.metadata.sort_order == 'coordinate'
+        bam_file = pysam.AlignmentFile(dataset.file_name, mode='rb',
+                                       index_filename=dataset.metadata.bam_index.file_name)
+        assert bam_file.has_index() is True
diff --git a/test/unit/datatypes/test_bcf.py b/test/unit/datatypes/test_bcf.py
new file mode 100644
index 00000000000..f455a50f670
--- /dev/null
+++ b/test/unit/datatypes/test_bcf.py
@@ -0,0 +1,25 @@
+from galaxy.datatypes.binary import (
+    Bcf,
+    BcfUncompressed
+)
+from .util import (
+    get_dataset,
+    get_input_files,
+)
+
+
+def test_bcf_sniff():
+    bcf = Bcf()
+    bcfu = BcfUncompressed()
+    with get_input_files('1.bcf', '1.bcf_uncompressed') as input_files:
+        compressed, uncompressed = input_files
+        assert bcf.sniff(compressed) is True
+        assert bcf.sniff(uncompressed) is False
+        assert bcfu.sniff(compressed) is False
+        assert bcfu.sniff(uncompressed) is True
+
+
+def test_bcf_set_meta():
+    bcf = Bcf()
+    with get_input_files('1.bcf') as input_files, get_dataset(input_files[0], index_attr='bcf_index') as dataset:
+        bcf.set_meta(dataset)
diff --git a/test/unit/datatypes/test_cram.py b/test/unit/datatypes/test_cram.py
new file mode 100644
index 00000000000..c117ff785d2
--- /dev/null
+++ b/test/unit/datatypes/test_cram.py
@@ -0,0 +1,20 @@
+import os
+
+import pysam
+
+from galaxy.datatypes.binary import CRAM
+from .util import (
+    get_dataset,
+    get_input_files,
+)
+
+
+def test_cram():
+    c = CRAM()
+    with get_input_files('2.cram') as input_files, get_dataset(input_files[0], index_attr='cram_index') as dataset:
+        assert os.path.exists(dataset.metadata.cram_index.file_name) is False
+        c.set_index_file(dataset=dataset, index_file=dataset.metadata.cram_index)
+        assert os.path.exists(dataset.metadata.cram_index.file_name) is True
+        c.set_meta(dataset)
+        pysam.AlignmentFile(dataset.file_name, index_filename=dataset.metadata.cram_index.file_name)
+    assert dataset.metadata.cram_version == '3.0'
diff --git a/test/unit/datatypes/test_vcf.py b/test/unit/datatypes/test_vcf.py
new file mode 100644
index 00000000000..e1219c3211f
--- /dev/null
+++ b/test/unit/datatypes/test_vcf.py
@@ -0,0 +1,29 @@
+import pysam
+
+from galaxy.datatypes.tabular import (
+    Vcf,
+    VcfGz,
+)
+from .util import (
+    get_dataset,
+    get_input_files,
+)
+
+
+def test_vcf_sniff():
+    vcf = Vcf()
+    vcf_gz = VcfGz()
+    with get_input_files('1.vcf.gz', '1.vcf') as input_files:
+        compressed, uncompressed = input_files
+        assert vcf_gz.sniff(compressed) is True
+        assert vcf_gz.sniff(uncompressed) is False
+        assert vcf.sniff(compressed) is False
+        assert vcf.sniff(uncompressed) is True
+
+
+def test_vcf_gz_set_meta():
+    vcf_gz = VcfGz()
+    with get_input_files('1.vcf.gz') as input_files, get_dataset(input_files[0], index_attr='tabix_index') as dataset:
+        vcf_gz.set_meta(dataset)
+        f = pysam.VariantFile(dataset.file_name, index_filename=dataset.metadata.tabix_index.file_name)
+        assert isinstance(f.index, pysam.libcbcf.TabixIndex) is True
diff --git a/test/unit/datatypes/util.py b/test/unit/datatypes/util.py
new file mode 100644
index 00000000000..47d5be832b9
--- /dev/null
+++ b/test/unit/datatypes/util.py
@@ -0,0 +1,51 @@
+import os
+import shutil
+import tempfile
+from contextlib import contextmanager
+
+from galaxy.datatypes.sniff import get_test_fname
+from galaxy.util.bunch import Bunch
+from galaxy.util.hash_util import md5_hash_file
+
+
+@contextmanager
+def get_dataset(filename, index_attr='bam_index', dataset_id=1, has_data=True):
+    dataset = Bunch()
+    dataset.has_data = lambda: True
+    dataset.id = dataset_id
+    dataset.metadata = Bunch()
+    with get_input_files(filename) as input_files, get_tmp_path() as index_path:
+        dataset.file_name = input_files[0]
+        index = Bunch()
+        index.file_name = index_path
+        setattr(dataset.metadata, index_attr, index)
+        yield dataset
+
+
+@contextmanager
+def get_tmp_path(should_exist=False, suffix=''):
+    _, path = tempfile.mkstemp(suffix=suffix)
+    if not should_exist:
+        os.remove(path)
+    yield path
+    try:
+        os.remove(path)
+    except Exception:
+        pass
+
+
+@contextmanager
+def get_input_files(*args):
+    temp_dir = tempfile.mkdtemp()
+    test_files = []
+    try:
+        for filename in args:
+            shutil.copy(get_test_fname(filename), temp_dir)
+            test_files.append(os.path.join(temp_dir, filename))
+        md5_sums = [md5_hash_file(f) for f in test_files]
+        yield test_files
+        new_md5_sums = [md5_hash_file(f) for f in test_files]
+        for old_hash, new_hash, f in zip(md5_sums, new_md5_sums, test_files):
+            assert old_hash == new_hash, 'Unexpected change of content for file %s' % f
+    finally:
+        shutil.rmtree(temp_dir, ignore_errors=True)
diff --git a/tools/data_source/upload.xml b/tools/data_source/upload.xml
index 8d3fe1959c5..6f4233dc66d 100644
--- a/tools/data_source/upload.xml
+++ b/tools/data_source/upload.xml
@@ -5,9 +5,6 @@
     from your computer  
   </description>
   <action module="galaxy.tools.actions.upload" class="UploadToolAction"/>
-  <requirements>
-      <requirement type="package">samtools</requirement>
-  </requirements>
   <command interpreter="python">
       upload.py $GALAXY_ROOT_DIR $GALAXY_DATATYPES_CONF_FILE $paramfile
     #set $outnum = 0