Content Comparison

...

Code Block

title	This is one command. Put it all on one line.

samtools mpileup -u -f NC_012967.1.fasta.fai SRR030257.sorted.bam > SRR030257.bcf

What are all the options doing? Try calling samtools mpileup without any options to see if you can figure it out before clicking below to

Expand

Option	purpose
-u	generates uncompressed BCF output
-f NC_012967.1.fasta.fai	faidx indexed reference sequence file
SRR030257.sorted.bam	BAM input file to calcluate pileups from
> SRR030257.bcf	Direct output to SRR030257.bcf

The samtools mpileup command will take a few minutes to run. You might consider As practice for a fairly common occurrence when working with the iDEV environment, once the command is running, you should try putting it in the background by pressing control-z and then typing the command bg so that you can do some other things in this terminal window at the same time – like starting to run variant calling on the BWA or bowtie mapping results. Remember, there are still many other processors available on this node!

Convert BCF to Once the mpileup command is complete, convert the BCF file to a "human-readable VCF:" VCF file using the bcftools command (the which command will tell you where this command is located and examination of that path should tell you how you have access to it).

Code Block

title	This is one command. Put it all on one line.

bcftools view -v -c -g SRR030257.bcf > SRR030257.vcf

What are these options doing?

Expand

title	See if you can start with the base command "bcftools" and figure out what each option is doing. Click here when you think you know.

Option	purpose
view	specific command to be executed by bcftools
-v	output potential variant sites only
-c	SNP calling
-g	call genotypes at variant sites
SRR030257.bcf	input bcf file
> SRR030257.vcf	output as a vcf file

Take a look at the samtools_bowtie2/SRR030257.vcf file using less. It has a nice header explaining what the columns mean. Below this are the rows of data describing potential genetic variants.

...

Follow the same directions to call variants in the BWA or Bowtie mapped reads.Just Just be sure you don't write over your old files. Maybe create new directories like samtools_bwa and samtools_bowtie for the output in each case.You You could also try running all of the commands from inside of the samtools_bwa directory, just for a change of pace.

...

VCF format has alternative Allele Frequency tags denoted by AF1. AF= Try the following command to see what values we have in our files.

Code Block
grep AF1 samtools_bowtie2/AF= SRR030257.vcf

Expand

title	Optional:

...

For the data we are dealing with, predictions with an allele frequency not equal to 1 are not really applicable. (The reference genome is haploid. There aren't any heterozygotes.) How can we remove these lines from the file?

...

What does the -v flag do in grep?

codeSolution...

Expand

Hint...

Try looking at grep --help to see what you can come up with.

Code Block

language	bash
title	Here for answer
collapse	true

grep -v *something*  # The -v flag inverts the match effecitvely showing you everything that does not match your input

Expand

Solution...

title	Going farther

Code Block
cat input.vcf \| grep

AF1

AF=1 > output.vcf

Is not practical, since we will lose vital VCF formatting and may not be able to use this file in the future.

Code Block
cat input.vcf \| grep -v

AF1

AF=0 > output.vcf

Will preserve all lines that don't have a

AF1

AF=0 value and is one way of doing this.

Code Block
sed -i '/AF1=0/ d' input.vcf

Is a way of doing it in-line and not requiring you to make another file. (But it writes over your existing file!)

Comparing the results of different mappers using bedtools

Often you want to compare the results of variant calling on different samples or using different pipelines. Bedtools is a suite of utility programs that work on a variety of file formats, one of which is conveniently VCF format. It provides many ways of slicing, dicing, and comparing the information in VCF files. For example, we can use it to find out what predictions are the same and which are different from the variant calling on reads mapped with different programs if you generated VCF files for each one.Set Set up a new output directory and copy the respective VCF files to it, renaming them so that we know where they came from:

Code Block

language	bash
title	If you have done any of the optional other mapping tutorials, consider the following comparisons. Remember the use of cp -n is useful to make sure you don't overwrite any existing files.

mkdir comparison
cp -n samtools_bowtie2/SRR030257.vcf comparison/bowtie2.vcf
cp -n samtools_bwa/SRR030257.vcf comparison/bwa.vcf
cp -n samtools_bowtie/SRR030257.vcf comparison/bowtie.vcf
cd comparison

Use the subcommands bedtools intersect and bedtools subtract we can find equal and different predictions between mappers. Try to figure out how to to do this on your own first. Hint: Remember that adding > output.vcf to the end of a command will pipe the output that is to the terminal into a file, so that you can save it.

Load Bedtools.

code

Code Block

I'm stuck? Show me the commands...

language

Expand

bash
title	load bedtools
collapse	true

module load bedtools

code

Code Block

collapse
language	bash
title	Finding common mutations.

true

bedtools intersect -a bowtie2.vcf -b bwa.vcf > common_bowtie2_bwa.vcf

code

Code Block

collapse
language	bash
title	Finding mutations that are unique for each mapper.

true

bedtools subtract -a bowtie2.vcf -b common_bowtie2_bwa.vcf > unique_bowtie2.vcf
bedtools subtract -a bwa.vcf -b common_bowtie2_bwa.vcf > unique_bwa.vcf

Further Optional Exercises

...

Version	Old Version 3	New Version 4
Changes made by	Deatherage, Daniel E	Deatherage, Daniel E
Saved on	May 25, 2015	May 25, 2015

Versions Compared

Key

Comparing the results of different mappers using bedtools

Further Optional Exercises