Page History

...

...

# 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 -s -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 .

launcher_creator.py -j aln_script.cmds -n aln_script -t 0201:00:00 -w 4 -a OTH21164 -q normal
sbatch --reservation=CoreNGSday4CoreNGS-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

While we're waiting for the job to complete, lets look at the aln_script.cmds file.

...

Exercise #4: Bowtie2

...

alignment - Vibrio cholerae RNA-seq

While we have focused on aligning eukaryotic data, the same tools can be used with prokaryotic data. The major differences are less about the underlying data and much more about the external/public databases that store and distribute reference data. If we want to study a prokaryote, the reference data is usually downloaded from a resource like GenBank.  

...

1. Prepare the vibCho reference index for bowtie2 from GenBank records
2. Align reads using bowtie2, producing a SAM file
3. Convert the SAM file to a BAM file (samtools view) 
4. Sort the BAM file by genomic location (samtools sort)
5. Index the BAM file (samtools index)
6. Gather simple alignment statistics (samtools flagstat and samtools idxstatidxstats)

Obtaining the GenBank records

First prepare a directory to work infor the vibCho fasta, and change to it:

mkdir -p $SCRATCH/core_ngs/references/vibChofasta
cd $SCRATCH/core_ngs/references/vibChofasta

V. cholerae has two chromosomes. We download each separately.

...

mkdir -p $SCRATCH/core_ngs/references/vibChofasta
cd $SCRATCH/core_ngs/references/vibChofasta
cp $CORENGS/ncbireferences/vibChofasta/NC_* .

Once you have the 4 files locally in your $SCRATCH/core_ngs/references/vibCho directory, combine them using cat:

cd $SCRATCH/core_ngs/references/vibChofasta
cat NC_01258[23].fa   > vibCho.O395.fa
cat NC_01258[23].gff3 > vibCho.O395.gff3

# verify there are 2 contigs in vibCho.O395.fa
grep -P '^>' vibCho.O395.fa

Now we have a reference sequence file that we can use with the bowtie2 reference builder, and ultimately align sequence data against.

...

idev -m 120 -A OTH21164 -N 1 -r CoreNGSday4
 CoreNGS-Thu
# or
idev -m 90 -A OTH21164 -N 1 -p development 

Go ahead and load the bowtie2 module so we can examine some help pages and options.

module load biocontainers
module load bowtie2

...

First create a directory specifically for the bowtie2 index, then build the index using bowtie-build.

mkdir -p $SCRATCH/core_ngs/references/fasta
cd $SCRATCH/core_ngs/references/fasta
cp $CORENGS/references/fasta/vibCho* .

mkdir -p $SCRATCH/core_ngs/references/bt2/vibCho
cd $SCRATCH/core_ngs/references/bt2/vibCho

# Symlink to the fasta file you created using relative path syntax
ln -sf $SCRATCH/core_ngs/references../../fasta/vibCho.O395.fa
# or, to catch up:
ln -sf $CORENGS/references/vibCho.O395.fa

bowtie2-build vibCho
bowtie2-build vibCho.O395.fa vibCho.O395

This should also go pretty fast. You can see the resulting files using ls like before.

Performing the bowtie2 alignment

We'll set up a new directory to perform the V. cholerae data alignment. But first make sure you have the FASTQ file to align and the vibCho bowtie2 indexMake sure you're in an idev session with the bowtie2 BioContainers module loaded:

# Get the FASTQ to align
mkdir -p $SCRATCH/core_ngs/alignment/fastq
cp $CORENGS/alignment/*fastq.gz $SCRATCH/core_ngs/alignment/fastq/

# Set up the bowtie2 index
mkdir -p $SCRATCH/core_ngs/references/bt2/vibCho
cp $CORENGS/idx/bt2/vibCho/*.*  $SCRATCH/core_ngs/references/bt2/vibCho/

idev -m 120 -A OTH21164  -N 1 -r CoreNGS-Thu
# or
idev -m 90 -A OTH21164 -N 1 -p development

module load biocontainers
module load bowtie2

We'll set up a new directory to perform the V. cholerae data alignment. But first make sure you have the FASTQ file to align and the vibCho bowtie2 index.Make sure you're in an idev session with the bowtie2 BioContainers module loaded:

idev# Get a pre-mbuilt 120vibCho -Aindex OTH21164if you -N 1 -r CoreNGSday4
module load biocontainers
module load bowtie2didn't already build one
mkdir -p $SCRATCH/core_ngs/references/bt2/vibCho
cp $CORENGS/references/bt2/vibCho/*.*  $SCRATCH/core_ngs/references/bt2/vibCho/

# Get the FASTQ to align
mkdir -p $SCRATCH/core_ngs/alignment/fastq
cp $CORENGS/alignment/*fastq.gz  $SCRATCH/core_ngs/alignment/fastq/

Now Now set up a directory to do this alignment, with symbolic links to the bowtie2 index directory and the directory containing the FASTQ to align:

...

bowtie2 [options]* -x <bt2-idx> {-1 <m1> -2 <m2> | -U <r>} [-S -U <r>} [-S <sam>]

...

<sam>]

So execute this bowtie2 global, single-end alignment command:

mkdir -p $SCRATCH/core_ngs/references/fasta
cp $CORENGS/references/fasta/vibCho*  $SCRATCH/core_ngs/references/fasta/

mkdir -p $SCRATCH/core_ngs/references/bt2/vibCho
cp $CORENGS/references/bt2/vibCho/*.*  $SCRATCH/core_ngs/references/bt2/vibCho/

mkdir -p $SCRATCH/core_ngs/alignment/vibCho
cd $SCRATCH/core_ngs/alignment/vibCho
ln -sf ../../references/bt2/vibCho
ln -sf ../../alignment/fastq fqmkdir -p $SCRATCH/core_ngs/alignment/vibCho

cd $SCRATCH/core_ngs/alignment/vibCho
bowtie2 -x vibCho/vibCho.O395 -U fq/cholera_rnaseq.fastq.gz \
  -S cholera_rnaseq.sam 2>&1 | tee aln_global.log

Notes:

• -x  vibCho vibCho/vibCho.O395.fa – prefix path of index files
• -U fq/cholera_rnaseq.fastq.gz – FASTQ file for single-end (Unpaired) alignment
• -S cholera_rnaseq.sam – tells bowtie2 to report alignments in SAM format to the specified file
• 2>&1 redirects standard error to standard output
  • while the alignment data is being written to the cholera_rnaseq.sam file, bowtie2 will report its progress to standard error.
• | tee aln.log takes the bowtie2 progress output and pipes it to the tee program
  • tee takes its standard input and writes it to the specified file and also to standard output
  • that way, you can see the progress output now, but also save it to review later (or supply to MultiQC)

Since the FASTQ file is not large, this should not take too long, and you will see progress output like this:

...

After bowtie2 came out with a local alignment option, it wasn't long before bwa developed its own local alignment algorithm called BWA-MEM (for Maximal Exact Matches), implemented by the bwa mem command. 

bwa mem has the following advantages:

• It provides the simplicity of using bwa without the complexities of local alignment
• It can align different portions of a read to different locations on the genome
  • In a total RNA-seq experiment, reads will (at some frequency) span a splice junction themselves
    • or a pair of reads in a paired-end library will fall on either side of a splice junction.
  • We want to be able to align these splice-adjacent reads for many reasons, from accurate transcript quantification to novel fusion transcript discovery.

This exercise will align a human total RNA-seq dataset composed (by design) almost exclusively of that includes numerous reads that cross splice junctions.

...

In the transcriptome-aware alignment above, reads that span splice junctions are reported in the SAM file with genomic coordinates that start in the first exon and end in the second exon (the CIGAR string uses the N operator, e.g. 30M1000N60M).

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 of the two will be marked as secondary (0x100 flag).

...

First set up our working directory for this alignment. Since it takes a long time to build a bwa index for a large genome (here human hg38/GRCh38), we'll use one that the BioITeam maintains in its /work/projects/BioITeam/ref_genome area.

# Make sure you're in an idev session
idev -m 120 -N 1 -A OTH21164 -r CoreNGSday4  CoreNGS-Thu
# or
idev -m 90 -N 1 -A OTH21164 -p development

# 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

...

• The bwa mem alignment
  • the program's progress output (on standard error) is redirected to a log file (2>hs_rna.bwamem.log)
  • its alignment records (on standard output) is piped to the next step (conversion to BAM)
• Conversion of bwa mem's SAM output to BAM format
  • recall that the -b option to samtools view says to output in BAM format
• Sorting the BAM file
  • samtools sort takes the binary output from samtools view and writes a sorted BAM file.

Because the progress output is being redirected to a log file, we won't see anything until the command completes. Then you should have a human_rnaseq.sort.bam file and an hs_rna.bwamem.log logfile.

...

cd $SCRATCH/core_ngs/alignment/bwamem
zcat ./fastq/human_rnaseq.fastq.gz | wc -l | awk '{print $1/4}'

cd $SCRATCH/core_ngs/alignment/bwamem
samtools flagstat human_rnaseq.sort.bam | tee hs_rnaseq.flagstat.txt

133570 + 0 in total (QC-passed reads + QC-failed reads)
33570 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
133450 + 0 mapped (99.91% : N/A)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (N/A : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (N/A : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

...

...