小生这厢有礼了(BioFaceBook Personal Blog) » 二代测序

Metagenomics standard operating procedure v2

szypanther — Sun, 23 Aug 2020 03:28:53 +0000

http://www.360doc.com/content/20/0823/10/71250389_931761136.shtml

https://github.com/LangilleLab/microbiome_helper/wiki/Metagenomics-standard-operating-procedure-v2

Note that this workflow is continually being updated. If you want to use the below commands be sure to keep track of them locally.

Last updated: 9 Oct 2018 (see revisions above for earlier minor versions and Metagenomics SOPs on the SOP for earlier major versions)

Important points to keep in mind:

The below options are not necessarily best for your data; it is important to explore different options, especially at the read quality filtering stage.
If you are confused about a particular step be sure to check out the pages under Metagenomic Resources on the right side-bar.
You should check that only the FASTQs of interest are being specified at each step (e.g. with the ls command). If you are re-running commands multiple times or using optional steps the input filenames and/or folders could differ.
The tool parallel comes up several times in this workflow. Be sure to check out our tutorial on this tool here.
You can always run parallel with the --dry-run option to see what commands are being run to double-check they are doing what you want!

This workflow starts with raw paired-end MiSeq or NextSeq data in demultiplexed FASTQ format assumed to be located within a folder called raw_data.

1. First Steps

1.1 Join lanes together (all Illumina, except MiSeq)

Concatenate multiple lanes of sequencing together (e.g. for NextSeq data). If you do NOT do this step, remember to change cat_lanes to raw_data in the further commands below that have assumed you are working with the most common type of metagenome data.

concat_lanes.pl raw_data/* -o cat_lanes -p 4

1.2 Inspect read quality

Run fastqc to allow manual inspection of the quality of sequences.

mkdir fastqc_out
fastqc -t 4 cat_lanes/* -o fastqc_out/

1.3 (Optional) Stitch F+R reads

Stitch paired-end reads together (summary of stitching results are written to “pear_summary_log.txt”). Note: it is important to check the % of reads assembled. It may be better to concatenate the forward and reverse reads together if the assembly % is too low (see step 2.2).

run_pear.pl -p 4 -o stitched_reads cat_lanes/*

If you don’t stitch your reads together at this step you will need to unzip your FASTQ files before continuing with the below commands.

2. Read Quality-Control and Contaminant Screens

2.1 Running KneadData

Use kneaddata to run pre-processing tools. First Trimmomatic is run to remove low quality sequences. Then Bowtie2 is run to screen out contaminant sequences. Below we are screening out reads that map to the human or PhiX genomes. Note KneadData is being run below on all unstitched FASTQ pairs with parallel, you can see our quick tutorial on this tool here. For a detailed breakdown of the options in the below command see this page. The forward and reverse reads will be specified by “_1″ and “_2″ in the output files, ignore the “R1″ in each filename. Note that the \ characters at the end of each line are just to split the command over multiple lines to make it easier to read.

parallel -j 1 --link 'kneaddata -i {1} -i {2} -o kneaddata_out/ \
-db /home/shared/bowtiedb/GRCh38_PhiX --trimmomatic /usr/local/prg/Trimmomatic-0.36/ \
-t 4 --trimmomatic-options "SLIDINGWINDOW:4:20 MINLEN:50" \
--bowtie2-options "--very-sensitive --dovetail" --remove-intermediate-output' \
 ::: cat_lanes/*_R1.fastq ::: cat_lanes/*_R2.fastq

Clean up the output directory (helps downstream commands) by moving the discarded sequences to a subfolder:

mkdir kneaddata_out/contam_seq

mv kneaddata_out/*_contam*.fastq kneaddata_out/contam_seq

You can produce a logfile summarizing the kneadData output with this command:

kneaddata_read_count_table --input kneaddata_out --output kneaddata_read_counts.txt

2.2 Concatenate unstitched output (Omit if data stitched above)

If you did not stitch your paired-end reads together with pear, then you can concatenate the forward and reverse FASTQs together now. Note this is done after quality filtering so that both reads in a pair are either discarded or retained. It’s important to only specify the “paired” FASTQs outputted by kneaddata in the below command.

concat_paired_end.pl -p 4 --no_R_match -o cat_reads kneaddata_out/*_paired_*.fastq

You should check over the commands that are printed to screen to make sure the correct FASTQs are being concatenated together.

3. Determine Functions with HUMAnN2

3.1 Run HUMAnN2

Run humann2 with parallel to calculate abundance of UniRef90 gene families and MetaCyc pathways. If you are processing environment data (e.g. soil samples) the vast majority of the reads may not map using this approach. Instead, you can try mapping against the UniRef50 database (which you can point to with the --protein-database option).

mkdir humann2_out

parallel -j 4 'humann2 --threads 1 --input {} --output humann2_out/{/.}' ::: cat_reads/*fastq

3.2 Merge individual sample data together

Join HUMAnN2 output per sample into one table.

mkdir humann2_final_out

humann2_join_tables -s --input humann2_out/ --file_name pathabundance --output humann2_final_out/humann2_pathabundance.tsv
humann2_join_tables -s --input humann2_out/ --file_name pathcoverage --output humann2_final_out/humann2_pathcoverage.tsv
humann2_join_tables -s --input humann2_out/ --file_name genefamilies --output humann2_final_out/humann2_genefamilies.tsv

3.3 Table output normalization

Re-normalize gene family and pathway abundances (so that all samples are in units of copies per million).

humann2_renorm_table --input humann2_final_out/humann2_pathabundance.tsv --units cpm --output humann2_final_out/humann2_pathabundance_cpm.tsv
humann2_renorm_table --input humann2_final_out/humann2_genefamilies.tsv --units cpm --output humann2_final_out/humann2_genefamilies_cpm.tsv

3.4 Separate out taxonomic contributions

Split HUMAnN2 output abundance tables in stratified and unstratified tables (stratified tables include the taxa associated with a functional profile).

humann2_split_stratified_table --input humann2_final_out/humann2_pathabundance_cpm.tsv --output humann2_final_out
humann2_split_stratified_table --input humann2_final_out/humann2_genefamilies_cpm.tsv --output humann2_final_out
humann2_split_stratified_table --input humann2_final_out/humann2_pathcoverage.tsv --output humann2_final_out

3.5 Format STAMP function file

Convert unstratified HUMAnN2 abundance tables to STAMP format by changing header-line. These commands remove the comment character and the spaces in the name of the first column. Trailing descriptions of the abundance datatype are also removed from each sample’s column name.

sed 's/_Abundance-RPKs//g' humann2_final_out/humann2_genefamilies_cpm_unstratified.tsv | sed 's/# Gene Family/GeneFamily/' > humann2_final_out/humann2_genefamilies_cpm_unstratified.spf
sed 's/_Abundance//g' humann2_final_out/humann2_pathabundance_cpm_unstratified.tsv | sed 's/# Pathway/Pathway/' > humann2_final_out/humann2_pathabundance_cpm_unstratified.spf
sed 's/_Coverage//g' humann2_final_out/humann2_pathcoverage_unstratified.tsv | sed 's/# Pathway/Pathway/' > humann2_final_out/humann2_pathcoverage_unstratified.spf

3.6 Extract MetaPhlAn2 taxonomic compositions

Since HUMAnN2 also runs MetaPhlAn2 as an initial step, we can use the output tables already created to get the taxonomic composition of our samples. First we need to gather all the output MetaPhlAn2 results per sample into a single directory and then merge them into a single table using MetaPhlAn2’s merge_metaphlan_tables.py command. After this file is created we can fix the header so that each column corresponds to a sample name without the trailing “_metaphlan_bugs_list” description. Note that MetaPhlAn2 works best for well-characterized environments, like the human gut, and has low sensitivity in other environments.

mkdir metaphlan2_out
cp humann2_out/*/*/*metaphlan_bugs_list.tsv metaphlan2_out/
/usr/local/metaphlan2/utils/merge_metaphlan_tables.py metaphlan2_out/*metaphlan_bugs_list.tsv > metaphlan2_merged.txt
sed -i 's/_metaphlan_bugs_list//g' metaphlan2_merged.txt

3.7 Format STAMP taxonomy file

Lastly we can convert this MetaPhlAn2 abundance table to STAMP format

metaphlan_to_stamp.pl metaphlan2_merged.txt > metaphlan2_merged.spf

Vipie: web pipeline for parallel characterization of viral populations from multiple NGS samples

szypanther — Mon, 16 Sep 2019 06:12:35 +0000

Background

Next generation sequencing (NGS) technology allows laboratories to investigate virome composition in clinical and environmental samples in a culture-independent way. There is a need for bioinformatic tools capable of parallel processing of virome sequencing data by exactly identical methods: this is especially important in studies of multifactorial diseases, or in parallel comparison of laboratory protocols.

https://sourceforge.net/projects/vipie/

Utilities / Create consensus sequence from a BAM file

szypanther — Mon, 09 Sep 2019 01:34:13 +0000

Utilities / Create consensus sequence from a BAM file

Description

Given an indexed BAM file and corresponding reference genome (in fasta format), this tool constructs a consensus sequence based on the alignment.

Parameters

None

Details

This tool uses SAMtools, bcftools and vcfutils.pl script to create a consensus sequence for the given alignment file. The actual command line executed is:

  samtools mpileup -uf reference.fa aligment.bam | bcftools view -cg - | vcfutils vcf2fq

Note that the input BAM file must be sorted before it can be used by this tool.

Output

Output is a fasta formatted sequence file.

References

This tool is based on the SAMtools package. Please cite the article The Sequence alignment/map (SAM) format and SAMtools by Li H., Handsaker B., Wysoker A., Fennell T., Ruan J., Homer N., Marth G., Abecasis G., Durbin R. and 1000 Genome Project Data Processing Subgroup (2009) Bioinformatics, 25, 2078-9. [PMID: 19505943].

Good software

szypanther — Wed, 24 Jul 2019 06:46:13 +0000

multiqc ranger ployly

https://github.com/MultiQC

https://github.com/ranger/ranger

https://zhuanlan.zhihu.com/p/34369349

Data Binning and Plotting

szypanther — Tue, 26 Feb 2019 02:27:48 +0000

In statistics, data binning is a way to categorize a number of continuous values into a smaller number of buckets (bins). Each bucket defines an numerical interval. For example, if there is a variable about house-based education levels which are measured by continuous values ranged between 0 and 19, data binning will place each value into one bucket if the value falls into the interval that the bucket covers. This post shows data binning in R as well as visualizing the bins.

The dataset contains 32038 observations for mean education level per house. Load the data into R.

data <- read.csv("data/meaneducation.csv")
x <- data[,1] #take 1st column of data frame to a vector

After loading the data, inspect the data by running

summary(x)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00   11.88   12.43   12.52   13.10   19.00

The summary data shows the range of the variable, $[0, 19] "> [0, 19]$ . The median education level per house is $12.43 "> 12.43$

As the attribute is continuous, the data distribution can be visualized in a probability distribution graph. The following script creates a density plot showing the probability distribution of x.

library(ggplot2)
ggplot(data=data, aes(x=data[,1])) +
  geom_density(fill='lightblue') +
  geom_rug() +
  labs(x='mean education per house')

While the density plot is informative, it might be too technical for less technical people to read. In this case, the better data visualization is binning the data into discrete categories and plotting the count of each bin. A sample script is provided below for creating bins and plot bin counts in a bar chart.

Creating bins

# set up boundaries for intervals/bins
breaks <- c(0,2,5,8,10,13,15,19,21)
# specify interval/bin labels
labels <- c("<2", "2-5)", "5-8)", "8-10)", "10-13)", "13-15)", "15-19)", ">=19")
# bucketing data points into bins
bins <- cut(x, breaks, include.lowest = T, right=FALSE, labels=labels)
# inspect bins
summary(bins)

    <2   2-5)   5-8)  8-10) 10-13) 13-15) 15-19)    19- 
   120      0     24    330  22595   7760   1200      9

The script above results in a factor bins with 8 levels. Each level is marked by a string in the vector labels. Each value in bins indicates the level that the observation falls into. The summary of bins lists the count for each bin.

Plotting the bins

After generating the bins, simply plot the bins by the following script:

plot(bins, main="Frequency", ylab="Number of Houses",col="bisque")

Using ggplot

y <- cbind(data, bins)
ggplot(data=y, aes(x=y$bins,fill=..count..)) + 
  geom_bar(color='black', alpha=0.9) + 
  stat_count(geom="text", aes(label=..count..), hjust=-0.1) +
  theme_bw() + 
  labs(y="Number of Houses",x="Mean Education Values") +
  coord_flip() +
  ylim(0,25000) +
  scale_x_discrete(drop=FALSE) # include the bins of length zero

ChIPseeker for ChIP peak annotation (转贴）

szypanther — Tue, 28 Nov 2017 06:35:13 +0000

https://guangchuangyu.github.io/2014/04/chipseeker-for-chip-peak-annotation/

ChIPpeakAnno WAS the only R package for ChIP peak annotation. I used it for annotating peak in my recent study.

I found it does not consider the strand information of genes. I reported the bug to the authors, but they are reluctant to change.

So I decided to develop my own package, ChIPseeker, and it’s now available in Bioconductor.

> require(TxDb.Hsapiens.UCSC.hg19.knownGene)
> txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
> require(ChIPseeker)
> peakfile = getSampleFiles()[[4]]
> peakfile
[1] "/usr/local/Cellar/r/3.1.0/R.framework/Versions/3.1/Resources/library/ChIPseeker/extdata//GEO_sample_data/GSM1295076_CBX6_BF_ChipSeq_mergedReps_peaks.bed.gz"
> peakAnno <- annotatePeak(peakfile, tssRegion=c(-3000, 1000), TxDb=txdb, annoDb="org.Hs.eg.db")
>> loading peak file...              2014-04-13 12:03:45 PM 
>> preparing features information...         2014-04-13 12:03:45 PM 
>> identifying nearest features...       2014-04-13 12:03:45 PM 
>> calculating distance from peak to TSS...  2014-04-13 12:03:46 PM 
>> assigning genomic annotation...       2014-04-13 12:03:46 PM 
>> adding gene annotation...             2014-04-13 12:03:50 PM 
>> assigning chromosome lengths          2014-04-13 12:03:50 PM 
>> done...                   2014-04-13 12:03:50 PM 
> head(as.GRanges(peakAnno))
GRanges object with 6 ranges and 14 metadata columns:
      seqnames             ranges strand |          V4        V5
                      |     
  [1]     chr1 [ 815092,  817883]      * | MACS_peak_1    295.76
  [2]     chr1 [1243287, 1244338]      * | MACS_peak_2     63.19
  [3]     chr1 [2979976, 2981228]      * | MACS_peak_3    100.16
  [4]     chr1 [3566181, 3567876]      * | MACS_peak_4    558.89
  [5]     chr1 [3816545, 3818111]      * | MACS_peak_5     57.57
  [6]     chr1 [6304864, 6305704]      * | MACS_peak_6     54.57
             annotation  geneChr geneStart   geneEnd geneLength geneStrand
                    
  [1] Distal Intergenic     chr1    803451    812182       8732          -
  [2]  Promoter (<=1kb)     chr1   1243994   1247057       3064          +
  [3]  Promoter (<=1kb)     chr1   2976181   2980350       4170          -
  [4]  Promoter (<=1kb)     chr1   3547331   3566671      19341          -
  [5]  Promoter (<=1kb)     chr1   3816968   3832011      15044          +
  [6]            3  UTR     chr1   6304252   6305638       1387          +
           geneId transcriptId distanceToTSS         ENSEMBL      SYMBOL
                   
  [1]      284593   uc001abt.4         -5701 ENSG00000230368      FAM41C
  [2]      126789   uc001aed.3             0 ENSG00000169972       PUSL1
  [3]      440556   uc001aka.3             0 ENSG00000177133   LINC00982
  [4]       49856   uc001ako.3             0 ENSG00000116213      WRAP73
  [5]   100133612   uc001alg.3             0 ENSG00000236423   LINC01134
  [6]      390992   uc009vly.2          1452 ENSG00000173673        HES3
                                          GENENAME
                                       
  [1] family with sequence similarity 41, member C
  [2]              pseudouridylate synthase-like 1
  [3]   long intergenic non-protein coding RNA 982
  [4]      WD repeat containing, antisense to TP73
  [5]  long intergenic non-protein coding RNA 1134
  [6]       hes family bHLH transcription factor 3
  ---
  seqlengths:
        chr1     chr10     chr11     chr12 ...      chr8      chr9      chrX
   249250621 135534747 135006516 133851895 ... 146364022 141213431 155270560

annotatePeak function can accept peak file or GRanges object that contained peak information. The sample peakfile shown above is a sample output from MACS. All the information contained in peakfile will be kept in the output of annotatePeak.

The annotation column annotates genomic features of the peak, that is whether peak is located in Promoter, Exon, UTR, Intron or Intergenic. If the peak is annotated by Exon or Intron, more detail information will be provided. For instance, Exon (38885 exon 3 of 11) indicates that the peak is located at the 3th exon of gene 38885 (entrezgene ID) which contain 11 exons in total.

All the names of column that contains information of nearest gene annotated will start with gene in camelCaps style.

If parameter annoDb=“org.Hs.eg.db” is provided, extra columns of ENSEMBL, SYMBOL and GENENAME will be provided in the final output.

For more detail, please refer to the package vignette. Comments or suggestions are welcome.

——- below was added at April 30, 2014 ———–

After two weeks of ChIPseeker release in Bioconductor, I am very glad to receive an email from Cambridge Systems Biology Centre that they are using ChIPseeker in annotating Drosophila ChIP data.

Although I only illustrate examples in human, ChIPseeker do supports all species only if genomic information was available.

> require(TxDb.Dmelanogaster.UCSC.dm3.ensGene)
> require(ChIPseeker)
> txdb <- TxDb.Dmelanogaster.UCSC.dm3.ensGene
>
> chr <- c("chrX", "chrX", "chrX")
> start <- c(71349, 77124, 86364)
> end <- c(72835, 77836, 88378)
> peak <- GRanges(seqnames = chr, ranges = IRanges(start, end))
> peakAnno <- annotatePeak(peak,
+                          tssRegion=c(-3000,100),
+                          TxDb = txdb,
+                          annoDb="org.Dm.eg.db")
>> preparing features information...         2014-04-30 21:28:27
>> identifying nearest features...       2014-04-30 21:28:30
>> calculating distance from peak to TSS...  2014-04-30 21:28:31
>> assigning genomic annotation...       2014-04-30 21:28:31
>> adding gene annotation...             2014-04-30 21:28:42
Loading required package: org.Dm.eg.db
Loading required package: DBI

>> assigning chromosome lengths           2014-04-30 21:28:43
>> done...                   2014-04-30 21:28:43
> as.GRanges(peakAnno)
GRanges object with 3 ranges and 12 metadata columns:
      seqnames         ranges strand |
                  |
  [1]     chrX [71349, 72835]      * |
  [2]     chrX [77124, 77836]      * |
  [3]     chrX [86364, 88378]      * |
                                           annotation  geneChr geneStart
                                            
  [1]                                Promoter (<=1kb)     chrX     72458
  [2] Intron (FBtr0100135/FBgn0029128, intron 1 of 7)     chrX     72458
  [3] Intron (FBtr0100135/FBgn0029128, intron 2 of 7)     chrX     72458
        geneEnd geneLength geneStrand      geneId transcriptId distanceToTSS
                   
  [1]     96056      23599          + FBgn0029128  FBtr0100135             0
  [2]     96056      23599          + FBgn0029128  FBtr0100135          5378
  [3]     96056      23599          + FBgn0029128  FBtr0100135         15920
         ENTREZID      SYMBOL    GENENAME
        
  [1]     3354987         tyn     trynity
  [2]     3354987         tyn     trynity
  [3]     3354987         tyn     trynity
  ---
  seqlengths:
       chrX
   22422827
>

Citation

Yu G, Wang LG and He QY^*. ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics2015, 31(14):2382-2383.

Mapping reads with bwa and bowtie

szypanther — Tue, 07 Nov 2017 09:00:20 +0000

In this tutorial, we’re going to take a set of Illumina reads from an inbred Drosophila melanogaster line, and map them back to the reference genome. (After these steps, we could do things like generate a list of SNPs at which this line differs from the reference strain, or generate a genome sequence for this fly strain, but we’ll get to that later on in the course.) We are also going to use two different (but popular) mapping tools, bwa and bowtie. Among their differences is that bowtie (while smokin’ fast) does not deal with “gapped” alignments, i.e. it does not handle insertion/deletions well.

Getting the data

If you just finished the FastQC tutorial, you can keep working on the same machine. Otherwise, launch an EC2 instance, make a volume out of the same snapshot that we did before, and mount it. The data for this tutorial are in /data/drosophila:

%% cd /data/drosophila/

For this example, you’ll also need the Drosophila reference genome:

%% curl -O ftp://ftp.flybase.net/genomes/Drosophila_melanogaster/dmel_r5.37_FB2011_05/fasta/dmel-all-chromosome-r5.37.fasta.gz
%% gunzip dmel-all-chromosome-r5.37.fasta.gz

Installing and running bwa

To actually do the mapping, we need to download and install bwa:

First we are going to grab the source files for bwa from sourceforge, using curl. It is important to know that we need to specify a few flags to let the program know that we want to leep the name (and filetype) for the file (-O) and that curl should follow relative hyperlinks (-L), to deal with redirection of the file site (or else curl won’t work with sourceforge).:

%% curl -O -L http://sourceforge.net/projects/bio-bwa/files/bwa-0.5.9.tar.bz2

Now we want to uncompress the tarball file using “tar”. x extracts, v is verbose (telling you what it is doing), f skips prompting for each individual file, and j tells it to unzip .bz2 files.:

%% tar xvfj bwa-0.5.9.tar.bz2
%% cd bwa-0.5.9

The make command calls a program that helps to automate the compiling process for the program.:

%% make

(Note: If your system doesn’t have make installed already, you’ll need to run “apt-get install make” before you can build bwa – but if you are using the right AMI, it should be pre-installed.)

Copy the executable for bwa to a directory for binaries which is in your shell search path:

%% cp bwa /usr/local/bin
%% cd ..

If you want to see what is in your shell search path (which can be modified later), run:

%% echo $PATH

Now there are several steps involved in mapping our sequence reads and getting the output into a usable form. First we need to tell bwa to make an index of the reference genome; this will take a few minutes, so we’ve already got the index already generated in the data directory, but if you want to try it yourself, you can run (but see the note below first!!):

%% bwa index dmel-all-chromosome-r5.37.fasta

Note: This step takes several minutes. If you run it, it will overwrite the index files we have already made, so don’t run the above line exactly; instead, create a copy of the reference genome and then index the copy instead, so that we can preserve our pre-computed reference index:

%% cp dmel-all-chromosome-r5.37.fasta dmel-all-chromosome-r5.37.copy.fasta
%% bwa index dmel-all-chromosome-r5.37.copy.fasta

Next, we do the actual mapping. These were paired-end reads, which means that for each DNA fragment, we have sequence data from both ends. The sequences are therefore stored in two separate files (one for the data from each end), so we have two mapping steps to perform. For now, we’ll use bwa’s default settings. The files you’ll be running this on are datasets that have been trimmed down to just the first 1 million sequence reads to speed things up, but at the end you’ll be able to work with the final product from an analysis of the full dataset that we ran earlier (some of these steps take upwards of an hour on the full dataset, but just a couple minutes on the trimmed dataset). Run:

%% bwa aln dmel-all-chromosome-r5.37.fasta RAL357_1.fastq > RAL357_1.sai
%% bwa aln dmel-all-chromosome-r5.37.fasta RAL357_2.fastq > RAL357_2.sai

These .sai files aren’t very useful to us, so we need to convert them into SAM files. In this step, bwa takes the information from the two separate ends of each sequence and combines everything together. Here’s how you do it (this may take around 10 minutes):

%% bwa sampe dmel-all-chromosome-r5.37.fasta RAL357_1.sai RAL357_2.sai RAL357_1.fastq RAL357_2.fastq > RAL357_bwa.sam

The SAM file is technically human-readable; take a look at it with:

%% more RAL357_bwa.sam

…It’s not very easy to understand (if you are really curious about the SAM format, there is a 12-page manual at http://samtools.sourceforge.net/SAM1.pdf). For now we’ll use bowtie to map the same reads, and we’ll use another tool to visualize these mappings in a more intuitive way.

bwa options

There are several options you can configure in bwa. Probably one of the most important is how many mismatches you will allow between a read and a potential mapping location for that location to be considered a match. The default is 4% of the read length, but you can set this to be either another proportion of the read length, or a fixed integer. For example, if you ran:

%% bwa aln -n 4 dmel-all-chromosome-r5.37.fasta RAL357_1.fastq > RAL357_1.sai

This would do almost the same thing as above, except this time, all locations in the reference genome that contain four or fewer mismatches to a given sequence read would be considered a match to that read.

Alternatively, you could do:

%% bwa aln -n 0.01 dmel-all-chromosome-r5.37.fasta RAL357_1.fastq > RAL357_1.sai

This would only allow reads to be mapped to locations at which the reference genome differs by 1% or less from a given read fragment.

If you want to speed things up, you can tell it to run the alignment on multiple threads at once (this will only work if your computer has a multi-core processor, which our Amazon image does). To do so, use the -t option to specify the number of threads. For example, the following line would run in two simultaneous threads:

%% bwa aln -t 2 dmel-all-chromosome-r5.37.fasta RAL357_1.fastq > RAL357_1.sai

bwa can also handle single-end reads. The only difference is that you would use samse instead of sampe to generate your SAM file:

%% bwa samse dmel-all-chromosome-r5.37.fasta RAL357_1.sai RAL357_1.fastq > RAL357_1.sam

Now let us align our reads using bowtie

(Note: For simplicity we are going to put all of the bowtie related files into the same directory. For your own work, you may want to organize your file structure better than we have).

Let’s get bowtie from Sourceforge:

%% curl -O -L http://sourceforge.net/projects/bowtie-bio/files/bowtie/0.12.7/bowtie-0.12.7-linux-x86_64.zip

unzip the file, and create a directory for bowtie. In this case, the program is precompiled so it comes as a binary executable:

%% unzip bowtie-0.12.7-linux-x86_64.zip

Change directory:

%% cd bowtie-0.12.7

Copy the bowtie files to a directory in you shell search path, and then move back to the parent directory (/data/drosophila):

%% cp bowtie bowtie-build bowtie-inspect /usr/local/bin

Let’s create a new directory, “drosophila_bowtie” where we are going to place all the bowtie results:

%% cd ..
%% mkdir drosophila_bowtie
%% cd drosophila_bowtie

Now we are going to build an index of the Drosophila genome using bowtie just like we did with bwa. The original Drosophila reference genome is in the same location as we used before. Again, we have already performed the indexing step (it takes about 7 minutes), so if you want to try it yourself, index a copy so you don’t over-write the one we’ve pre-run for you:

%%  bowtie-build /data/drosophila/dmel-all-chromosome-r5.37.fasta  drosophila_bowtie

Now we get to map! We are going to use the default options for bowtie for the moment. Let’s go through this. there are a couple of flags that we have set, since We have paired end reads for these samples, and multiple processors. The general format for bowtie is (don’t run this):

%% bowtie indexFile fastqFile outputFile

However we have some more details we want to include, so there are a couple of flags that we have to set. -S means that we want the output in SAM format. -p 2 is for multithreading (using more than one processor). In this case we have two to use. -1 -2 tells bowtie that these are paired end reads (the .fastq), and specifies which one is which.

This should take 35-40 minutes to run on the full dataset so we’ll run it on a trimmed version (should take about 3 minutes; later we’ll give you pre-computed results for the full set.):

%% bowtie -S -p 2 drosophila_bowtie -1 /data/drosophila/RAL357_1.fastq -2 /data/drosophila/RAL357_2.fastq RAL357_bowtie.sam

You may see warning messages like:

Warning: Exhausted best-first chunk memory for read SRR018286.1830915 USI-EAS034_2_PE_FC304DDAAXX:8:21:450:1640 length=45/1 (patid 1830914); skipping read

We will talk about some options you can set to deal with this.

The bowtie manual can be found here: http://bowtie-bio.sourceforge.net/manual.shtml

Some additional useful arguments/options (at least for me) -m # Suppresses all alignments for a particular read if more than m reportable alignments exist. -v # no more than v mismatches in the entire length of the read -n -l # max number of mismatches in the high quality “seed”, which is the the first l base pairs of a read. -chunkmbs # number of mb of memory a thread is given to store path. Useful when you get warnings like above –best # make Bowtie “guarantee” that reported singleton alignments are “best” given the options –tryhard # try hard to find valid alignments, when they exit. VERY SLOW.

Processing the output for use with Samtools

Even the SAM file isn’t very useful unless we can get it into a program that generates more readable output or lets us visualize things in a more intuitive way. For now, we’ll get the output into a sorted BAM file so we can look at it using Samtools later.

Download and install Samtools:

%% cd /data/drosophila
%% curl -O -L http://sourceforge.net/projects/samtools/files/samtools/0.1.16/samtools-0.1.16.tar.bz2
%% tar xvfj samtools-0.1.16.tar.bz2
%% cd samtools-0.1.16
%% make
%% cp samtools /usr/local/bin
%% cd misc/
%% cp *.pl maq2sam-long maq2sam-short md5fa md5sum-lite wgsim /usr/local/bin/
%% cd /data/drosophila

Like bwa, Samtools also requires us to go through several steps before we have our data in usable form. First, we need to have Samtools generate its own index of the reference genome:

%% samtools faidx dmel-all-chromosome-r5.37.fasta

Next, we need to convert the SAM file into a BAM file. (A BAM file is just a binary version of a SAM file.):

%% samtools import dmel-all-chromosome-r5.37.fasta.fai RAL357_bwa.sam RAL357_bwa.bam

Now, we need to sort the BAM file:

%% samtools sort RAL357_bwa.bam RAL357_bwa.sorted

And last, we need Samtools to index the BAM file:

%% samtools index RAL357_bwa.sorted.bam

Let us do this again for the bowtie output. We move back to the drosophila_bowtie directory (you could do this all from the other directory, but it gets harder to read the command with long pathnames):

%% cd drosophila_bowtie
%% samtools import ../dmel-all-chromosome-r5.37.fasta.fai RAL357_bowtie.sam RAL357_bowtie.bam

Now, we need to sort the BAM file (also slow):

%% samtools sort RAL357_bowtie.bam RAL357_bowtie.sorted

And last, we need Samtools to index the BAM file:

%% samtools index RAL357_bowtie.sorted.bam

All done! Now we can use the sorted BAM file in Samtools to visualize our mappings, generate lists of SNPs, and call consensus sequences. We’ll get to all of that later on today and in the rest of the course.

Viewing the output with TView

Before we can use TView to compare Bowtie and BWA mappings, we need to sort the Bowtie BAM file, and generate an index for it.:

%% cd drosophila_bowtie
%% samtools sort RAL357_full_bowtie.bam RAL357_full_bowtie.sorted
%% samtools index RAL357_full_bowtie.sorted.bam

Now that we’ve generated the files, we can view the output with TView. We’ll compare two different sorted:

%% cd ..
%% samtools tview RAL357_full_bwa.sorted.bam

Now open an additional terminal window, and load the Bowtie mapping file there as well.:

%% cd /data/drosophila/drosophila_bowtie
%% samtools tview RAL357_full_bowtie.sorted.bam ../dmel-all-chromosome-r5.37.fasta

BACTERIAL GENOMICS TUTORIAL （repost）

szypanther — Thu, 08 Jun 2017 04:00:22 +0000

[Originally posted by Kat on her BacPathGenomics blog, April 2013]

This is a shameless plug for an article and accompanying tutorial I’ve just published together with David Edwards, my excellent MSc Bioinformatics student from the University of Melbourne. It’s currently available as a PDF pre-pub from BMC Microbial Informatics and Experimentation, but the web version will be available soon. The accompanying tutorial is available here.

The idea for this came from discussions at last year’s ASM (Australian Society of Microbiology) meeting, where it was highlighted that there was a lack of courses and tutorials available for biologists to learn the basics of genomic analysis so that they can make use of next gen sequencing. Michael Wise, a founding editor of BMC Microbial Informatics and Experimentation based at UWA in Perth, suggested the new journal would be an ideal home for such a tutorial… so here we are:

Beginner’s guide to comparative bacterial genome analysis using next-generation sequence data

http://www.microbialinformaticsj.com/content/3/1/2/

High throughput sequencing is now fast and cheap enough to be considered part of the toolbox for investigating bacteria, and there are thousands of bacterial genome sequences available for comparison in the public domain. Bacterial genome analysis is increasingly being performed by diverse groups in research, clinical and public health labs alike, who are interested in a wide array of topics related to bacterial genetics and evolution. Examples include outbreak analysis and the study of pathogenicity and antimicrobial resistance. In this beginner’s guide, we aim to provide an entry point for individuals with a biology background who want to perform their own bioinformatics analysis of bacterial genome data, to enable them to answer their own research questions. We assume readers will be familiar with genetics and the basic nature of sequence data, but do not assume any computer programming skills. The main topics covered are assembly, ordering of contigs, annotation, genome comparison and extracting common typing information. Each section includes worked examples using publicly available E. coli data and free software tools, all which can be performed on a desktop computer.

Four great tools

In the paper and tutorial, we introduce the four tools which we rely on most for basic analysis of bacterial genome assemblies: Velvet, ACT, Mauve and BRIG. All except ACT were developed as part of a PhD project, and have endured well beyond the original PhD to become well-known bioinformatics tools. New students take note!

In the paper, each tool is highlighted in its own figure, which includes some basic instructions. This is reproduced below, but is covered in much more detail in the tutorial that comes with the paper (link at the bottom).

1. Velvet for genome assembly

Possibly the most popular and widely used short read assembler, developed by the amazing Dan Zerbino during his PhD at EBI in Cambridge. Quite a PhD project!

[ Download | Paper | Protocol ]

Reads are assembled into contigs using Velvet and VelvetOptimiser in two steps, (1) velveth converts reads to k-mers using a hash table, and (2) velvetg assembles overlapping k-mers into contigs via a de Bruijn graph. VelvetOptimiser can be used to automate the optimisation of parameters for velveth and velvetg and generate an optimal assembly. To generate an assembly of E. coli O104:H4 using the command-line toolVelvet:

• Download Velvet [23] (we used version 1.2.08 on Mac OS X, compiled with a maximum k-mer length of 101 bp)

• Download the paired-end Illumina reads for E. coli O104:H4 strain TY-2482 (ENA accession SRR292770)

• Convert the reads to k-mers using this command:

velveth out_data_35 35 -fastq.gz -shortPaired -separate SRR292770_1.fastq.gz SRR292770_2.fastq.gz

• Then, assemble overlapping k-mers into contigs using this command:

velvetg out_data_35 -clean yes -exp_cov 21 -cov_cutoff 2.81 -min_contig_lgth 200

This will produce a set of contigs in multifasta format for further analysis. See Additional file 1: Tutorial for further details, including help with downloading reads and using VelvetOptimiser.

2. ACT for pairwise genome comparison

Part of the Sanger Institute’s Artemis suite of tools. Also look at Artemis (single genome viewer), DNA Plotter (which can draw circular diagrams of your genomes) and BAMView (which can display mapped reads overlaid on a reference genome), they are all available here.

[ Download | Paper | Manual ]

Artemis and ACT are free, interactive genome browsers (we used ACT 11.0.0 on Mac OS X).

• Open the assembled E. coli O104:H4 contigs in Artemis and write out a single, concatenated sequence using File -> Write -> All Bases -> FASTA Format.

• Generate a comparison file between the concatenated contigs and 2 alternative reference genomes using the website WebACT.

• Launch ACT and load in the reference sequences, contigs and comparison files, to get a 3-way comparison like the one shown here.

Here, the E. coli O104:H4 contigs are in the middle row, the enteroaggregative E. coli strain Ec55989 is on top and the enterohaemorrhagic E. coli strain EDL933 is below. Details of the comparison can be viewed by zooming in, to the level of genes or DNA bases.

3. Mauve for contig ordering and multiple genome comparison

Developed by the wonderful Aaron Darling during his PhD, he is now Associate Professor at University of Technology Sydney. Also see Mauve Assembly Metrics, an optional plugin for assessing assembly quality which was developed for the Assemblathon.

[ Download | Paper | User Guide ]

Mauve is a free alignment tool with an interactive browser for visualising results (we used Mauve 2.3.1 on Mac OS X).

• Launch Mauve and select File -> Align with progressiveMauve

• Click ‘Add Sequence…’ to add your genome assembly (e.g. annotated E. coli O104:H4 contigs) and other reference genomes for comparison.

• Specify a file for output, then click ‘Align…’

• When the alignment is finished, a visualization of the genome blocks and their homology will be displayed, as shown here. E. coli O104:H4 is on the top, red lines indicate contig boundaries within the assembly. Sequences outside coloured blocks do not have homologs in the other genomes.

4. BRIG (BLAST Ring Image Generator) for multiple genome comparison

From Nabil-Fareed Alikhan at the University of Queensland, also as part of a graduate project, which I believe is still in progress…

[ Download | Download BLAST | Paper | Tutorial ]

BRIG is a free tool that requires a local installation of BLAST (we used BRIG 0.95 on Mac OS X). The output is a static image.

• Launch BRIG and set the reference sequence (EHEC EDL933 chromosome) and the location of other E. coli sequences for comparison. If you include reference sequences for the Stx2 phage and LEE pathogenicity island, it will be easy to see where these sequences are located.

• Click ‘Next’ and specify the sequence data and colour for each ring to be displayed in comparison to the reference.

• Click ‘Next’ and specify a title for the centre of the image and an output file, then click ‘Submit’ to run BRIG.

• BRIG will create an output file containing a circular image like the one shown here. It is easy to see that the Stx2 phage is present in the EHEC chromosomes (purple) and the outbreak genome (black), but not the EAEC or EPEC chromosomes.

Tutorial

The tutorial accompanying the article is available here. To give you an idea of what’s covered, here is the table of contents:

1. Genome assembly and annotation…………………………………………………………… 2

1.1 Downloading E. coli sequences for assembly…………………………………………….. 2

1.2 Examining quality of reads (FastQC)………………………………………………………… 2

1.3 Velvet – assembling reads into contigs………………………………………………………. 4

1.3.1 Using VelvetOptimiser to optimise de novo assembly with Velvet………….. 6

1.4 Ordering contigs against a reference using Mauve………………………………………. 7

1.4.1 Viewing the ordered contigs (Mauve)………………………………………………… 10

1.4.2 Viewing the ordered contigs (ACT)……………………………………………………. 13

1.5 Mauve Assembly Metrics – Statistical View of the Contigs………………………… 15

1.6 Annotation with RAST……………………………………………………………………………. 15

1.6.1 Alternatives to RAST………………………………………………………………………. 19

2. Comparative genome analysis……………………………………………………………….. 20

2.1 Downloading E. coli genome sequences for comparative analysis………………. 20

2.2 Mauve – for multiple genome alignment……………………………………………………. 21

2.3 ACT – for detailed pairwise genome comparisons……………………………………… 24

2.3.1 Generating comparison files for ACT…………………………………………………. 24

2.3.2 Viewing genome comparisons in ACT……………………………………………….. 27

2.4 BRIG – Visualizing reference-based comparisons of multiple sequences……… 29

3. Typing and specialist tools……………………………………………………………………. 34

3.1 PHAST – for identification of phage sequences…………………………………………. 34

3.2 ResFinder – for identification of resistance gene sequences………………………… 34

3.3 Multilocus sequence typing…………………………………………………………………….. 34

3.4 PATRIC – online genome comparison tool………………………………………………… 34

POPULATION GENOMICS OF KLEBSIELLA (Repost)

szypanther — Thu, 08 Jun 2017 03:58:56 +0000

https://holtlab.net/2015/06/23/population-genomics-of-klebsiella/

Well, after almost 6 years, our Klebsiella pneumoniae genomics paper is finally out!

It’s a beast of a thing and there are still a million and one questions to address just from this one data set. For those interested in looking at the data for themselves, the raw reads are available under accessionERP000165, the assemblies are in Sylvain Brisse’s Klebsiella pneumoniae BIGSdb at the Pasteur Institute, and the tree + metadata are available for your interactive viewing pleasure in MicroReact.

The paper itself is open access in PNAS, you can read it here.

Genomic analysis of diversity, population structure, virulence, and antimicrobial resistance in Klebsiella pneumoniae, an urgent threat to public health

There have been lots of really nice Klebs genomics papers out in the last 18 months or so, examining the evolution of the ST258 clone that carries the KPC gene (K. pneumoniae carbapaenemase) and is wreaking havoc in hospitals all over the place (including recently in Melbourne), and also several hospital-based studies tracking transmission and evolution of local drug-resistant outbreaks.

But that is just the tip of the K. pneumoniae iceberg.

Our paper asks a completely different set of questions, which you could basically sum up as “what the hell is Klebsiella pneumoniae anyway?”

To do this, we sequenced ~300 genomes of really diverse K. pneumoniae strains. We didn’t have much information about genetic diversity to go on, so we chose strains with different phenotypes (antimicrobial resistance patterns, capsular serotypes or sequence types where known), from different sources (human and animal, asymptomatic carriage and infections of various kinds), and from different geographical locations.

This was done by an international group of collaborators who pooled their resources, not only sharing their precious strain collections but also digging through hospital and other records to find as much information about the strains as possible.

You can view the tree and associated metadata, including geographical origin and source information, over on Microreact.

We found out some pretty interesting things about Klebsiella pneumoniae, including the fact that what’s identified as K. pneumoniae using standard tests is actually a mixed bag of three related species, that now have their own names: K. pneumoniae (KpI group, which includes the majority of clinical isolates and all the stuff you might have heard of like the clone that causes rhinoscleromatis, and the KPC clone ST258, and the hypervirulent clone ST23); K. quasipneumoniae; and K. variicola (plant associated and usually nitrogen-fixing).

By now, this species stuff has been nutted out (mainly by co-author Sylvain Brisse from Institut Pasteur) by analysing marker gene sequences, but it’s really important to be able to show that those patterns hold at the whole-genome level, and we found some interesting things about the distribution of the rarer species (see the paper for details).

Importantly, we did the whole pan-genome analysis thing and found that as a population, K. pneumoniaehas more genes than humans. Almost 30,000 in fact. Each individual strain has ~5,500 genes, but <2,000 of those are core genes that are common to all K. pneumoniae. The rest are accessory genes that can come and go, helping the bug to adapt to new environments.

One of the cool things we were able to do with our data set, which you just can’t do with genomic studies focused on specific clones or outbreaks, was to look at statistical associations between accessory genes and phenotypes. Admittedly our available phenotypes were pretty limited, but we found a few important things.

1. VIRULENCE

We screened for genes associated with virulence in humans by focusing in on invasive infections, and comparing gene frequencies in human isolates from invasive community-acquired infections (i.e. the kind of infections that land you in hospital) vs. those in human carriage isolates or hospital acquired infections (i.e. the kind of infections that get you when you are already in hospital for something else and are particularly vulnerable to infection).

The only genes that were significantly associated with invasive infection in humans were rmpA and rmpA2, which upregulate capsule production, and genes related to iron acquisition (specifically acquired siderophore systems that can help to steal iron from animal hosts – see paper for details). These genes have been known about for some time, based on mouse models and knowledge of other pathogens, however we were able to show that these genes are significantly associated with invasive K. pneumoniaedisease in humans, which is not something that can be proven directly using experimental systems. (The siderophore story actually goes a bit deeper than the iron issue… it’s a bit too complex to go into here but I recommend reading Michael Bachman’s work e.g. “Interaction of lipocalin 2, transferrin, and siderophores determines the replicative niche of Klebsiella pneumoniae during pneumonia” in MBio, 2012).

Interestingly, doing the same test in bovine isolates showed that the story is very different: we had a lot of isolates from dairy herds, including clinical and subclinical mastitis; asymptomatic carriage isolates and strains from the farm environment… and found that an acquired lactose operon was almost perfectly associated with mastitis in cows! Something similar has been observed before in Streptococcus agalactiae.

2. ANTIBIOTIC RESISTANCE

Resistance genes were associated with human hospital isolates and human carriage isolates. This is far from an ideal study design to test this, as we had different types of collections from different geographical regions; however, even when you look within different local collections you see the same patterns: (a) comparing bovine and human isolates from NY state, the resistance genes were all in human isolates not cow isolates; (b) comparing human carriage and infection isolates (both nosocomial and community acquired) in Vietnam, the resistance genes were mainly in human carriage and hospital isolates, not in community infections; (c) in the remaining countries, isolates from infections acquired in hospital had more resistance genes than those that were considered nosocomial (diagnosed within 48 hours of admission).

What’s really interesting is that while resistance genes and virulence genes are both highly mobile components of the accessory genome, they were essentially orthogonal in their distribution. The resistance genes were mainly in hospital acquired infections and carriage isolates, whereas the virulence strains were mainly found in isolates from community acquired infections.

So far, this has resulted in the emergence of two very different kinds of K. pneumoniae clones of importance to human health: hypervirulent clones, and multidrug resistant clones. This is pretty lucky, as it means the hypervirulent clones are generally sensitive to antibiotics (although antimicrobial treatment is difficult for some conditions, like liver abscess), and the problem of untreatable highly drug resistant Klebs infections has not spread outside of hospitals.

Unfortunately, our luck appears to be runnning out and we are already starting to see the convergence of virulence and resistance. Hypervirulent ST23 strains, which have all four of the acquired siderophore systems, are accumulating antibiotic resistance genes. And about half of the KPC Klebs ST258 strains causing problems in hospitals globally have one of the siderophore gene clusters, yersiniabactin, which has been shown in clinical ST258 isolates to confer enhanced ability to cause pneumonia. How long till the other virulence genes creep in? We need to be watching!

Also, our data indicates that there are plenty of other hypervirulent or multidrug resistant Klebsclones emerging out there… convergence of virulence and resistance could happen in any one of them, so we need to be thinking and monitoring beyond the well-known ST23 and ST258 strains.

In any case, genomic surveillance is going to become really important for Klebsiella…

TOOLS FOR BACTERIAL COMPARATIVE GENOMICS

szypanther — Thu, 08 Jun 2017 03:51:17 +0000

Yesterday I spoke at a workshop for JAMS TOAST (Sydney’s Joint Academic Microbiology Seminars – bioinformatics workshop)… I was asked to cover tools for comparative genomics, so I put together a list of the tried and tested programs that I find most useful for this kind of analysis. So here is the list.

First, a few caveats…

These are mostly tools with a graphical user interface (mostly Java based)… this means they should be pretty accessible to most users, however if you want to do analyses that are a bit more custom or niche, you will have to get your hands dirty and use the commandline (which you should learn to do anyway!!)

These tools are useful for small-ish scale genomic comparisons, in the order of 2-20 genomes.

Most of these tools are for assembled data, hence we start with how to assemble your data… this will become less of an issue as we move to long read sequencing with PacBio and MinION etc, but for the moment most of the data I work with is from large scale sequencing projects with Illumina (100s-1000s) so we use mapping-based approaches for a lot of tasks… so I have included a few comments about this at the end.

Beginner’s guide with walk-through tutorial

Some of these tools, particularly the visualisation of whole genome comparisons (using Artemis & ACT, Mauve, and BRIG) are covered the in the tutorial from our 2013 “Beginner’s guide to comparative bacterial genome analysis using next-generation sequence data“. So if you want a walk-through, that’s a good place to start. Please note that the links to the read sets used in the tutorial have changed, but you can still find the reads at NCBI or ENA under the accession provided (SRR292770).

First things first – Are my reads good quality?

FastQC – Generate graphical reports of read quality from the fastq files.

Assembly

SPAdes – de Bruijn graph assembly, incorporating multiple kmers and read pairing information in the building of the graph. Think of this as a more sophisticated version of Velvet… in my experience, it nearly always provides better assemblies than Velvet, except on the rare occasion (1-5% of read sets) where it fails to get a good assembly at all. In which case, try Velvet!

Velvet – The first and most widely used de Bruijn graph assembler built to tackle the problem of short reads. Graphs are built using a single kmer value, and read pairing information used for scaffolding only (unlike SPAdes, where multiple kmers are incorporated into a single graph and read pairing is also used directly in building the graph). How do you know what kmer to use? Use Velvet Optimiser. Hate the command line? Try Vague, a GUI wrapper for Velvet.

How do I judge if I have a good assembly? Try QUAST

What other assemblers are there? What’s best for what task? Take a look at Nucelotid.es and Assemblathon.

How can I view my assembly graphs? Try Bandage – freshly released from Ryan Wick, a MSc (Bioinformatics) student in my lab. Bandage allows you to view and manipulate de Druijn graphs output by Velvet or SPAdes… lots of super cool features and useful applications, see the github site for examples.

Working with assembled data

Now you have a nice set of assembled contigs – where are all the genes?

Whole genome annotation

RAST – Web tool (upload contigs), uses the subsystems in the SEED database and provides detailed annotation and pathway analysis. Takes several hours per genome but I think this is the best way to get a high quality annotation (if you have only a few genomes to annotate).

Prokka – Standalone command line tool, takes just a few minutes per genome. This is the best way to get good quality annotation in a flash, which is particularly useful if you have loads of genomes or need to annotate a pangenome or metagenome. Note however that the quality of functional information is not as good as RAST, and you will need several extra steps if you want to do functional profiling and pathway analysis of your genome(s)… which is in-built in RAST.

Annotating specific types of features

Resistance genes

CARD – best combination of easy interface + pretty good database
ARG-Annot – best quality database (in my experience, focusing on Enterobacteriaceae)
ResFinder – easy interface, database needs ongoing development

Virulence genes

PATRIC – for certain bugs only, but has good online tools for genome comparisons.
VFDB – broader range of species, but varying levels of comprehensiveness and you need to do more of the work yourself.

Insertion sequences

IS saga – Upload your genome and have IS saga find all the transposes in your genome using their IS finder database

Phage

PHAST – Upload your genome and this will identify likely prophage regions, summarising these at the level of whole phage and also individual genes.

Viewing your genome – The Artemis Genome Browser

There are zillions of genome browsers out there, but I still love Artemis… and not just because I’m from the Sanger Institute. Unlike most genome browsers, Artemis was custom-built for bacterial genomes, which let’s face it are really quite different from humans and other eukaryotes.

The default view shows you your sequence and annotation, with 6 frame translation and allows you to easily edit or create features in the annotation, graph sequence-based functions like GC content and GC skew, and do all manner of other useful things. It’s been around for a zillion years (well, at least 10 or so) and is very well developed and supported.

Artemis has lots of cool features built in, including the ‘BamView’ feature that allows you to view BAM files that show the alignment of reads mapped to your genome, zoomed in to the base level or zoomed out to look at coverage and SNP distributions… this is also super handy for viewing RNAseq data, as you can easily see the stacks of reads derived from coding regions.

Artemis also has DNA Plotter built in, which you can use to generate those pretty circular figures of your genome sequences and their features.

Plus, when you’ve got used to using Artemis to get to know your shiny new genome, you can move on toviewing comparisons against other genomes using ACT – the Artemis Comparison Tool.

Comparing whole genome assemblies

NOTE: Walk-throughs of these tools, using examples from the 2011 E. coli outbreak in Germany, are covered in the “Beginner’s guide to comparative bacterial genome analysis using next-generation sequence data“.

ACT (Artemis Comparison Tool) – Visualises BLAST (or similar) comparisons of genomes. This is most useful for comparisons of two or a few genomes, and makes it easy to spot and zoom in to regions of difference.

Mauve – Whole genome alignment and viewer that can output SNPs, regions of difference, homologous blocks, etc. It can also be used to assess assembly quality against a reference, using Mauve Contig Metrics.

BRIG (BLAST Ring Image Generator) – Gives a global view of whole genome comparisons by visualising BLAST comparisons via pretty circular figures. This is suitable for comparing lots of genomes, although because you have to enter each one through the GUI, it’s tricky to do more than a dozen or so.

Whole genome SNP-based phylogenies (from assembled data)

You can’t go past Adam Phiippy’s Harvest Suite

Parsnp – Compare genomes to a reference (using MUMmer) to identify core genome SNPs and build a phylogeny

Gingr – View the phylogeny and associated SNP calls (VCF format)… also useful for visualising tree + VCF that you have created in other ways, e.g. from mapping.

Detecting recombination in whole genome comparisons

Gubbins – A new implementation of the approach first used in Nick Croucher’s 2011 Science paper on Streptococcus pneumoniae. Command-line driven and runs pretty fast (<2 hours usually on our data).

BRAT NextGen – Uses a similar idea to Gubbins but using Bayesian clustering is GUI-driven… sounds nice, but actually I find it less convenient than Gubbins as there are manual steps required and then you need to run lots of iterations to get significance values.

Mapping based analyses

Why?

If you have specific questions to answer, where precise variant detection is important (e.g. allele calling, MLST, SNP detection, typing, mutation detection), mapping provides greater sensitivity and specificity than assembled data. Basically, if you want to be really sure about a variant call, you should be using the full information available in the reads rather than relying on the assembler and consensus base caller to get things right every time. See our SRST2 paper if you don’t believe me.

Also, if you need quick answers to specific questions, this is almost always going to be achieved faster and more accurately if you work direct from reads without attempting to generate high quality assemblies first.

The basics

For mapping our go-to is BWA or Bowtie2 (getting from fastq -> BAM). For processing of BAMs we use: SAMtools and BAMtools for variant calling, and BAMstats and BEDtools for summarising coverage and other information from the alignments.

Pipelines for specific tasks

There are loads of pipelines around the place that use the basic tools above to do specific tasks. A few of ours are:

SRST2 – MLST, resistance genes, virulence genes
ISMapper – IS (insertion sequence / tranposase) insertions
RedDog – Whole genome SNP-based phylogenies

https://holtlab.net/2015/02/25/tools-for-bacterial-comparative-genomics/