Repository navigation
Mapping II: Aligners
Remember that depending on the type of data, there will be a better suited mapper. For instance, for DNA sequencing, bowtie2 and bwa are two of the most used short read aligners. While bowtie2 is faster, bwa is more sensitive.
Create a directory called Mappers and copy chr21_dna_R1.fastq.gz, chr21_dna_R2.fastq, chr21_rna_R1.fastq.gz, chr21_rna_R2.fastq from /home/teacher3/Fastq.
We could use the human genome as our reference, but since we are only working with chr21, let's just use this chromosome. This way we can practice generating the needed files for the different tools. Create a directory called db and copy the nucleotide sequence chr21.fasta and its annotation genes.chr21.gtf, from /home/teacher3/db.
Note: Add the time command at the beginning of your cmds, to track the time each process takes. Remember to fill in the table for our evaluation
Let's start by looking at the DNA data.
This is a software package for mapping low-divergent sequences against a large reference genome. BWA-MEM, one of the three algorithms, is generally recommended for its higher speed and accuracy.
The first step
is to create the index for our reference. This has to be done only once for each new reference genome and sometimes it is just a matter of downloading it from iGenomes, for example. Use bwa index.
Then use bwa mem to do the mapping, you will get a SAM file as an output. Convert it to BAM format, sort it and index it. Use samtools to do these manipulations (view, sort and index).
Now generate some statistics using samtools flagstat, samtools idxstats . Save the output files so you can come back to them later.
This is an ultrafast and memory-efficient tool for aligning sequencing reads to long reference sequences. It keeps its memory footprint small and supports gapped, local, and paired-end alignment modes.
Create the index for our reference by using bowtie2-build
Run bowtie2 to map the reads, use samtools again to convert your alignment into BAM, then sort it, index it and generate flagstat and idxstats statistics.
Q. What are your conclusions when comparing these 2 algorithms?
Now let's map the RNA data.
This is a fast splice junction mapper for RNA-Seq reads, that uses Bowtie as a read aligner, then analyzes the mapping results to identify splice junctions between exons.
Here we don't need to create the genome index, since it uses Bowtie2, and we have already created it. So go ahead and run tophat2 using the GTF file with known transcripts. A tophat_out directory will be created, where you can find the alignment file already in BAM format and some statistics.
Let's use picard tools for sorting and indexing (SortSam and BuildBamIndex) the alignment file. But before, there is an extra step when using picard, and that is to add some metadata to the alignment (AddOrReplaceReadGroups). Nowadays more tools are requiring this annotation to assure correct merging, processing and splitting of the data. So it is a good practice to have your data correctly annotated.
Note: picard is not installed globally, otherwise you would use picard.jar
To run picard:
java -jar $picard program_name
Now, use CollectAlignmentSummaryMetrics and BamIndexStats to collect some metrics.
This tool uses a novel strategy for spliced alignments that detects novel splice junctions, thus it makes it more accurate. It is faster than other algorithms but it requires more RAM.
Generate the index files with the following cmd line (change the paths to your files accordingly):
mkdir Star
STAR --runMode genomeGenerate --genomeDir Star --genomeFastaFiles chr21.fasta --sjdbGTFfile genes.chr21.gtf --sjdbOverhang 124
Run the mapping using the GTF file with annotations. Several files will be created, the alignment metrics are stored in Log.final.out. Besides, there are several flags within STAR that can help with the post processing of the alignment file (sort and index) on the fly. But you can also use picard for that, so you choose.
Collect some metrics using picard tools.
Q. What are your conclusions when comparing these 2 algorithms? Are the metrics from the program similar to those you calculated with picard tools?
Q. What are some differences between DNA and RNA mappers
Previous: Mapping I: Intro
Created by Marcela Dávila