Table of Contents |
---|
Overview
Most approaches for predicting structural variants require you to have paired-end or mate-pair reads. They use the distribution of distances separating these reads to find outliers and also look at pairs with incorrect orientations. As mentioned during several of the presentations, many researchers choose to ignore these types of mutations and combined with the increased difficulty of accurately identifying them, the community is less settled on the "best" way to analyze them. Here we present a tutorial on SVDetect based on the quality of its instructions, and easy of installation despite its use of relatively hefty configuration files.
...
Example: E. coli genome with structural variation
Here's an E. coli genome re-sequencing sample where a key mutation producing a new structural variant was responsible for a new phenotype involving citrate, one of Dacia's favorite topicssomething the Barrick lab has studied.
Code Block | ||||
---|---|---|---|---|
| ||||
cds cp -r $BI/gva_course/structural_variation/data BDIBGVA_sv_tutorial cd BDIBGVA_sv_tutorial |
This is Illumina mate-paired data (having a larger insert size than paired-end data) from genome re-sequencing of an E. coli clone.
File Name | Description | Sample |
---|---|---|
| Paired-end Illumina, First of mate-pair, FASTQ format | Re-sequenced E. coli genome |
| Paired-end Illumina, Second of mate-pair, FASTQ format | Re-sequenced E. coli genome |
| Reference Genome in FASTA format | E. coli B strain REL606 |
NC_012967.1.lengths | Simple tab delimtered file based on the size of the reference needed for SVDetect so you don't have to create it yourself |
Map data using bowtie2
First we need to (surprise!) map the data. This will hopefully reinforce the bowtie2 tutorial you just completed, but if you are feeling adventurous you could use BWA as optional reinforcement.
Warning | |||||||
---|---|---|---|---|---|---|---|
| |||||||
Make sure you are on an idev node using the command: showq -u
|
Code Block |
---|
bowtie2-build NC_012967.1.fasta NC_012967.1 bowtie2 -p 48 -X 5000 --rf -x NC_012967.1 -1 61FTVAAXX_2_1.fastq -2 61FTVAAXX_2_2.fastq -S 61FTVAAXX.sam |
...
--rf
tells bowtie2 that your read pairs are in the "reverse-forward" orientation of a mate-pair library-X 5000
tells bowtie2 to not mark read pairs as discordant unless their insert size is greater than 5000 bases.
...
You may notice that these commands complete pretty quickly. Always remember speed is not necessarily representative of how taxing something is for TACC's head node, and always try to be a good TACC citizen and do as much as you can on idev nodes or as job submissions
Expand | ||
---|---|---|
| ||
bowtie doesn't output discordantly mapped pairs! |
Run SVDetect
The first step is to look at all mapped read pairs and whittle down the list only to those that have an unusual insert sizes (distances between the two reads in a pair).
...
Expand | ||
---|---|---|
| ||
|
SVDetect demonstrates a common strategy in some programs with complex input where instead of including a lot of options on the command line, it reads in a simple text file that sets all of the required options. Lets look at how to create a configuration file:
Note |
---|
You'll need to substitute the output of the pwd command on lines 7 and 8 below in front of "/full/path/to/" , but you will need to leave "61FTVAAXX.ab.sam" and "NC_012967.1.lengths". This is often a source of problems in this tutorial |
Code Block | ||||
---|---|---|---|---|
| ||||
<general> input_format=sam sv_type=all mates_orientation=RF read1_length=35 read2_length=35 mates_file=/full/path/to/61FTVAAXX.ab.sam cmap_file=/full/path/to/NC_012967.1.lengths num_threads=48 </general> <detection> split_mate_file=0 window_size=2000 step_length=1000 </detection> <filtering> split_link_file=0 nb_pairs_threshold=3 strand_filtering=1 </filtering> <bed> <colorcode> 255,0,0=1,4 0,255,0=5,10 0,0,255=11,100000 </colorcode> </bed> |
...
You also need to make sure you have a copy of the tab-delimited file of chromosome lengths named NC_012967.1.lengths. YOU CAN NOT COPY PASTE THIS COMMAND into a new file for 2 reasons! The first reason you can't copy the command is the tab characters don't translate correctly. The second is the nano text editor is not available on the compute nodes, and learning to use the vim editor for a single line is way more work than it is worth. Make sure the the NC_012967.1.lengths file you copied has the following structure up to the comment, and that the <tab> is replaced with an actual tab character.
...
Take a look at the resulting file: 61FTVAAXX.ab.sam.links.filtered.sv.txt
. Another downside of command line applications is that while you can print files to the screen, the formatting is not always the nicest. On the plus side in 95% of cases, you can directly copy the output from the terminal window to excel and make better sense of what the columns actually are
WeI've highlighted a few lines below:
...
Expand | |||||||||
---|---|---|---|---|---|---|---|---|---|
| |||||||||
Optional: Install SVDetectWe have installed SVdetect for you already as installation is a bit difficult (though still much easier than the alternatives listed in the introduction). You can verify it's location using which SVDetect in your Install SVDetect scriptsNavigate to the SVDetect project page More information: Download the code onto TACC.
Move the Perl scripts and make them executable
Install required Perl modulesSVdetect requires a few Perl modules to be installed. In the default TACC environment, you can use the cpan shell to install most well-behaved Perl modules (with the exception of some complicated ones that require other libraries to be installed or things to compile). Here's how:
|