| Tip | ||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| ||||||||||||||
Use our summer school reservation (core-ngs-class-0605) for today when submitting batch jobs to get higher priority on the ls6 normal queue.
Note that today's reservation name (core-ngs-class-0605) is different from the TACC allocation/project for this class, which is OTH21164. |
...
| Reference | Species | Base Length | Contig Number | Source | Download |
|---|---|---|---|---|---|
| hg19 | Human | 3.1 Gbp | 25 (really 93) | UCSC | UCSC GoldenPath |
| sacCer3 | Yeast | 12.2 Mbp | 17 | UCSC | UCSC GoldenPath |
| mirbase v20 | Human subset | 160 Kbp | 1908 | miRBase | miRBase Downloads |
| vibCho (O395ASM836960v1) | Vibrio cholerae | ~4.2 Mbp2 | 3 | GenBank | GenBank Downloads |
Searching genomes is computationally hard work and takes a long time if done on linear genomic sequence. So aligners require that references first be indexed to accelerate lookup. The aligners we are using each require a different index, but use the same method (the Burrows-Wheeler Transform) to get the job done.
...
| Expand | ||
|---|---|---|
| ||
Examine the sub-command's usage:
If you just type The pipe ( | ) only connects standard output from bwa aln to standard input of the more command – but no standard output is generated. And the standard error text that is generated just goes to your Terminal, bypassing more. |
...
| Tip | ||
|---|---|---|
| ||
Double check that output was written by doing When redirection ( > ) to a file is used, the target file is created whether or not the command produced any output! |
...
Next we use the bwa sampe command to pair the reads and output SAM format data. Just type bwa sampe 2>&1 | more to see its usage.
For this command you provide the same reference index prefix as for bwa aln, along with the two .sai files and the two original FASTQ files. Also, bwa writes its output to standard output by default, so redirect that to a .sam file.
...
Exercise: How many lines does the SAM file have? How does this compare to the number of input sequences (R1+R2)?
| Expand | ||
|---|---|---|
| ||
The alignment SAM file will contain records for both R1 and R2 reads, so we need to count sequences in both files.
So the SAM file has 18 more lines than the R1+R2 total. These are the header records that appear before any alignment records. |
...
The samtools view utility provides a way of converting between SAM (text) and BAM (binary, compressed) format. It also provides many, many other functions which we will discuss lster. To get a preview, execute samtools view 2>&1 | more. You should see:
| Code Block | ||
|---|---|---|
| ||
Usage: samtools view [options] <in.bam>|<in.sam>|<in.cram> [region ...]
Options:
-b output BAM
-C output CRAM (requires -T)
-1 use fast BAM compression (implies -b)
-u uncompressed BAM output (implies -b)
-h include header in SAM output
-H print SAM header only (no alignments)
-c print only the count of matching records
-o FILE output file name [stdout]
-U FILE output reads not selected by filters to FILE [null]
-t FILE FILE listing reference names and lengths (see long help) [null]
-X include customized index file
-L FILE only include reads overlapping this BED FILE [null]
-r STR only include reads in read group STR [null]
-R FILE only include reads with read group listed in FILE [null]
-d STR:STR
only include reads with tag STR and associated value STR [null]
-D STR:FILE
only include reads with tag STR and associated values listed in
FILE [null]
-q INT only include reads with mapping quality >= INT [0]
-l STR only include reads in library STR [null]
-m INT only include reads with number of CIGAR operations consuming
query sequence >= INT [0]
-f INT only include reads with all of the FLAGs in INT present [0]
-F INT only include reads with none of the FLAGS in INT present [0]
-G INT only EXCLUDE reads with all of the FLAGs in INT present [0]
-s FLOAT subsample reads (given INT.FRAC option value, 0.FRAC is the
fraction of templates/read pairs to keep; INT part sets seed)
-M use the multi-region iterator (increases the speed, removes
duplicates and outputs the reads as they are ordered in the file)
-x STR read tag to strip (repeatable) [null]
-B collapse the backward CIGAR operation
-? print long help, including note about region specification
-S ignored (input format is auto-detected)
--no-PG do not add a PG line
--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
of OPTION or OPTION=VALUE
-T, --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 |
...
| Expand | ||
|---|---|---|
| ||
Note that samtools (like bwa) writes its help to standard error, but less and more only accept input on standard input. So the syntax redirecting standard error to standard input must be used before the pipe to less or more.
then search for "header" ( /header ) |
...
| Code Block | ||||
|---|---|---|---|---|
| ||||
samtools view yeast_pe.bam | cut -f 1-4 | more |
(e.g. samtools view yeast_pe.bam | cut -f 1-4 | more), you You will notice that read names appear in adjacent pairs (for the R1 and R2), in the same order they appeared in the original FASTQ file. This means the records are in name order and, searching through the file is very inefficient. samtools sort re-orders entries in the SAM file either by locus (contig name + coordinate position) or by read name.
If you execute samtools sort 2>&1 | more, you see its help page:
| Code Block | ||
|---|---|---|
| ||
Usage: samtools sort [options...] [in.bam]
Options:
-l INT Set compression level, from 0 (uncompressed) to 9 (best)
-m INT Set maximum memory per thread; suffix K/M/G recognized [768M]
-n Sort by read name
-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
--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
of OPTION or OPTION=VALUE
--reference FILE
Reference sequence FASTA FILE [null]
-@, --threads INT
Number of additional threads to use [0] |
...