...
| Code Block | ||
|---|---|---|
| ||
align_bwa_illumina.sh 2022_0506_0510 Align Illumina SE or PE data with bwa. Produces a sorted, indexed, duplicate-marked BAM file and various statistics files. Usage: align_bwa_illumina.sh <aln_mode> <in_file> <out_pfx> <assembly> [ paired trim_sz trim_sz2 seq_fmt qual_fmt ] Required arguments: aln_mode Alignment mode, either global (bwa aln) or local (bwa mem). in_file For single-end alignments, path to input sequence file. For paired-end alignments using fastq, path to the the R1 fastq file which must contain the string 'R1' in its name. The corresponding 'R2' must have the same path except for 'R1'. out_pfx Desired prefix of output files in the current directory. assembly One of hg38, hg19, hg38, mm10, mm9, sacCer3, sacCer1, ce11, ce10, danRer7, hs_mirbase, mm_mirbase, or reference index prefix. Optional arguments: paired 0 = single end alignment (default); 1 = paired end. trim_sz Size to trim reads to. Default 0 (no trimming) trim_sz2 Size to trim R2 reads to for paired end alignments. Defaults to trim_sz seq_fmt Format of sequence file (fastq, bam or scarf). Default is fastq if the input file has a '.fastq' extension; scarf if it has a '.sequence.txt' extension. qual_type Type of read quality scores (sanger, illumina or solexa). Default is sanger for fastq, illumina for scarf. Environment variables: show_only 1 = only show what would be done (default not set) aln_args other bowtie2 options (e.g. '-T 20' for mem, '-l 20' for aln) no_markdup 1 = don't mark duplicates (default 0, mark duplicates) run_fastqc 1 = run fastqc (default 0, don't run). Note that output will be in the directory containing the fastq files. keep 1 = keep unsorted BAM (default 0, don't keep) bwa_bin BWA binary to use. Default bwa 0.7.x. Note that bwa 0.6.2 or earlier should be used for scarf and other short reads. also: NUM_THREADS, BAM_SORT_MEM, SORT_THREADS, JAVA_MEM_ARG Examples: align_bwa_illumina.sh local ABC_L001_R1.fastq.gz my_abc hg38 1 align_bwa_illumina.sh global ABC_L001_R1.fastq.gz my_abc hg38 1 50 align_bwa_illumina.sh global sequence.txt old sacCer3 0 '' '' scarf solexa |
...
| Code Block | ||||
|---|---|---|---|---|
| ||||
# Make sure you're not in an idev session by looking at the hostname hostname # If the hostname looks like "c455-004.ls6.tacc.utexas.edu", exit the idev session # Copy over the Yeast data if needed mkdir -p $SCRATCH/core_ngs/alignment/fastq cp $CORENGS/alignment/Sample_Yeast*.gz $SCRATCH/core_ngs/alignment/fastq/ # Make a new alignment directory for running these scripts mkdir -p $SCRATCH/core_ngs/alignment/bwa_script cd $SCRATCH/core_ngs/alignment/bwa_script ln -ssf -f ../fastq # Copy the alignment commands file and submit the batch job cd $SCRATCH/core_ngs/alignment/bwa_script cp $CORENGS/tacc/aln_script.cmds . # Use -a TRA23004 below if OTH21164 hasn't been working for you... launcher_creator.py -j aln_script.cmds -n aln_script -t 01:00:00 -w 4 -a OTH21164 -q normal sbatch --reservation=CoreNGS-Thu aln_script.slurm # or launcher_creator.py -j aln_script.cmds -n aln_script -t 01:00:00 -w 4 -a OTH21164 -q development sbatch aln_script.slurm showq -u |
...
Since this was a paired end alignment there is paired-end specific information reported. Note that this Yeast dataset was of poor quality, especially the R2 reads. You can see this by the relatively low R1 alignment rate (~63%) and the really low R2 alignment rate (~28%).
You can also view statistics on insert sizes for properly paired reads in the bwa_global.iszinfo.txt file. This tells you the average (mean) insert size, standard deviation, mode (most common value), and fivenum values (minimum, 1st quartile, median, 3rd quartile, maximum).
...
| Code Block | ||||
|---|---|---|---|---|
| ||||
idev -m 120 -A OTH21164 -N 1 -r CoreNGS-Thu # or -A TRA23004 # or idev -m 90 -A OTH21164 -N 1 -p development # or -A TRA23004 |
Go ahead and load the bowtie2 module so we can examine some help pages and options.
...
| Code Block | ||
|---|---|---|
| ||
idev -m 120 -A OTH21164 -N 1 -r CoreNGS-Thu # or -A TRA23004 # or idev -m 90 -A OTH21164 -N 1 -p development # or -A TRA23004 module load biocontainers module load bowtie2 |
...
| Expand | |||||||
|---|---|---|---|---|---|---|---|
| |||||||
Reports these alignment statistics:
Interestingly, the local alignment rate here is lower than we saw with the global alignment. Usually local alignments have higher alignment rates than corresponding global ones. |
...
Using bwa mem for RNA-seq alignment is sort of a "poor man's" RNA-seq alignment method. Real splice-aware aligners like tophat2STAR, hisat2 or STAR tophat have more complex algorithms (as shown below) – and take a lot more time!
...
BWA MEM does not know about the exon structure of the genome. But it can align different sub-sections of a read to two different locations, producing two alignment records from one input read (one . One of the two will be marked as secondary (0x100 flag).
...
| Code Block | ||||
|---|---|---|---|---|
| ||||
# Make sure you're in an idev session idev -m 120 -N 1 -A OTH21164 -r CoreNGS-Thu # or -A TRA23004 # or idev -m 90 -N 1 -A OTH21164 -p development # or -A TRA23004 # Load the modules we'll need module load biocontainers module load bwa module load samtools # Copy over the FASTQ data if needed mkdir -p $SCRATCH/core_ngs/alignment/fastq cp $CORENGS/alignment/*.gz $SCRATCH/core_ngs/alignment/fastq/ # Make a new alignment directory for running these scripts cds mkdir -p core_ngs/alignment/bwamem cd core_ngs/alignment/bwamem ln -sf ../fastq ln -sf /work/projects/BioITeam/ref_genome/bwa/bwtsw/hg38 |
Now take a look at bwa mem usage (type bwa mem with no arguments, or bwa mem 2>&1 | more). The most important parameters are the following:
...
| Code Block | ||
|---|---|---|
| ||
cd $SCRATCH/core_ngs/alignment/bwamem
bwa mem -M hg38/hg38.fa fastq/human_rnaseq.fastq.gz 2>hs_rna.bwamem.log | \
samtools view -b | \
samtools sort -O BAM -T human_rnaseq.tmp > human_rnaseq.sort.bam |
...