.\" Automatically generated by Pod::Man 4.14 (Pod::Simple 3.40) .\" .\" Standard preamble: .\" ======================================================================== .de Sp \" Vertical space (when we can't use .PP) .if t .sp .5v .if n .sp .. .de Vb \" Begin verbatim text .ft CW .nf .ne \\$1 .. .de Ve \" End verbatim text .ft R .fi .. .\" Set up some character translations and predefined strings. \*(-- will .\" give an unbreakable dash, \*(PI will give pi, \*(L" will give a left .\" double quote, and \*(R" will give a right double quote. \*(C+ will .\" give a nicer C++. Capital omega is used to do unbreakable dashes and .\" therefore won't be available. \*(C` and \*(C' expand to `' in nroff, .\" nothing in troff, for use with C<>. .tr \(*W- .ds C+ C\v'-.1v'\h'-1p'\s-2+\h'-1p'+\s0\v'.1v'\h'-1p' .ie n \{\ . ds -- \(*W- . ds PI pi . if (\n(.H=4u)&(1m=24u) .ds -- \(*W\h'-12u'\(*W\h'-12u'-\" diablo 10 pitch . if (\n(.H=4u)&(1m=20u) .ds -- \(*W\h'-12u'\(*W\h'-8u'-\" diablo 12 pitch . ds L" "" . ds R" "" . ds C` "" . ds C' "" 'br\} .el\{\ . ds -- \|\(em\| . ds PI \(*p . ds L" `` . ds R" '' . ds C` . ds C' 'br\} .\" .\" Escape single quotes in literal strings from groff's Unicode transform. .ie \n(.g .ds Aq \(aq .el .ds Aq ' .\" .\" If the F register is >0, we'll generate index entries on stderr for .\" titles (.TH), headers (.SH), subsections (.SS), items (.Ip), and index .\" entries marked with X<> in POD. Of course, you'll have to process the .\" output yourself in some meaningful fashion. .\" .\" Avoid warning from groff about undefined register 'F'. .de IX .. .nr rF 0 .if \n(.g .if rF .nr rF 1 .if (\n(rF:(\n(.g==0)) \{\ . if \nF \{\ . de IX . tm Index:\\$1\t\\n%\t"\\$2" .. . if !\nF==2 \{\ . nr % 0 . nr F 2 . \} . \} .\} .rr rF .\" .\" Accent mark definitions (@(#)ms.acc 1.5 88/02/08 SMI; from UCB 4.2). .\" Fear. Run. Save yourself. No user-serviceable parts. . \" fudge factors for nroff and troff .if n \{\ . ds #H 0 . ds #V .8m . ds #F .3m . ds #[ \f1 . ds #] \fP .\} .if t \{\ . ds #H ((1u-(\\\\n(.fu%2u))*.13m) . ds #V .6m . ds #F 0 . ds #[ \& . ds #] \& .\} . \" simple accents for nroff and troff .if n \{\ . ds ' \& . ds ` \& . ds ^ \& . ds , \& . ds ~ ~ . ds / .\} .if t \{\ . ds ' \\k:\h'-(\\n(.wu*8/10-\*(#H)'\'\h"|\\n:u" . ds ` \\k:\h'-(\\n(.wu*8/10-\*(#H)'\`\h'|\\n:u' . ds ^ \\k:\h'-(\\n(.wu*10/11-\*(#H)'^\h'|\\n:u' . ds , \\k:\h'-(\\n(.wu*8/10)',\h'|\\n:u' . ds ~ \\k:\h'-(\\n(.wu-\*(#H-.1m)'~\h'|\\n:u' . ds / \\k:\h'-(\\n(.wu*8/10-\*(#H)'\z\(sl\h'|\\n:u' .\} . \" troff and (daisy-wheel) nroff accents .ds : \\k:\h'-(\\n(.wu*8/10-\*(#H+.1m+\*(#F)'\v'-\*(#V'\z.\h'.2m+\*(#F'.\h'|\\n:u'\v'\*(#V' .ds 8 \h'\*(#H'\(*b\h'-\*(#H' .ds o \\k:\h'-(\\n(.wu+\w'\(de'u-\*(#H)/2u'\v'-.3n'\*(#[\z\(de\v'.3n'\h'|\\n:u'\*(#] .ds d- \h'\*(#H'\(pd\h'-\w'~'u'\v'-.25m'\f2\(hy\fP\v'.25m'\h'-\*(#H' .ds D- D\\k:\h'-\w'D'u'\v'-.11m'\z\(hy\v'.11m'\h'|\\n:u' .ds th \*(#[\v'.3m'\s+1I\s-1\v'-.3m'\h'-(\w'I'u*2/3)'\s-1o\s+1\*(#] .ds Th \*(#[\s+2I\s-2\h'-\w'I'u*3/5'\v'-.3m'o\v'.3m'\*(#] .ds ae a\h'-(\w'a'u*4/10)'e .ds Ae A\h'-(\w'A'u*4/10)'E . \" corrections for vroff .if v .ds ~ \\k:\h'-(\\n(.wu*9/10-\*(#H)'\s-2\u~\d\s+2\h'|\\n:u' .if v .ds ^ \\k:\h'-(\\n(.wu*10/11-\*(#H)'\v'-.4m'^\v'.4m'\h'|\\n:u' . \" for low resolution devices (crt and lpr) .if \n(.H>23 .if \n(.V>19 \ \{\ . ds : e . ds 8 ss . ds o a . ds d- d\h'-1'\(ga . ds D- D\h'-1'\(hy . ds th \o'bp' . ds Th \o'LP' . ds ae ae . ds Ae AE .\} .rm #[ #] #H #V #F C .\" ======================================================================== .\" .IX Title "Bio::DB::HTS 3pm" .TH Bio::DB::HTS 3pm "2020-11-09" "perl v5.32.0" "User Contributed Perl Documentation" .\" For nroff, turn off justification. Always turn off hyphenation; it makes .\" way too many mistakes in technical documents. .if n .ad l .nh .SH "LICENSE" .IX Header "LICENSE" Copyright [2015\-2018] EMBL-European Bioinformatics Institute .PP Licensed under the Apache License, Version 2.0 (the \*(L"License\*(R"); you may not use this file except in compliance with the License. You may obtain a copy of the License at .PP .Vb 1 \& http://www.apache.org/licenses/LICENSE\-2.0 .Ve .PP Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an \*(L"\s-1AS IS\*(R" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND,\s0 either express or implied. See the License for the specific language governing permissions and limitations under the License. .SH "NAME" Bio::DB::HTS \-\- Read files using HTSlib including BAM/CRAM, Tabix and BCF database files .SH "SYNOPSIS" .IX Header "SYNOPSIS" .Vb 1 \& use Bio::DB::HTS; \& \& # high level API \& # Note that the high level API does not reset the CRAM file pointer to the start \& # of the file as the method to do so is (at time or writing) not easily accessible. \& # Therefore a new HTS object may be needed to repeat a query. \& my $hts = Bio::DB::HTS\->new(\-bam =>"data/ex1.bam", \& \-fasta=>"data/ex1.fa", \& ); \& \& my @targets = $hts\->seq_ids; \& my @alignments = $hts\->get_features_by_location(\-seq_id => \*(Aqseq2\*(Aq, \& \-start => 500, \& \-end => 800); \& for my $a (@alignments) { \& \& # where does the alignment start in the reference sequence \& my $seqid = $a\->seq_id; \& my $start = $a\->start; \& my $end = $a\->end; \& my $strand = $a\->strand; \& my $cigar = $a\->cigar_str; \& my $paired = $a\->get_tag_values(\*(AqPAIRED\*(Aq); \& \& # where does the alignment start in the query sequence \& my $query_start = $a\->query\->start; \& my $query_end = $a\->query\->end; \& \& my $ref_dna = $a\->dna; # reference sequence bases \& my $query_dna = $a\->query\->dna; # query sequence bases \& \& my @scores = $a\->qscore; # per\-base quality scores \& my $match_qual= $a\->qual; # quality of the match \& } \& \& my @pairs = $hts\->get_features_by_location(\-type => \*(Aqread_pair\*(Aq, \& \-seq_id => \*(Aqseq2\*(Aq, \& \-start => 500, \& \-end => 800); \& \& for my $pair (@pairs) \& { \& my $length = $pair\->length; # insert length \& my ($first_mate,$second_mate) = $pair\->get_SeqFeatures; \& my $f_start = $first_mate\->start; \& my $s_start = $second_mate\->start; \& } \& \& # low level API \& my $hfile = Bio::DB::HTSfile\->open(\*(Aq/path/to/alignment_file\*(Aq); \& my $header = $hfile\->header_read; \& my $target_count = $header\->n_targets; \& my $target_names = $header\->target_name; \& while (my $align = $hfile\->read1($header)) \& { \& my $seqid = $target_names\->[$align\->tid]; \& my $start = $align\->pos+1; \& my $end = $align\->calend; \& my $cigar = $align\->cigar_str; \& } \& \& Bio::DB::HTSfile\->index_build($bamfile); \& my $index = Bio::DB::HTSfile\->index_load($hfile); \& my $index = Bio::DB::HTSfile\->index_open_in_safewd($hfile); \& \& my $callback = sub { \& my $alignment = shift; \& my $start = $alignment\->start; \& my $end = $alignment\->end; \& my $seqid = $target_names\->[$alignment\->tid]; \& print $alignment\->qname," aligns to $seqid:$start..$end\en"; \& } \& my $header = $index\->header; \& $index\->fetch($hfile,$header\->parse_region(\*(Aqseq2\*(Aq),$callback); .Ve .SH "DESCRIPTION" .IX Header "DESCRIPTION" This module provides a Perl interface to the HTSlib library for indexed and unindexed \s-1SAM/BAM/CRAM\s0 sequence alignment databases. It provides support for retrieving information on individual alignments, read pairs, and alignment coverage information across large regions. It also provides callback functionality for calling SNPs and performing other base-by-base functions. .SS "The high-level \s-1API\s0" .IX Subsection "The high-level API" The high-level \s-1API\s0 provides a BioPerl-compatible interface to indexed \&\s-1BAM\s0 and \s-1CRAM\s0 files. The alignment file database is treated as a collection of Bio::SeqFeatureI features, and can be searched for features by name, location, type and combinations of feature tags such as whether the alignment is part of a mate-pair. .PP When opening a alignment database using the high-level \s-1API,\s0 you provide the pathnames of two files: the \s-1FASTA\s0 file that contains the reference genome sequence, and the \s-1BAM\s0 file that contains the query sequences and their alignments. If either of the two files needs to be indexed, the indexing will need to be built. You can then query the database for alignment features by combinations of name, position, type, and feature tag. .PP The high-level \s-1API\s0 provides access to up to four feature \*(L"types\*(R": .PP .Vb 2 \& * "match": The "raw" unpaired alignment between a read and the \& reference sequence. \& \& * "read_pair": Paired alignments; a single composite \& feature that contains two subfeatures for the alignments of each \& of the mates in a mate pair. \& \& * "coverage": A feature that spans a region of interest that contains \& numeric information on the coverage of reads across the region. \& \& * "region": A way of retrieving information about the reference \& sequence. Searching for features of type "region" will return a \& list of chromosomes or contigs in the reference sequence, rather \& than read alignments. \& \& * "chromosome": A synonym for "region". .Ve .PP \&\fBFeatures\fR can be en masse in a single call, retrieved in a memory-efficient streaming basis using an iterator, or interrogated using a filehandle that return a series of SAM-format lines. .PP \&\fB\s-1SAM\s0 alignment flags\fR can be retrieved using BioPerl's feature \*(L"tag\*(R" mechanism. For example, to interrogate the \s-1FIRST_MATE\s0 flag, one fetches the \*(L"\s-1FIRST_MATE\*(R"\s0 tag: .PP .Vb 1 \& warn "aye aye captain!" if $alignment\->get_tag_values(\*(AqFIRST_MATE\*(Aq); .Ve .PP The Bio::SeqFeatureI interface has been extended to retrieve all flags as a compact human-readable string, and to return the \s-1CIGAR\s0 alignment in a variety of formats. .PP \&\fBSplit alignments\fR, such as reads that cover introns, are dealt with in one of two ways. The default is to leave split alignments alone: they can be detected by one or more \*(L"N\*(R" operations in the \s-1CIGAR\s0 string. Optionally, you can choose to have the \s-1API\s0 split these alignments across two or more subfeatures; the \s-1CIGAR\s0 strings of these split alignments will be adjusted accordingly. .PP \&\fBInterface to the pileup routines\fR The \s-1API\s0 provides you with access to the samtools \*(L"pileup\*(R" \s-1API.\s0 This gives you the ability to write a callback that will be invoked on every column of the alignment for the purpose of calculating coverage, quality score metrics, or \s-1SNP\s0 calling. .PP \&\fBAccess to the reference sequence\fR When you create the Bio::DB::HTS object, you can pass the path to a \s-1FASTA\s0 file containing the reference sequence. Alternatively, you may pass an object that knows how to retrieve \s-1DNA\s0 sequences across a range via the \fBseq()\fR or \fBfetch_seq()\fR methods, as described under \fBnew()\fR. .PP If the \s-1SAM/BAM\s0 file has \s-1MD\s0 tags, then these tags will be used to reconstruct the reference sequence when necessary, in which case you can completely omit the \-fasta argument. Note that not all \s-1SAM/BAM\s0 files have \s-1MD\s0 tags, and those that do may not use them correctly due to the newness of this part of the \s-1SAM\s0 spec. You may wish to populate these tags using samtools' \*(L"calmd\*(R" command. .PP If the \-fasta argument is omitted and no \s-1MD\s0 tags are present, then the reference sequence will be returned as 'N'. .PP The \fBmain object classes\fR that you will be dealing with in the high-level \s-1API\s0 are as follows: .PP .Vb 8 \& * Bio::DB::HTS \-\- A collection of alignments and reference sequences. \& * Bio::DB::HTS::Alignment \-\- The alignment between a query and the reference. \& * Bio::DB::HTS::Query \-\- An object corresponding to the query sequence in \& which both (+) and (\-) strand alignments are \& shown in the reference (+) strand. \& * Bio::DB::HTS::Target \-\- An interface to the query sequence in which \& (\-) strand alignments are shown in reverse \& complement .Ve .PP You may encounter other classes as well. These include: .PP .Vb 9 \& * Bio::DB::HTS::Segment \-\- This corresponds to a region on the reference \& sequence. \& * Bio::DB::HTS::Constants \-\- This defines CIGAR symbol constants and flags. \& * Bio::DB::HTS::AlignWrapper \-\- An alignment helper object that adds split \& alignment functionality. See Bio::DB::HTS::Alignment \& for the documentation on using it. \& * Bio::DB::HTS::ReadIterator \-\- An iterator that mediates the one\-feature\-at\-a\-time \& retrieval mechanism. \& * Bio::DB::HTS::FetchIterator \-\- Another iterator for feature\-at\-a\-time retrieval. .Ve .SS "The low-level \s-1API\s0" .IX Subsection "The low-level API" The low-level \s-1API\s0 closely mirrors that of the HTSlib library. It provides the ability to open and read \s-1SAM, BAM\s0 and \s-1CRAM\s0 files, build indexes, and perform searches across them. .PP The classes you will be interacting with in the low-level \s-1API\s0 are as follows: .PP .Vb 6 \& * Bio::DB::HTS \-\- Methods that read and write SAM, BAM and CRAM files. \& * Bio::DB::HTS::Header \-\- Methods for manipulating the BAM file header. \& * Bio::DB::HTS::Alignment \-\- Methods for manipulating alignment data. \& * Bio::DB::HTS::Pileup \-\- Methods for manipulating the pileup data structure. \& * Bio::DB::HTS::Fai \-\- Methods for creating and reading from indexed Fasta \& files. .Ve .SH "METHODS" .IX Header "METHODS" We cover the high-level \s-1API\s0 first. The high-level \s-1API\s0 code can be found in the files Bio/DB/HTS.pm and Bio/DB/HTS/*.pm. .SS "Bio::DB::HTS Constructor and basic accessors" .IX Subsection "Bio::DB::HTS Constructor and basic accessors" .ie n .IP "$sam = Bio::DB::HTS\->new(%options)" 4 .el .IP "\f(CW$sam\fR = Bio::DB::HTS\->new(%options)" 4 .IX Item "$sam = Bio::DB::HTS->new(%options)" The Bio::DB::HTS object combines a Fasta file of the reference sequences with an \s-1SAM/BAM/CRAM\s0 alignment file to allow for convenient retrieval of human-readable sequence IDs and reference sequences. The \fBnew()\fR constructor accepts a \-name=>value style list of options as follows: .Sp .Vb 2 \& Option Description \& \-\-\-\-\-\- \-\-\-\-\-\-\-\-\-\-\-\-\- \& \& \-bam Path to the SAM/BAM/CRAM alignment file that contains the \& alignments (required). A http: or ftp: URL is accepted. \& \& \-fasta Path to the Fasta file that contains \& the reference sequences (optional). Alternatively, \& you may pass any object that supports a seq() \& or fetch_seq() method and takes the three arguments \& ($seq_id,$start,$end). \& \& \-expand_flags A boolean value. If true then the standard \& alignment flags will be broken out as \& individual tags such as \*(AqM_UNMAPPED\*(Aq (default false). \& \& \-split_splices A boolean value. If true, then alignments that \& are split across splices will be broken out \& into a single alignment containing two sub\- \& alignments (default false). \& \& \-split The same as \-split_splices. \& \& \-force_refseq Always use the reference sequence file to derive the \& reference sequence, even when the sequence can be \& derived from the MD tag. This is slower, but safer \& when working with BAM files derived from buggy aligners \& or when the reference contains non\-canonical (modified) \& bases. \& \& \-autoindex Create an alignment index file if one does not exist \& or the current one has a modification date \& earlier than the alignment file. .Ve .Sp An example of a typical \fBnew()\fR constructor invocation is: .Sp .Vb 4 \& $hts = Bio::DB::HTS\->new(\-fasta => \*(Aq/home/projects/genomes/hu17.fa\*(Aq, \& \-bam => \*(Aq/home/projects/alignments/ej88.bam\*(Aq, \& \-expand_flags => 1, \& \-split_splices => 1); .Ve .Sp If the \fB\-fasta\fR argument is present, then you will be able to use the interface to fetch the reference sequence's bases. Otherwise, calls that return the reference sequence will return sequences consisting entirely of \*(L"N\*(R". .Sp \&\fB\-expand_flags\fR option, if true, has the effect of turning each of the standard \s-1SAM\s0 flags into a separately retrievable \fBtag\fR in the Bio::SeqFeatureI interface. Otherwise, the standard flags will be concatenated in easily parseable form as a tag named \*(L"\s-1FLAGS\*(R".\s0 See \&\fBget_all_tags()\fR and \fBget_tag_values()\fR for more information. .Sp Any two-letter extension flags, such as H0 or H1, will always appear as separate tags regardless of the setting. .Sp \&\fB\-split_splices\fR has the effect of breaking up alignments that contain an \*(L"N\*(R" operation into subparts for more convenient manipulation. For example, if you have both paired reads and spliced alignments in the \s-1BAM\s0 file, the following code shows the subpart relationships: .Sp .Vb 4 \& $pair = $hts\->get_feature_by_name(\*(AqE113:01:01:23\*(Aq); \& @mates = $pair\->get_SeqFeatures; \& @mate1_parts = $mates[0]\->get_SeqFeatures; \& @mate2_parts = $mates[1]\->get_SeqFeatures; .Ve .Sp Because there is some overhead to splitting up the spliced alignments, this option is false by default. .Sp \&\fBRemote access\fR to alignment files located on an \s-1HTTP\s0 or \s-1FTP\s0 server is possible. Simply replace the path to the \s-1BAM\s0 file with the appropriate \&\s-1URL.\s0 Note that incorrect URLs may lead to a core dump. .Sp It is not currently possible to refer to a remote \s-1FASTA\s0 file. These will have to be downloaded locally and indexed before using. .ie n .IP "$flag = $hts\->expand_flags([$new_value])" 4 .el .IP "\f(CW$flag\fR = \f(CW$hts\fR\->expand_flags([$new_value])" 4 .IX Item "$flag = $hts->expand_flags([$new_value])" Get or set the expand_flags option. This can be done after object creation and will have an immediate effect on all alignments fetched from the alignment file. .ie n .IP "$flag = $hts\->split_splices([$new_value])" 4 .el .IP "\f(CW$flag\fR = \f(CW$hts\fR\->split_splices([$new_value])" 4 .IX Item "$flag = $hts->split_splices([$new_value])" Get or set the split_splices option. This can be done after object creation and will affect all alignments fetched from the alignment file \&\fBsubsequently.\fR .ie n .IP "$header = $hts\->header" 4 .el .IP "\f(CW$header\fR = \f(CW$hts\fR\->header" 4 .IX Item "$header = $hts->header" Return the Bio::DB::HTS::Header object associated with the alignment file. You can manipulate the header using the low-level \s-1API.\s0 .ie n .IP "$hts_path = $hts\->hts_path" 4 .el .IP "\f(CW$hts_path\fR = \f(CW$hts\fR\->hts_path" 4 .IX Item "$hts_path = $hts->hts_path" Return the path of the alignment file used to create the hts object. This makes the object more portable. .ie n .IP "$hts_file = $hts\->$hts_file" 4 .el .IP "\f(CW$hts_file\fR = \f(CW$hts\fR\->$hts_file" 4 .IX Item "$hts_file = $hts->$hts_file" Returns the low-level Bio::DB::HTSfile object associated with the opened file. .ie n .IP "$fai = $hts\->fai" 4 .el .IP "\f(CW$fai\fR = \f(CW$hts\fR\->fai" 4 .IX Item "$fai = $hts->fai" Returns the Bio::DB::HTS::Fai object associated with the Fasta file. You can then manipulate this object with the low-level \s-1API.\s0 .Sp \&\fBThe index can be built automatically for you if it does not already exist.\fR If index building is necessarily, the process will need write privileges to the same directory in which the Fasta file resides.> If the process does not have write permission, then the call will fail. .ie n .IP "$hts_idx = $hts\->hts_index" 4 .el .IP "\f(CW$hts_idx\fR = \f(CW$hts\fR\->hts_index" 4 .IX Item "$hts_idx = $hts->hts_index" Return the Bio::DB::HTS::Index object associated with the alignment file. .Sp The index is not automatically built. .ie n .IP "$hts\->clone" 4 .el .IP "\f(CW$hts\fR\->clone" 4 .IX Item "$hts->clone" Bio::DB::HTS objects are not stable across \fBfork()\fR operations. If you fork, you must call \fBclone()\fR either in the parent or the child process before attempting to call any methods. .SS "Getting information about reference sequences" .IX Subsection "Getting information about reference sequences" The Bio::DB::HTS object provides the following methods for getting information about the reference sequence(s) contained in the associated Fasta file. .ie n .IP "@seq_ids = $hts\->seq_ids" 4 .el .IP "\f(CW@seq_ids\fR = \f(CW$hts\fR\->seq_ids" 4 .IX Item "@seq_ids = $hts->seq_ids" Returns an unsorted list of the IDs of the reference sequences (known elsewhere in this document as seq_ids). This is the same as the identifier following the \*(L">\*(R" sign in the Fasta file (e.g. \*(L"chr1\*(R"). .ie n .IP "$num_targets = $hts\->n_targets" 4 .el .IP "\f(CW$num_targets\fR = \f(CW$hts\fR\->n_targets" 4 .IX Item "$num_targets = $hts->n_targets" Return the number of reference sequences. .ie n .IP "$length = $hts\->length('seqid')" 4 .el .IP "\f(CW$length\fR = \f(CW$hts\fR\->length('seqid')" 4 .IX Item "$length = $hts->length('seqid')" Returns the length of the reference sequence named \*(L"seqid\*(R". .ie n .IP "$seq_id = $hts\->target_name($tid)" 4 .el .IP "\f(CW$seq_id\fR = \f(CW$hts\fR\->target_name($tid)" 4 .IX Item "$seq_id = $hts->target_name($tid)" Translates a numeric target \s-1ID\s0 (\s-1TID\s0) returned by the low-level \s-1API\s0 into a seq_id used by the high-level \s-1API.\s0 .ie n .IP "$length = $hts\->target_len($tid)" 4 .el .IP "\f(CW$length\fR = \f(CW$hts\fR\->target_len($tid)" 4 .IX Item "$length = $hts->target_len($tid)" Translates a numeric target \s-1ID\s0 (\s-1TID\s0) from the low-level \s-1API\s0 to a sequence length. .ie n .IP "$dna = $hts\->seq($seqid,$start,$end)" 4 .el .IP "\f(CW$dna\fR = \f(CW$hts\fR\->seq($seqid,$start,$end)" 4 .IX Item "$dna = $hts->seq($seqid,$start,$end)" Returns the \s-1DNA\s0 across the region from start to end on reference seqid. Note that this is a string, not a Bio::PrimarySeq object. If no \-fasta path was passed when the sam object was created, then you will receive a series of N nucleotides of the requested length. .SS "Creating and querying segments" .IX Subsection "Creating and querying segments" Bio::DB::HTS::Segment objects refer regions on the reference sequence. They can be used to retrieve the sequence of the reference, as well as alignments that overlap with the region. .ie n .IP "$segment = $hts\->segment($seqid,$start,$end);" 4 .el .IP "\f(CW$segment\fR = \f(CW$hts\fR\->segment($seqid,$start,$end);" 4 .IX Item "$segment = $hts->segment($seqid,$start,$end);" .PD 0 .ie n .IP "$segment = $hts\->segment(\-seq_id=>'chr1',\-start=>5000,\-end=>6000);" 4 .el .IP "\f(CW$segment\fR = \f(CW$hts\fR\->segment(\-seq_id=>'chr1',\-start=>5000,\-end=>6000);" 4 .IX Item "$segment = $hts->segment(-seq_id=>'chr1',-start=>5000,-end=>6000);" .PD Segments are created using the Bio:DB::HTS\->\fBsegment()\fR method. It can be called using one to three positional arguments corresponding to the seq_id of the reference sequence, and optionally the start and end positions of a subregion on the sequence. If the start and/or end are undefined, they will be replaced with the beginning and end of the sequence respectively. .Sp Alternatively, you may call \fBsegment()\fR with named \-seq_id, \-start and \&\-end arguments. .Sp All coordinates are 1\-based. .ie n .IP "$seqid = $segment\->seq_id" 4 .el .IP "\f(CW$seqid\fR = \f(CW$segment\fR\->seq_id" 4 .IX Item "$seqid = $segment->seq_id" Return the segment's sequence \s-1ID.\s0 .ie n .IP "$start = $segment\->start" 4 .el .IP "\f(CW$start\fR = \f(CW$segment\fR\->start" 4 .IX Item "$start = $segment->start" Return the segment's start position. .ie n .IP "$end = $segment\->end" 4 .el .IP "\f(CW$end\fR = \f(CW$segment\fR\->end" 4 .IX Item "$end = $segment->end" Return the segment's end position. .ie n .IP "$strand = $segment\->strand" 4 .el .IP "\f(CW$strand\fR = \f(CW$segment\fR\->strand" 4 .IX Item "$strand = $segment->strand" Return the strand of the segment (always 0). .ie n .IP "$length = $segment\->length" 4 .el .IP "\f(CW$length\fR = \f(CW$segment\fR\->length" 4 .IX Item "$length = $segment->length" Return the length of the segment. .ie n .IP "$dna = $segment\->dna" 4 .el .IP "\f(CW$dna\fR = \f(CW$segment\fR\->dna" 4 .IX Item "$dna = $segment->dna" Return the \s-1DNA\s0 string for the reference sequence under this segment. .ie n .IP "$seq = $segment\->seq" 4 .el .IP "\f(CW$seq\fR = \f(CW$segment\fR\->seq" 4 .IX Item "$seq = $segment->seq" Return a Bio::PrimarySeq object corresponding to the sequence of the reference under this segment. You can get the actual \s-1DNA\s0 string in this redundant-looking way: .Sp .Vb 1 \& $dna = $segment\->seq\->seq .Ve .Sp The advantage of working with a Bio::PrimarySeq object is that you can perform operations on it, including taking its reverse complement and subsequences. .ie n .IP "@alignments = $segment\->features(%args)" 4 .el .IP "\f(CW@alignments\fR = \f(CW$segment\fR\->features(%args)" 4 .IX Item "@alignments = $segment->features(%args)" Return alignments that overlap the segment in the associated alignment file. The optional \f(CW%args\fR list allows you to filter features by name, tag or other attributes. See the documentation of the Bio::DB::HTS\->\fBfeatures()\fR method for the full list of options. Here are some typical examples: .Sp .Vb 2 \& # get all the overlapping alignments \& @all_alignments = $segment\->features; \& \& # get an iterator across the alignments \& my $iterator = $segment\->features(\-iterator=>1); \& while (my $align = $iterator\->next_seq) { do something } \& \& # get a SAM filehandle across the alignments \& my $fh = $segment\->features(\-fh=>1); \& while (<$fh>) { print } \& \& # get only the alignments with unmapped mates \& my @unmapped = $segment\->features(\-flags=>{M_UNMAPPED=>1}); \& \& # get coverage across this region \& my ($coverage) = $segment\->features(\*(Aqcoverage\*(Aq); \& my @data_points = $coverage\->coverage; \& \& # grep through features using a coderef \& my @reverse_alignments = $segment\->features( \& \-filter => sub { \& my $a = shift; \& return $a\->strand < 0; \& }); .Ve .ie n .IP "$tag = $segment\->primary_tag" 4 .el .IP "\f(CW$tag\fR = \f(CW$segment\fR\->primary_tag" 4 .IX Item "$tag = $segment->primary_tag" .PD 0 .ie n .IP "$tag = $segment\->source_tag" 4 .el .IP "\f(CW$tag\fR = \f(CW$segment\fR\->source_tag" 4 .IX Item "$tag = $segment->source_tag" .PD Return the strings \*(L"region\*(R" and \*(L"sam/bam\*(R" respectively. These methods allow the segment to be passed to BioPerl methods that expect Bio::SeqFeatureI objects. .ie n .IP "$segment\->name, $segment\->display_name, $segment\->get_SeqFeatures, $segment\->get_tag_values" 4 .el .IP "\f(CW$segment\fR\->name, \f(CW$segment\fR\->display_name, \f(CW$segment\fR\->get_SeqFeatures, \f(CW$segment\fR\->get_tag_values" 4 .IX Item "$segment->name, $segment->display_name, $segment->get_SeqFeatures, $segment->get_tag_values" These methods are provided for Bio::SeqFeatureI compatibility and don't do anything of interest. .SS "Retrieving alignments, mate pairs and coverage information" .IX Subsection "Retrieving alignments, mate pairs and coverage information" The \fBfeatures()\fR method is an all-purpose tool for retrieving alignment information from the \s-1SAM/BAM/CRAM\s0 alignment file database. In addition, the methods \&\fBget_features_by_name()\fR, \fBget_features_by_location()\fR and others provide convenient shortcuts to \fBfeatures()\fR. .PP These methods either return a list of features, an iterator across a list of features, or a filehandle opened on a pseudo-SAM file. .ie n .IP "@features = $hts\->features(%options)" 4 .el .IP "\f(CW@features\fR = \f(CW$hts\fR\->features(%options)" 4 .IX Item "@features = $hts->features(%options)" .PD 0 .ie n .IP "$iterator = $hts\->features(\-iterator=>1,%more_options)" 4 .el .IP "\f(CW$iterator\fR = \f(CW$hts\fR\->features(\-iterator=>1,%more_options)" 4 .IX Item "$iterator = $hts->features(-iterator=>1,%more_options)" .ie n .IP "$filehandle = $hts\->features(\-fh=>1,%more_options)" 4 .el .IP "\f(CW$filehandle\fR = \f(CW$hts\fR\->features(\-fh=>1,%more_options)" 4 .IX Item "$filehandle = $hts->features(-fh=>1,%more_options)" .ie n .IP "@features = $hts\->features('type1','type2'...)" 4 .el .IP "\f(CW@features\fR = \f(CW$hts\fR\->features('type1','type2'...)" 4 .IX Item "@features = $hts->features('type1','type2'...)" .PD This is the all-purpose interface for fetching alignments and other types of features from the database. Arguments are a \-name=>value option list selected from the following list of options: .Sp .Vb 2 \& Option Description \& \-\-\-\-\-\- \-\-\-\-\-\-\-\-\-\-\-\-\- \& \& \-type Filter on features of a given type. You may provide \& either a scalar typename, or a reference to an \& array of desired feature types. Valid types are \& "match", "read_pair", "coverage" and "chromosome." \& \& See below for a full explanation of feature types. \& \& \-name Filter on reads with the designated name. Note that \& this can be a slow operation unless accompanied by \& the feature location as well. \& \& \-seq_id Filter on features that align to seq_id between start \& \-start and end. \-start and \-end must be used in conjunction \& \-end with \-seq_id. If \-start and/or \-end are absent, they \& will default to 1 and the end of the reference \& sequence, respectively. \& \& \-flags Filter features that match a list of one or more \& flags. See below for the format. \& \& \-attributes The same as \-flags, for compatibility with other \& \-tags APIs. \& \& \-filter Filter on features with a coderef. The coderef will \& receive a single argument consisting of the feature \& and should return true to keep the feature, or false \& to discard it. \& \& \-iterator Instead of returning a list of features, return an \& iterator across the results. To retrieve the results, \& call the iterator\*(Aqs next_seq() method repeatedly \& until it returns undef to indicate that no more \& matching features remain. \& \& \-fh Instead of returning a list of features, return a \& filehandle. Read from the filehandle to retrieve \& each of the results in TAM format, one alignment \& per line read. This only works for features of type \& "match." .Ve .Sp The high-level \s-1API\s0 introduces the concept of a \fBfeature \*(L"type\*(R"\fR in order to provide several convenience functions. You specify types by using the optional \fB\-type\fR argument. The following types are currently supported: .Sp \&\fBmatch\fR. The \*(L"match\*(R" type corresponds to the unprocessed \s-1SAM\s0 alignment. It will retrieve single reads, either mapped or unmapped. Each match feature's \fBprimary_tag()\fR method will return the string \*(L"match.\*(R" The features returned by this call are of type Bio::DB::HTS::AlignWrapper. .Sp \&\fBread_pair\fR. The \*(L"paired_end\*(R" type causes the sam interface to find and merge together mate pairs. Fetching this type of feature will yield a series of Bio::SeqFeatureI objects, each as long as the total distance on the reference sequence spanned by the mate pairs. The top-level feature is of type Bio::SeqFeature::Lite; it contains two Bio::DB::HTS::AlignWrapper subparts. .Sp Call \fBget_SeqFeatures()\fR to get the two individual reads. Example: .Sp .Vb 6 \& my @pairs = $hts\->features(\-type=>\*(Aqread_pair\*(Aq); \& my $p = $pairs[0]; \& my $i_length = $p\->length; \& my @ends = $p\->get_SeqFeatures; \& my $left = $ends[0]\->start; \& my $right = $ends[1]\->end; .Ve .Sp \&\fBcoverage\fR. The \*(L"coverage\*(R" type causes the sam interface to calculate coverage across the designated region. It only works properly if accompanied by the desired location of the coverage graph; \-seq_id is a mandatory argument for coverage calculation, and \-start and \-end are optional. The call will return a single Bio::SeqFeatureI object whose \&\fBprimary_tag()\fR is \*(L"coverage.\*(R" To recover the coverage data, call the object's \fBcoverage()\fR method to obtain an array (list context) or arrayref (scalar context) of coverage counts across the region of interest: .Sp .Vb 5 \& my ($coverage) = $hts\->features(\-type=>\*(Aqcoverage\*(Aq,\-seq_id=>\*(Aqseq1\*(Aq); \& my @data = $coverage\->coverage; \& my $total; \& for (@data) { $total += $_ } \& my $average_coverage = $total/@data; .Ve .Sp By default the coverage graph will be at the base pair level. So for a region 5000 bp wide, \fBcoverage()\fR will return an array or arrayref with exactly 5000 elements. However, you also have the option of calculating the coverage across larger bins. Simply append the number of intervals you are interested to the \*(L"coverage\*(R" typename. For example, fetching \*(L"coverage:500\*(R" will return a feature whose \&\fBcoverage()\fR method will return the coverage across 500 intervals. .Sp \&\fBchromosome\fR or \fBregion\fR. The \*(L"chromosome\*(R" or \*(L"region\*(R" type are interchangeable. They ask the sam interface to construct Bio::DB::HTS::Segment representing the reference sequences. These two calls give similar results: .Sp .Vb 3 \& my $segment = $hts\->segment(\*(Aqseq2\*(Aq,1=>500); \& my ($seg) = $hts\->features(\-type=>\*(Aqchromosome\*(Aq, \& \-seq_id=>\*(Aqseq2\*(Aq,\-start=>1,\-end=>500); .Ve .Sp Due to an unresolved bug, you cannot fetch chromosome features in the same call with matches and other feature types call. Specifically, this works as expected: .Sp .Vb 1 \& my @chromosomes = $hts\->features (\-type=>\*(Aqchromosome\*(Aq); .Ve .Sp But this doesn't (as of 18 June 2009): .Sp .Vb 1 \& my @chromosomes_and_matches = $hts\->features(\-type=>[\*(Aqmatch\*(Aq,\*(Aqchromosome\*(Aq]); .Ve .Sp If no \-type argument is provided, then \fBfeatures()\fR defaults to finding features of type \*(L"match.\*(R" .Sp You may call \fBfeatures()\fR with a plain list of strings (positional arguments, not \-type=>value arguments). This will be interpreted as a list of feature types to return: .Sp .Vb 1 \& my ($coverage) = $hts\->features(\*(Aqcoverage\*(Aq) .Ve .Sp For a description of the methods available in the features returned from this call, please see Bio::SeqfeatureI and Bio::DB::HTS::Alignment. .Sp You can \fBfilter\fR \*(L"match\*(R" and \*(L"read_pair\*(R" features by name, location and/or flags. The name and flag filters are not very efficient. Unless they are combined with a location filter, they will initiate an exhaustive search of the \s-1BAM\s0 database. .Sp Name filters are case-insensitive, and allow you to use shell-style \&\*(L"*\*(R" and \*(L"?\*(R" wildcards. Flag filters created with the \fB\-flag\fR, \&\fB\-attribute\fR or \fB\-tag\fR options have the following syntax: .Sp .Vb 4 \& \-flag => { FLAG_NAME_1 => [\*(Aqlist\*(Aq,\*(Aqof\*(Aq,\*(Aqpossible\*(Aq,\*(Aqvalues\*(Aq], \& FLAG_NAME_2 => [\*(Aqlist\*(Aq,\*(Aqof\*(Aq,\*(Aqpossible\*(Aq,\*(Aqvalues\*(Aq], \& ... \& } .Ve .Sp The value of \fB\-flag\fR is a hash reference in which the keys are flag names and the values are array references containing lists of acceptable values. The list of values are \s-1OR\s0'd with each other, and the flag names are \s-1AND\s0'd with each other. .Sp The \fB\-filter\fR option provides a completely generic filtering interface. Provide a reference to a subroutine. It will be called once for each potential feature. Return true to keep the feature, or false to discard it. Here is an example of how to find all matches whose alignment quality scores are greater than 80. .Sp .Vb 1 \& @features = $hts\->features(\-filter=>sub {shift\->qual > 80} ); .Ve .Sp By default, \fBfeatures()\fR returns a list of all matching features. You may instead request an iterator across the results list by passing \&\-iterator=>1. This will give you an object that has a single method, \&\fBnext_seq()\fR: .Sp .Vb 5 \& my $high_qual = $hts\->features(\-filter => sub {shift\->qual > 80}, \& \-iterator=> 1 ); \& while (my $feature = $high_qual\->next_seq) { \& # do something with the alignment \& } .Ve .Sp Similarly, by passing a true value to the argument \fB\-fh\fR, you can obtain a filehandle to a virtual \s-1SAM\s0 file. This only works with the \&\*(L"match\*(R" feature type: .Sp .Vb 6 \& my $high_qual = $hts\->features(\-filter => sub {shift\->qual > 80}, \& \-fh => 1 ); \& while (my $tam_line = <$high_qual>) { \& chomp($tam_line); \& # do something with it \& } .Ve .ie n .IP "@features = $hts\->get_features_by_name($name)" 4 .el .IP "\f(CW@features\fR = \f(CW$hts\fR\->get_features_by_name($name)" 4 .IX Item "@features = $hts->get_features_by_name($name)" Convenience method. The same as calling \f(CW$hts\fR\->features(\-name=>$name); .ie n .IP "$feature = $hts\->get_feature_by_name($name)" 4 .el .IP "\f(CW$feature\fR = \f(CW$hts\fR\->get_feature_by_name($name)" 4 .IX Item "$feature = $hts->get_feature_by_name($name)" Convenience method. The same as ($hts\->features(\-name=>$name))[0]. .ie n .IP "@features = $hts\->get_features_by_location($seqid,$start,$end)" 4 .el .IP "\f(CW@features\fR = \f(CW$hts\fR\->get_features_by_location($seqid,$start,$end)" 4 .IX Item "@features = $hts->get_features_by_location($seqid,$start,$end)" Convenience method. The same as calling \&\f(CW$hts\fR\->features(\-seq_id=>$seqid,\-start=>$start,\-end=>$end). .ie n .IP "@features = $hts\->get_features_by_flag(%flags)" 4 .el .IP "\f(CW@features\fR = \f(CW$hts\fR\->get_features_by_flag(%flags)" 4 .IX Item "@features = $hts->get_features_by_flag(%flags)" Convenience method. The same as calling \&\f(CW$hts\fR\->features(\-flags=>\e%flags). This method is also called \&\fBget_features_by_attribute()\fR and \fBget_features_by_tag()\fR. Example: .Sp .Vb 1 \& @features = $hts\->get_features_by_flag(H0=>1) .Ve .ie n .IP "$feature = $hts\->get_feature_by_id($id)" 4 .el .IP "\f(CW$feature\fR = \f(CW$hts\fR\->get_feature_by_id($id)" 4 .IX Item "$feature = $hts->get_feature_by_id($id)" The high-level \s-1API\s0 assigns each feature a unique \s-1ID\s0 composed of its read name, position and strand and returns it when you call the feature's \fBprimary_id()\fR method. Given that \s-1ID,\s0 this method returns the feature. .ie n .IP "$iterator = $hts\->get_seq_stream(%options)" 4 .el .IP "\f(CW$iterator\fR = \f(CW$hts\fR\->get_seq_stream(%options)" 4 .IX Item "$iterator = $hts->get_seq_stream(%options)" Convenience method. This is the same as calling \&\f(CW$hts\fR\->features(%options,\-iterator=>1). .ie n .IP "$fh = $hts\->get_seq_fh(%options)" 4 .el .IP "\f(CW$fh\fR = \f(CW$hts\fR\->get_seq_fh(%options)" 4 .IX Item "$fh = $hts->get_seq_fh(%options)" Convenience method. This is the same as calling \&\f(CW$hts\fR\->features(%options,\-fh=>1). .ie n .IP "$fh = $hts\->tam_fh" 4 .el .IP "\f(CW$fh\fR = \f(CW$hts\fR\->tam_fh" 4 .IX Item "$fh = $hts->tam_fh" Convenience method. It is the same as calling \f(CW$hts\fR\->features(\-fh=>1). .ie n .IP "@types = $hts\->types" 4 .el .IP "\f(CW@types\fR = \f(CW$hts\fR\->types" 4 .IX Item "@types = $hts->types" This method returns the list of feature types (e.g. \*(L"read_pair\*(R") returned by the current version of the interface. .SS "The generic \fBfetch()\fP and \fBpileup()\fP methods" .IX Subsection "The generic fetch() and pileup() methods" Lastly, the high-level \s-1API\s0 supports two methods for rapidly traversing indexed \s-1BAM\s0 databases. .ie n .IP "$hts\->fetch($region,$callback)" 4 .el .IP "\f(CW$hts\fR\->fetch($region,$callback)" 4 .IX Item "$hts->fetch($region,$callback)" This method traverses the indicated region and invokes a callback code reference on each match. Specify a region using the syntax \&\*(L"seqid:start\-end\*(R", or either of the alternative syntaxes \&\*(L"seqid:start..end\*(R" and \*(L"seqid:start,end\*(R". If start and end are absent, then the entire reference sequence is traversed. If end is absent, then the end of the reference sequence is assumed. .Sp The callback will be called repeatedly with a Bio::DB::HTS::AlignWrapper on the argument list. .Sp Example: .Sp .Vb 5 \& $hts\->fetch(\*(Aqseq1:600\-700\*(Aq, \& sub { \& my $a = shift; \& print $a\->display_name,\*(Aq \*(Aq,$a\->cigar_str,"\en"; \& }); .Ve .Sp Note that the \fBfetch()\fR operation works on reads that \fBoverlap\fR the indicated region. Therefore the callback may be called for reads that align to the reference at positions that start before or end after the indicated region. .ie n .IP "$hts\->pileup($region,$callback [,$keep_level])" 4 .el .IP "\f(CW$hts\fR\->pileup($region,$callback [,$keep_level])" 4 .IX Item "$hts->pileup($region,$callback [,$keep_level])" This method, which is named after the native \fBbam_lpileupfile()\fR function in the C interfaces, traverses the indicated region and generates a \*(L"pileup\*(R" of all the mapped reads that cover it. The user-provided callback function is then invoked on each position of the alignment along with a data structure that provides access to the individual aligned reads. .Sp As with \fBfetch()\fR, the region is specified as a string in the format \&\*(L"seqid:start\-end\*(R", \*(L"seqid:start..end\*(R" or \*(L"seqid:start,end\*(R". .Sp The callback is a coderef that will be invoked with three arguments: the seq_id of the reference sequence, the current position on the reference (in 1\-based coordinates!), and a reference to an array of Bio::DB::HTS::Pileup objects. Here is the typical call signature: .Sp .Vb 4 \& sub { \& my ($seqid,$pos,$pileup) = @_; \& # do something \& } .Ve .Sp For example, if you call pileup on the region \*(L"seq1:501\-600\*(R", then the callback will be invoked for all reads that overlap the indicated region. The first invocation of the callback will typically have a \&\f(CW$pos\fR argument somewhat to the left of the desired region and the last call will be somewhat to the right. You may wish to ignore positions that are outside of the requested region. Also be aware that the reference sequence position uses 1\-based coordinates, which is different from the low-level interface, which use 0\-based coordinates. .Sp The size of the \f(CW$pileup\fR array reference indicates the read coverage at that position. Here is a simple average coverage calculator: .Sp .Vb 10 \& my $depth = 0; \& my $positions = 0; \& my $callback = sub { \& my ($seqid,$pos,$pileup) = @_; \& next unless $pos >= 501 && $pos <= 600; \& $positions++; \& $depth += @$pileup; \& } \& $hts\->pileup(\*(Aqseq1:501\-600\*(Aq,$callback); \& print "coverage = ",$depth/$positions; .Ve .Sp Each Bio::DB::HTS::Pileup object describes the position of a read in the alignment. Briefly, Bio::DB::HTS::Pileup has the following methods: .Sp .Vb 2 \& $pileup\->alignment The alignment at this level (a \& Bio::DB::HTS::AlignWrapper object). \& \& $pileup\->qpos The position of the read base at the pileup site, \& in 0\-based coordinates. \& \& $pileup\->pos The position of the read base at the pileup site, \& in 1\-based coordinates; \& \& $pileup\->level The level of the read in the multiple alignment \& view. Note that this field is only valid when \& $keep_level is true, so it may not be relevant post \& htslib move. \& \& $pileup\->indel Length of the indel at this position: 0 for no indel, positive \& for an insertion (relative to the reference), negative for a \& deletion (relative to the reference.) \& \& $pileup\->is_del True if the base on the padded read is a deletion. \& \& $pileup\->is_refskip True if the base on the padded read is a gap relative to the reference (denoted as < or > in the pileup) \& \& $pileup\->is_head True if this is the first base in the query sequence. \& \& $pileup\->is_tail True if this is the last base in the query sequence. .Ve .Sp See \*(L"Examples\*(R" for a very simple \s-1SNP\s0 caller. .ie n .IP "$hts\->fast_pileup($region,$callback [,$keep_level])" 4 .el .IP "\f(CW$hts\fR\->fast_pileup($region,$callback [,$keep_level])" 4 .IX Item "$hts->fast_pileup($region,$callback [,$keep_level])" This is identical to \fBpileup()\fR except that the pileup object returns low-level Bio::DB::HTS::Alignment objects rather than the higher-level Bio::DB::HTS::AlignWrapper objects. This makes it roughly 50% faster, but you lose the align objects' \fBseq_id()\fR and \fBget_tag_values()\fR methods. As a compensation, the callback receives an additional argument corresponding to the Bio::DB::HTS object. You can use this to create AlignWrapper objects on an as needed basis: .Sp .Vb 8 \& my $callback = sub { \& my($seqid,$pos,$pileup,$hts) = @_; \& for my $p (@$pileup) { \& my $alignment = $p\->alignment; \& my $wrapper = Bio::DB::HTS::AlignWrapper\->new($alignment,$hts); \& my $has_mate = $wrapper\->get_tag_values(\*(AqPAIRED\*(Aq); \& } \& }; .Ve .IP "Bio::DB::HTS\->max_pileup_cnt([$new_cnt])" 4 .IX Item "Bio::DB::HTS->max_pileup_cnt([$new_cnt])" .PD 0 .ie n .IP "$hts\->max_pileup_cnt([$new_cnt])" 4 .el .IP "\f(CW$hts\fR\->max_pileup_cnt([$new_cnt])" 4 .IX Item "$hts->max_pileup_cnt([$new_cnt])" .PD The HTSlib library caps pileups at a set level, defaulting to 8000. The callback will not be invoked on a single position more than the level set by the cap, even if there are more reads. Called with no arguments, this method returns the current cap value. Called with a numeric argument, it changes the cap. There is currently no way to specify an unlimited cap. .Sp This method can be called as an instance method or a class method. .ie n .IP "$hts\->coverage2BedGraph([$fh])" 4 .el .IP "\f(CW$hts\fR\->coverage2BedGraph([$fh])" 4 .IX Item "$hts->coverage2BedGraph([$fh])" This special-purpose method will compute a four-column \s-1BED\s0 graph of the coverage across the entire alignment file and print it to \s-1STDOUT.\s0 You may provide a filehandle to redirect output to a file or pipe. .PP The next sections correspond to the low-level \s-1API,\s0 which let you create and manipulate Perl objects that correspond directly to data structures in the C interface. A major difference between the high and low level APIs is that in the high-level \s-1API,\s0 the reference sequence is identified using a human-readable seq_id. However, in the low-level \&\s-1API,\s0 the reference is identified using a numeric target \s-1ID\s0 (\*(L"tid\*(R"). The target \s-1ID\s0 is established during the creation of the alignment file and is a small 0\-based integer index. The Bio::DB::HTS::Header object provides methods for converting from seq_ids to tids. .SS "Indexed Fasta Files" .IX Subsection "Indexed Fasta Files" These methods relate to the indexed Fasta (\*(L".fai\*(R") files. .ie n .IP "$fai = Bio::DB::HTS::Fai\->load('/path/to/file.fa')" 4 .el .IP "\f(CW$fai\fR = Bio::DB::HTS::Fai\->load('/path/to/file.fa')" 4 .IX Item "$fai = Bio::DB::HTS::Fai->load('/path/to/file.fa')" Load an indexed Fasta file and return the object corresponding to it. If the index does not exist, it will be created automatically. Note that you pass the path to the Fasta file, not the index. .Sp For consistency with Bio::DB::HTS\->\fBopen()\fR this method is also called \&\fBopen()\fR. .ie n .IP "$dna_string = $fai\->fetch(""seqid:start\-end"")" 4 .el .IP "\f(CW$dna_string\fR = \f(CW$fai\fR\->fetch(``seqid:start\-end'')" 4 .IX Item "$dna_string = $fai->fetch(seqid:start-end)" Given a sequence \s-1ID\s0 contained in the Fasta file and optionally a subrange in the form \*(L"start-end\*(R", finds the indicated subsequence and returns it as a string. .SS "Alignment Files" .IX Subsection "Alignment Files" These methods provide interfaces to alignment files in \s-1SAM/BAM/CRAM\s0 format. .ie n .IP "$hts_file = Bio::DB::HTSfile\->open('/path/to/file.bam' [,$mode])" 4 .el .IP "\f(CW$hts_file\fR = Bio::DB::HTSfile\->open('/path/to/file.bam' [,$mode])" 4 .IX Item "$hts_file = Bio::DB::HTSfile->open('/path/to/file.bam' [,$mode])" Open the alignment file at the indicated path. Mode, if present, must be one of the file stream open flags (\*(L"r\*(R", \*(L"w\*(R", \*(L"wb\*(R", \*(L"wc\*(R", \*(L"a\*(R", \*(L"r+\*(R", etc.). If absent, mode defaults to \*(L"r\*(R". [write formats: w = \s-1SAM,\s0 wb = \s-1BAM,\s0 wc = \s-1CRAM\s0] .Sp Note that Bio::DB::HTS objects are not stable across \fBfork()\fR operations. If you fork, and intend to use the object in both parent and child, you must reopen the Bio::DB::HTS in either the child or the parent (but not both) before attempting to call any of the object's methods. .Sp The path may be an http: or ftp: \s-1URL,\s0 in which case a copy of the index file will be downloaded to the current working directory (see below) and all accesses will be performed on the remote \s-1BAM\s0 file. .Sp Example: .Sp .Vb 1 \& $hfile = Bio::DB::HTSfile\->open(\*(Aqhttp://some.site.com/nextgen/chr1_bowtie.bam\*(Aq); .Ve .ie n .IP "$header = $hfile\->\fBheader_read()\fR" 4 .el .IP "\f(CW$header\fR = \f(CW$hfile\fR\->\fBheader_read()\fR" 4 .IX Item "$header = $hfile->header_read()" Given an open alignment file, return a Bio::DB::HTS::Header object containing information about the reference sequence(s). Note that you must invoke \fBheader_read()\fR at least once before calling \fBread1()\fR. .ie n .IP "$status_code = $hfile\->header_write($header, [$reference])" 4 .el .IP "\f(CW$status_code\fR = \f(CW$hfile\fR\->header_write($header, [$reference])" 4 .IX Item "$status_code = $hfile->header_write($header, [$reference])" Given a Bio::DB::HTSfile::Header object and a \s-1BAM\s0 file opened in write mode, write the header to the file. If the write fails the process will be terminated at the C layer. If \f(CW$hfile\fR is \s-1CRAM\s0 formated a second argument \f(CW$reference\fR, which is the path to the reference Fasta file, must be passed. The result code is (currently) always zero. .ie n .IP "$alignment = $hfile\->read1($header)" 4 .el .IP "\f(CW$alignment\fR = \f(CW$hfile\fR\->read1($header)" 4 .IX Item "$alignment = $hfile->read1($header)" Read one alignment from the alignment file and return it as a Bio::DB::HTS::Alignment object. The \f(CW$header\fR parameter is returned by invoking \fBheader()\fR. .ie n .IP "$bytes = $hfile\->write1($header, $alignment)" 4 .el .IP "\f(CW$bytes\fR = \f(CW$hfile\fR\->write1($header, \f(CW$alignment\fR)" 4 .IX Item "$bytes = $hfile->write1($header, $alignment)" Given a \s-1BAM\s0 file that has been opened in write mode and a Bio::DB::HTS::Alignment object, write the alignment to the \s-1BAM\s0 file and return the number of bytes successfully written. .SS "Index methods" .IX Subsection "Index methods" The Bio::DB::HTS::Index object provides access to index (.bai|.csi, .crai) files. .ie n .IP "$status_code = Bio::DB::HTS\->index_build('/path/to/file.?am')" 4 .el .IP "\f(CW$status_code\fR = Bio::DB::HTS\->index_build('/path/to/file.?am')" 4 .IX Item "$status_code = Bio::DB::HTS->index_build('/path/to/file.?am')" Given the path to an alignment file, this function attempts to build an index. The process in which the alignment file exists must be writable by the current process and there must be sufficient disk space for the operation or the process will be terminated in the C library layer. The result code is currently always zero, but in the future may return a negative value to indicate failure. .Sp The index file built will depend on the alignment file type specified. For \s-1CRAM\s0 this will be a .crai file, for \s-1BAM\s0 .bai. .ie n .IP "$index = Bio::DB::HTS\->index('/path/to/file.?am',$reindex)" 4 .el .IP "\f(CW$index\fR = Bio::DB::HTS\->index('/path/to/file.?am',$reindex)" 4 .IX Item "$index = Bio::DB::HTS->index('/path/to/file.?am',$reindex)" Attempt to open the index for the indicated alignment file. If \f(CW$reindex\fR is true, and the index either does not exist or is out of date with respect to the alignment file (by checking modification dates), then attempt to rebuild the index. Will throw an exception if the index does not exist or if attempting to rebuild the index was unsuccessful. .ie n .IP "$index = Bio::DB::HTS\->index_load('/path/to/file.?am')" 4 .el .IP "\f(CW$index\fR = Bio::DB::HTS\->index_load('/path/to/file.?am')" 4 .IX Item "$index = Bio::DB::HTS->index_load('/path/to/file.?am')" Attempt to open the index file for an alignment file, returning a Bio::DB::HTS::Index object. The filename path to use is the alignment file, not the index file (i.e. .bam or .cram, not .bai|.csi or .crai) .ie n .IP "$index = Bio::DB::HTS\->index_open_in_safewd('/path/to/file.?am' [,$mode])" 4 .el .IP "\f(CW$index\fR = Bio::DB::HTS\->index_open_in_safewd('/path/to/file.?am' [,$mode])" 4 .IX Item "$index = Bio::DB::HTS->index_open_in_safewd('/path/to/file.?am' [,$mode])" When opening a remote alignmentfile, you may not wish for the index to be downloaded to the current working directory. This version of index_open copies the index into the directory indicated by the \s-1TMPDIR\s0 environment variable or the system-defined /tmp directory if not present. You may change the environment variable just before the call to change its behavior. .ie n .IP "$code = $index\->fetch($hfile,$tid,$start,$end,$callback [,$callback_data])" 4 .el .IP "\f(CW$code\fR = \f(CW$index\fR\->fetch($hfile,$tid,$start,$end,$callback [,$callback_data])" 4 .IX Item "$code = $index->fetch($hfile,$tid,$start,$end,$callback [,$callback_data])" This is the low-level equivalent of the \f(CW$hts\fR\->\fBfetch()\fR function described for the high-level \s-1API.\s0 Given a open \s-1BAM\s0 file object, the numeric \s-1ID\s0 of the reference sequence, start and end ranges on the reference, and a coderef, this function will traverse the region and repeatedly invoke the coderef with each Bio::DB::HTS::Alignment object that overlaps the region. .Sp Arguments: .Sp .Vb 2 \& Argument Description \& \-\-\-\-\-\-\-\- \-\-\-\-\-\-\-\-\-\-\- \& \& $hts_file The Bio::DB::HTSfile object that corresponds to the \& index object. \& \& $tid The target ID of the reference sequence. This can \& be obtained by calling $header\->parse_region() with \& an appropriate opened Bio::DB::HTS::Header object. \& \& $start The start and end positions of the desired range on \& the reference sequence given by $tid, in 0\-based \& $end coordinates. Like the $tid, these can be obtained from \& $header\->parse_region(). \& \& $callback A coderef that will be called for each read overlapping \& the designated region. \& \& $callback_data Any arbitrary Perl data that you wish to pass to the \& $callback (optional). .Ve .Sp The coderef's call signature should look like this: .Sp .Vb 4 \& my $callback = sub { \& my ($alignment,$data) = @_; \& ... \& } .Ve .Sp The first argument is a Bio::DB::HTS::Alignment object. The second is the callback data (if any) passed to \fBfetch()\fR. .Sp \&\fBFetch()\fR returns an integer code, but its meaning is not described in the \s-1SAM/BAM C\s0 library documentation. .ie n .IP "$index\->pileup($htsfile,$tid,$start,$end,$callback [,$callback_data])" 4 .el .IP "\f(CW$index\fR\->pileup($htsfile,$tid,$start,$end,$callback [,$callback_data])" 4 .IX Item "$index->pileup($htsfile,$tid,$start,$end,$callback [,$callback_data])" This is the low-level version of the \fBpileup()\fR method, which allows you to invoke a coderef for every position in a \s-1BAM\s0 alignment. Arguments are: .Sp .Vb 2 \& Argument Description \& \-\-\-\-\-\-\-\- \-\-\-\-\-\-\-\-\-\-\- \& \& $hts_file The Bio::DB::HTSfile object that corresponds to the \& index object. \& \& $tid The target ID of the reference sequence. This can \& be obtained by calling $header\->parse_region() with \& an appropriate opened Bio::DB::HTS::Header object. \& \& $start The start and end positions of the desired range on \& the reference sequence given by $tid, in 0\-based \& $end coordinates. Like the $tid, these can be obtained from \& $header\->parse_region(). \& \& $callback A coderef that will be called for each position of the \& alignment across the designated region. \& \& $callback_data Any arbitrary Perl data that you wish to pass to the \& $callback (optional). .Ve .Sp The callback will be invoked with four arguments corresponding to the numeric sequence \s-1ID\s0 of the reference sequence, the \fBzero-based\fR position on the alignment, an arrayref of Bio::DB::HTS::Pileup objects, and the callback data, if any. A typical call signature will be this: .Sp .Vb 5 \& $callback = sub { \& my ($tid,$pos,$pileups,$callback_data) = @_; \& for my $pileup (@$pileups) { \& # do something \& }; .Ve .Sp Note that the position argument is zero-based rather than 1\-based, as it is in the high-level \s-1API.\s0 .Sp The Bio::DB::HTS::Pileup object was described earlier in the description of the high-level \fBpileup()\fR method. .ie n .IP "$coverage = $index\->coverage($hfile,$tid,$start,$end [,$bins [,maxcnt]])" 4 .el .IP "\f(CW$coverage\fR = \f(CW$index\fR\->coverage($hfile,$tid,$start,$end [,$bins [,maxcnt]])" 4 .IX Item "$coverage = $index->coverage($hfile,$tid,$start,$end [,$bins [,maxcnt]])" Calculate coverage for the region on the target sequence given by \f(CW$tid\fR between positions \f(CW$start\fR and \f(CW$end\fR (zero-based coordinates). This method will return an array reference equal to the size of the region (by default). Each element of the array will be an integer indicating the number of reads aligning over that position. If you provide an option binsize in \f(CW$bins\fR, the array will be \f(CW$bins\fR elements in length, and each element will contain the average coverage over that region as a floating point number. .Sp By default, the underlying Samtools library caps coverage counting at a fixed value of 8000. You may change this default by providing an optional numeric sixth value, which changes the cap for the duration of the call, or by invoking Bio::DB::HTS\->max_pileup_cnt($new_value), which changes the cap permanently. Unfortunately there is no way of specifying that you want an unlimited cap. .SS "\s-1HTS\s0 header methods" .IX Subsection "HTS header methods" The Bio::DB::HTS::Header object contains information regarding the reference sequence(s) used to construct the corresponding alignment file. It is most frequently used to translate between numeric target IDs and human-readable seq_ids. Headers can be created by reading from a \s-1BAM\s0 file using Bio::DB::HTS\->\fBheader()\fR. You can also create header objects from scratch, although there is not much that you can do with such objects at this point. .ie n .IP "$header = Bio::DB::HTS::Header\->\fBnew()\fR" 4 .el .IP "\f(CW$header\fR = Bio::DB::HTS::Header\->\fBnew()\fR" 4 .IX Item "$header = Bio::DB::HTS::Header->new()" Return a new, empty, header object. .ie n .IP "$n_targets = $header\->n_targets" 4 .el .IP "\f(CW$n_targets\fR = \f(CW$header\fR\->n_targets" 4 .IX Item "$n_targets = $header->n_targets" Return the number of reference sequences in the database. .ie n .IP "$name_arrayref = $header\->target_name" 4 .el .IP "\f(CW$name_arrayref\fR = \f(CW$header\fR\->target_name" 4 .IX Item "$name_arrayref = $header->target_name" Return a reference to an array of reference sequence names, corresponding to the high-level \s-1API\s0's seq_ids. .Sp To convert from a target \s-1ID\s0 to a seq_id, simply index into this array: .Sp .Vb 1 \& $seq_id = $header\->target_name\->[$tid]; .Ve .ie n .IP "$length_arrayref = $header\->target_len" 4 .el .IP "\f(CW$length_arrayref\fR = \f(CW$header\fR\->target_len" 4 .IX Item "$length_arrayref = $header->target_len" Return a reference to an array of reference sequence lengths. To get the length of the sequence corresponding to \f(CW$tid\fR, just index into the array returned by \fBtarget_len()\fR: .Sp .Vb 1 \& $length = $header\->target_len\->[$tid]; .Ve .ie n .IP "$text = $header\->text" 4 .el .IP "\f(CW$text\fR = \f(CW$header\fR\->text" 4 .IX Item "$text = $header->text" .PD 0 .ie n .IP "$header\->text(""new value"")" 4 .el .IP "\f(CW$header\fR\->text(``new value'')" 4 .IX Item "$header->text(new value)" .PD Read the text portion of the header. The text can be replaced by providing the replacement string as an argument. Note that you should follow the header conventions when replacing the header text. No parsing or other error-checking is performed. .ie n .IP "($tid,$start,$end) = $header\->parse_region(""seq_id:start\-end"")" 4 .el .IP "($tid,$start,$end) = \f(CW$header\fR\->parse_region(``seq_id:start\-end'')" 4 .IX Item "($tid,$start,$end) = $header->parse_region(seq_id:start-end)" Given a string in the format \*(L"seqid:start\-end\*(R" (using a human-readable seq_id and 1\-based start and end coordinates), parse the string and return the target \s-1ID\s0 and start and end positions in 0\-based coordinates. If the range is omitted, then the start and end coordinates of the entire sequence is returned. If only the end position is omitted, then the end of the sequence is assumed. .ie n .IP "$header\->view1($alignment)" 4 .el .IP "\f(CW$header\fR\->view1($alignment)" 4 .IX Item "$header->view1($alignment)" This method will accept a Bio::DB::HTS::Alignment object, convert it to a line of \s-1TAM\s0 output, and write the output to \s-1STDOUT.\s0 In the low-level \s-1API\s0 there is currently no way to send the output to a different filehandle or capture it as a string. .SS "Bio::DB::HTS::Pileup methods" .IX Subsection "Bio::DB::HTS::Pileup methods" An array of Bio::DB::HTS::Pileup object is passed to the \fBpileup()\fR callback for each position of a multi-read alignment. Each pileup object contains information about the alignment of a single read at a single position. .ie n .IP "$alignment = $pileup\->alignment" 4 .el .IP "\f(CW$alignment\fR = \f(CW$pileup\fR\->alignment" 4 .IX Item "$alignment = $pileup->alignment" Return the Bio::DB::HTS::Alignment object at this level. This provides you with access to the aligning read. .ie n .IP "$alignment = $pileup\->b" 4 .el .IP "\f(CW$alignment\fR = \f(CW$pileup\fR\->b" 4 .IX Item "$alignment = $pileup->b" An alias for \fBalignment()\fR, provided for compatibility with the C \s-1API.\s0 .ie n .IP "$pos = $pileup\->qpos" 4 .el .IP "\f(CW$pos\fR = \f(CW$pileup\fR\->qpos" 4 .IX Item "$pos = $pileup->qpos" The position of the aligning base in the read in zero-based coordinates. .ie n .IP "$pos = $pileup\->pos" 4 .el .IP "\f(CW$pos\fR = \f(CW$pileup\fR\->pos" 4 .IX Item "$pos = $pileup->pos" The position of the aligning base in 1\-based coordinates. .ie n .IP "$level = $pileup\->level" 4 .el .IP "\f(CW$level\fR = \f(CW$pileup\fR\->level" 4 .IX Item "$level = $pileup->level" The \*(L"level\*(R" of the read in the BAM-generated text display of the alignment. .ie n .IP "$indel = $pileup\->indel" 4 .el .IP "\f(CW$indel\fR = \f(CW$pileup\fR\->indel" 4 .IX Item "$indel = $pileup->indel" Length of the indel at this position: 0 for no indel, positive for an insertion (relative to the reference), negative for a deletion (relative to the reference sequence.) .ie n .IP "$flag = $pileup\->is_del" 4 .el .IP "\f(CW$flag\fR = \f(CW$pileup\fR\->is_del" 4 .IX Item "$flag = $pileup->is_del" True if the base on the padded read is a deletion. .ie n .IP "$flag = $pileup\->is_refskip" 4 .el .IP "\f(CW$flag\fR = \f(CW$pileup\fR\->is_refskip" 4 .IX Item "$flag = $pileup->is_refskip" True if the base on the padded read is a gap relative to the reference (denoted as < or > in the pileup) .ie n .IP "$flag = $pileup\->is_head" 4 .el .IP "\f(CW$flag\fR = \f(CW$pileup\fR\->is_head" 4 .IX Item "$flag = $pileup->is_head" True if this is the first base in the query sequence. .ie n .IP "$flag = $pileup\->is_tail" 4 .el .IP "\f(CW$flag\fR = \f(CW$pileup\fR\->is_tail" 4 .IX Item "$flag = $pileup->is_tail" True if this is the last base in the query sequence. .SS "The alignment objects" .IX Subsection "The alignment objects" Please see Bio::DB::HTS::Alignment for documentation of the Bio::DB::HTS::Alignment and Bio::DB::HTS::AlignWrapper objects. .SH "DEPENDENCIES" .IX Header "DEPENDENCIES" Module::Build, Carp, Bio::Perl (>=1.006001), Test::More .SH "EXPORT" .IX Header "EXPORT" None .SH "AUTHORS" .IX Header "AUTHORS" Rishi Nag , original author. .PP Alessandro Vullo \f(CW\*(C`\*(C'\fR, the current developer and maintainer. .SH "CONTRIBUTORS" .IX Header "CONTRIBUTORS" Andy Yates, Keiran Raine, John Marshall, Zhicheng Liu, Can Wood, Dietmar Rieder, Chris Fields, David Jones, James Gilbert, Alex Hodgkins (Congenica Ltd.), Rob Aganrab .SH "KNOWN BUGS" .IX Header "KNOWN BUGS" .IP "\(bu" 4 \&\s-1SAM\s0 file reading and iterating over alignments does not work with older htslib versions (<1.5) .IP "\(bu" 4 The \fBpadded_alignment()\fR function with \s-1CRAM\s0 files may produce invalid output: unequal lenght of the strings that specify the pairwise alignment .PP Please report any bugs or feature requests to \f(CW\*(C`bug\-bio\-db\-hts at rt.cpan.org\*(C'\fR, or through the web interface at . I will be notified, and then you'll automatically be notified of progress on your bug as I make changes. .SH "TESTING AND CONTRIBUTING" .IX Header "TESTING AND CONTRIBUTING" You can obtain the most recent development version of this module via the GitHub repository at https://github.com/Ensembl/Bio\-DB\-HTS. Please feel free to submit bug reports, patches etc. .SH "SUPPORT" .IX Header "SUPPORT" You can find documentation for this module with the perldoc command. .PP .Vb 1 \& perldoc Bio::DB::HTS .Ve .PP You can also look for information at: .IP "\(bu" 4 \&\s-1RT: CPAN\s0's request tracker (report bugs here) .Sp .IP "\(bu" 4 AnnoCPAN: Annotated \s-1CPAN\s0 documentation .Sp .IP "\(bu" 4 \&\s-1CPAN\s0 Ratings .Sp .IP "\(bu" 4 Search \s-1CPAN\s0 .Sp .SH "EXAMPLES" .IX Header "EXAMPLES" For illustrative purposes only, here is an extremely stupid \s-1SNP\s0 caller that tallies up bases that are q>20 and calls a \s-1SNP\s0 if there are at least 4 non\-N/non\-indel bases at the position and at least 25% of them are a non-reference base. .PP .Vb 8 \& my @SNPs; # this will be list of SNPs \& my $snp_caller = sub { \& my ($seqid,$pos,$p) = @_; \& my $refbase = $hts\->segment($seqid,$pos,$pos)\->dna; \& my ($total,$different); \& for my $pileup (@$p) { \& my $b = $pileup\->alignment; \& next if $pileup\->indel or $pileup\->is_refskip; # don\*(Aqt deal with these ;\-) \& \& my $qbase = substr($b\->qseq,$pileup\->qpos,1); \& next if $qbase =~ /[nN]/; \& \& my $qscore = $b\->qscore\->[$pileup\->qpos]; \& next unless $qscore > 25; \& \& $total++; \& $different++ if $refbase ne $qbase; \& } \& if ($total >= 4 && $different/$total >= 0.25) { \& push @SNPs,"$seqid:$pos"; \& } \& }; \& \& $hts\->pileup(\*(Aqseq1\*(Aq,$snp_caller); \& print "Found SNPs: @SNPs\en"; .Ve .SH "GBrowse Compatibility" .IX Header "GBrowse Compatibility" The Bio::DB::HTS interface can be used as a backend to GBrowse (gmod.sourceforge.net/gbrowse). GBrowse can calculate and display coverage graphs across large regions, alignment cartoons across intermediate size regions, and detailed base-pair level alignments across small regions. .PP Here is a typical configuration for a \s-1BAM\s0 database that contains information from a shotgun genomic sequencing project. Some notes: .PP .Vb 3 \& * It is important to set "search options = none" in order to avoid \& GBrowse trying to scan through the BAM database to match read \& names. This is a time\-consuming operation. \& \& * The callback to "bgcolor" renders pairs whose mates are unmapped in \& red. \& \& * The callback to "balloon hover" causes a balloon to pop up with the \& read name when the user hovers over each paired read. Otherwise the \& default behavior would be to provide information about the pair as \& a whole. \& \& * When the user zooms out to 1001 bp or greaterp, the track switches \& to a coverage graph. \& \& [bamtest:database] \& db_adaptor = Bio::DB::HTSfile \& db_args = \-bam /var/www/gbrowse2/databases/bamtest/ex1.bam \& search options= default \& \& [Pair] \& feature = read_pair \& glyph = segments \& database = bamtest \& draw_target = 1 \& show_mismatch = 1 \& bgcolor = sub { \& my $f = shift; \& return $f\->get_tag_values(\*(AqM_UNMAPPED\*(Aq) ? \*(Aqred\*(Aq : \*(Aqgreen\*(Aq; \& } \& fgcolor = green \& height = 3 \& label = sub {shift\->display_name} \& label density = 50 \& bump = fast \& connector = dashed \& balloon hover = sub { \& my $f = shift; \& return \*(Aq\*(Aq unless $f\->type eq \*(Aqmatch\*(Aq; \& return \*(AqRead: \*(Aq.$f\->display_name.\*(Aq : \*(Aq.$f\->flag_str; \& } \& key = Read Pairs \& \& [Pair:1000] \& feature = coverage:1001 \& glyph = wiggle_xyplot \& height = 50 \& min_score = 0 \& autoscale = local .Ve .PP To show alignment data correctly when the user is zoomed in, you should also provide a pointer to the \s-1FASTA\s0 file containing the reference genome. In this case, modify the db_args line to read: .PP .Vb 2 \& db_args = \-bam /var/www/gbrowse2/databases/bamtest/ex1.bam \& \-fasta /var/www/gbrowse2/databases/bamtest/ex1.fa .Ve .SH "SEE ALSO" .IX Header "SEE ALSO" Bio::Perl, Bio::DB::HTS::Alignment, Bio::DB::HTS::Constants