Content Comparison

The BioITeam maintains a set of reference indexes for many common organisms and aligners. They can be found in aligner-specific sub-directories of the /work/projects/BioITeam/ref_genome area. E.g.:

/work/projects/BioITeam/ref_genome/
   bowtie2/
   bwa/
   hisat2/
   kallisto/
   star/
   tophat/

cds
mkdir -p core_ngs/references/fasta
cd core_ngs/references/fasta
cp $CORENGS/references/fasta/*.fa .

idev -m 180120 -N 1 -A OTH21164 -r core-ngs-class-0605  idev# Thursday
idev -m 120 -N 1 -A OTH21164 -r core-ngs-class-0606  # Friday

# or, without the reservation
idev -m 120 -N 1 -A OTH21164 -p development  

mkdir -p $SCRATCH/core_ngs/references/bwa/sacCer3
cd $SCRATCH/core_ngs/references/bwa/sacCer3

# Create a symbolic link to the sacCer3 FASTA file
# so it can
# be accessed directly from this directory
ln -sf ../../fasta/sacCer3.fa
ls -l

bwa aln -t 60 sacCer3/sacCer3.fa \
  fastq/Sample_Yeast_L005_R2.cat.fastq.gz > yeast_pe_R2.sai

The SAM file must contain both mapped and un-mapped reads, because there were 1,184,360 R1+R2 reads and the same number of alignment records.

Alignment records occur in the same read-name order as they did in the FASTQ, except that they come in pairs. The R1 read comes 1st, then the corresponding R2. This is called read name ordering.

The command below uses cut to extract the 1st Tabtab-delimited field (-f 1) from the first few alignment records:

grep -Pv '^@' yeast_pe.sam | head | cut -f 1

Next we'll see how to tell which read is the R1 and which is the R2.

Program: samtools (Tools for alignments in the SAM format)
Version: 1.920 (using htslib 1.920)

Usage:   samtools <command> [options]

Commands:
 
-- Indexing
 
   dict           create a sequence dictionary file
 
   faidx          index/extract FASTA
 
   fqidx          index/extract FASTQ
 
   index          index alignment



-- Editing
 
   calmd          recalculate MD/NM tags and '=' bases
 
   fixmate        fix mate information
 
   reheader       replace BAM header
 
   targetcut      cut fosmid regions (for fosmid pool only)
 
   addreplacerg   adds or replaces RG tags
 
   markdup        mark duplicates
   ampliconclip   clip oligos from the end of reads

-- File operations

    collate        shuffle and group alignments by name
     cat            concatenate BAMs
   consensus      produce a consensus Pileup/FASTA/FASTQ
   merge          merge sorted alignments
 
   mpileup        multi-way pileup
 
   sort           sort alignment file
 
   split          splits a file by read group
 
   quickcheck     quickly check if SAM/BAM/CRAM file appears intact
     fastq          converts a BAM to a FASTQ
 
   fasta          converts a BAM to a FASTA
   --import Statistics      bedcov  Converts FASTA or FASTQ files to SAM/BAM/CRAM
read depth per BEDreference region     Generates coveragea reference from aligned data
  alignment depthreset and percent coverage      depth Reverts aligner changes in reads

-- Statistics
   bedcov         read depth per BED region
   coverage       alignment depth and percent coverage
   depth          compute the depth

    flagstat       simple stats
     idxstats       BAM index stats
   cram-size      list CRAM Content-ID and Data-Series sizes
   phase          phase heterozygotes
 
   stats          generate stats (former bamcheck)
   ampliconstats  generate amplicon specific stats

-- Viewing
   flags explain flags BAM flags
   head    explain BAM flags     header tviewviewer
   tview  text    text alignment viewer
     view           SAM<->BAM<->CRAM conversion

    depad          convert padded BAM to unpadded BAM
   samples        list the samples in a set of SAM/BAM/CRAM files

-- Misc
   help [cmd]     display this help message or help for [cmd]
   version        detailed version information

There are two main "eras" of SAMtools development:

• "original" samtools, controlled by author Heng Li
  • v 0.1.19 is the last stable version
• "modern" samtools, controlled by a group of programmers
  • v 1.0, 1.1, 1.2 – avoid these (very buggy!)
  • v 1.3+ – finally stable!

Unfortunately, some functions with the same name in both version eras have different options and arguments! So be sure you know which version you're using. (The samtools version is usually reported at the top of its usage listing).

TACC BioContainers also offers the original samtools version: samtools/ctr-0.1.19--3.

Usage: samtools sort [options...] [in.bam]
Options:
  -l INT     Set compression level, from 0 (uncompressed) to 9 (best)
  -mu INT     Set maximum memory perOutput thread;uncompressed suffix K/M/G data (equivalent to -l 0)
  -m INT     Set maximum memory per thread; suffix K/M/G recognized [768M]
  -n         Sort by read nameM         Use minimiser for clustering unaligned/unplaced reads
  -R         Do not use reverse strand (only compatible with -M)
  -K INT     Kmer size to use for minimiser [20]
  -I FILE    Order minimisers by their position in FILE FASTA
  -w INT     Window size for minimiser indexing via -I ref.fa [100]
  -H         Squash homopolymers when computing minimiser
  -n         Sort by read name (natural): cannot be used with samtools index
  -N         Sort by read name (ASCII): cannot be used with samtools index
  -t TAG     Sort by value of TAG. Uses position as secondary index (or read name if -n is set)
  -o FILE    Write final output to FILE rather than standard output
  -T PREFIX  Write temporary files to PREFIX.nnnn.bam
      --no-PG
               Do not add a PG line
      --template-coordinate
               Sort by template-coordinate
      --input-fmt-option OPT[=VAL]
               Specify a single input file format option in the form
              
of OPTION or OPTION=VALUE
  -O, --output-fmt FORMAT[,OPT[=VAL]]...
               Specify output format (SAM, BAM, CRAM)
      --output-fmt-option OPT[=VAL]
     
         Specify a single output file format option in the form      Specify a single output file format option in the form
of OPTION or OPTION=VALUE
      --reference FILE
               Reference sequence FASTA FILE [null]
  -@, --threads INT
               Number of additional threads to use [0]
      --write-index
               Automatically index the output files [off]
      --verbosity INT
               Set level of verbosity

Usage: samtools index -M [-bc] [-m INT] <in1.bam> <in2.bam>...
   or: samtools index [-bc] [-m INT] <in.bam> [out.index]
Options:
  -b, --bai            Generate BAI-format index for BAM files [default]
  -c, --csi            Generate CSI-format index for BAM files
  -m, --min-shift INT   Set minimum interval size for CSI indices to 2^INT [14]
  -M                   Interpret all filename arguments as files to be indexed
  -o, --output FILE    Write index to FILE [alternative to <out.index> in args]
  -@, --threads INT    Sets the number of threads [none]