A de novo assembly with spades should work as well as a mapping with bowtie2.
While performing a de novo assembly you should also be aware of possible contamination which you may need to remove. By using a mapping this may not be necessary.
When you use a mapping, you must be aware that you may miss relevant and diverged regions because of the mapping to a reference strain. So a mappping would be preferrable if your target genome is closely related to your reference genome.
I'd like to give you the advice to use some kind of distribution system like [Anaconda][1] which helps you to keep your system clean and make it easier to reproduce your steps.
Here are some possible steps for both ways. I am using an e.coli k12 sample, sra number SRR6940089 :
# De novo assembly #
## 1. Quality Check ( you already did this ) ##
## 2. Run Spades ##
spades.py -o k12 -m 8 -t 4 --pe1-1 SRR6940089_1.fastq --pe1-2 SRR6940089_2.fastq &> k12.log
This command will start an assembly using 4 threads and 8 GB RAM. You can scale it to whatever your system is capable of or what you'd like.
You now have a folder called `k12` which holds the assembly information.
To get some information about the contigs you can use the [emboss][2] toolkit.
infoseq k12/contigs.fasta
fasta::contigs.fasta:NODE_859_length_78_cov_1.000000 - NODE_859_length_78_cov_1.000000 - N 78 43.59
fasta::contigs.fasta:NODE_860_length_78_cov_1.000000 - NODE_860_length_78_cov_1.000000 - N 78 44.87
fasta::contigs.fasta:NODE_861_length_78_cov_1.000000 - NODE_861_length_78_cov_1.000000 - N 78 43.59
This is the end of the contigs.fasta. So you see contigs.fasta contains 861 contigs. And you also see thate ther length and coverage is quiet low. So you need to remove this low coverage regions.
## 3. Remove low quality contigs ##
There are multiple ways how you could achieve this. I'd like to filter sequence smaller than 500 bp. the tricky part is to identify the coverage threshold. Contigs with a low coverage could be some artifacts or contamination.
## 4. Check for possible contamination ##
So to check for possible contamination you can use [kraken2][3]. But it needs some disk space. I am using here the minikraken database
kraken2 --db minikraken --threads 4 --output ecoli.kraken --report ecoli_kraken.report k12/contigs.fasta &> kraken.log
The `ecoli_kraken.report` shows an overview of how the contigs could be classified to availible species from the database.
Based on this results you can also remove some contigs.
## 5. Run Prokka
Now after you removed unwanted contigs, you can run [prokka][4] to get the annotation.
prokka --cpus 4 --outdir prokka k12/contigs.fasta &> prokka.log
And you will find the information in the prokka folder
# Mapping #
The mapping process is a little bit easier. You can use [bowtie2][5] and [samtools][6]
bowtie2 -x ecoli -1 SRR6940089_1.fastq -2 SRR6940089_2.fastq -p 4 | samtools view -@ 4 -O BAM -o ecoli.bam
158216 reads; of these:
158216 (100.00%) were paired; of these:
35696 (22.56%) aligned concordantly 0 times
120047 (75.88%) aligned concordantly exactly 1 time
2473 (1.56%) aligned concordantly >1 times
----
35696 pairs aligned concordantly 0 times; of these:
4458 (12.49%) aligned discordantly 1 time
----
31238 pairs aligned 0 times concordantly or discordantly; of these:
62476 mates make up the pairs; of these:
56512 (90.45%) aligned 0 times
5365 (8.59%) aligned exactly 1 time
599 (0.96%) aligned >1 times
82.14% overall alignment rate
You now also get some information about the similarity to your reference sequence.
Now you need to work a little bit more on the BAM file. You need to sort and index this.
samtools sort ecoli.bam -o ecoli_sort.bam
samtools index ecoli_sort.bam
Now you can create a consensus sequence using [bcftools][7] and seqtk
bcftools mpileup -f sequence.fasta ecoli_sorted.bam | bcftools call -c | vcfutils.pl vcf2fq | seqtk seq -aQ64 -q 20 -n N > ecoli_consensus.fasta
This gives you now a sequence where unknown bases are marked as N.
As I wrote before, using anaconda would make the installation a lot easier. So you can write the following lines into the file `bio.yaml`
name: bio
channels:
- bioconda
- defaults
dependencies:
- aragorn=1.2.38=h470a237_2
- barrnap=0.9=2
- bcftools=1.9=h4da6232_0
- bedtools=2.27.1=he941832_2
- blast=2.6.0=boost1.60_0
- bowtie2=2.3.4.3=py36h2d50403_0
- emboss=6.6.0=h6debe1e_0
- hmmer=3.2.1=hfc679d8_0
- infernal=1.1.2=h470a237_1
- kraken2=2.0.7_beta=pl526h2d50403_0
- libdeflate=1.0=h470a237_0
- libidn=7.45.0=2
- minced=0.3.2=0
- parallel=20160622=1
- perl-app-cpanminus=1.7044=pl526_1
- perl-bioperl=1.6.924=4
- perl-carp=1.38=pl526_1
- perl-constant=1.33=pl526_1
- perl-data-dumper=2.161=pl526_2
- perl-encode=2.88=pl526_1
- perl-exporter=5.72=pl526_1
- perl-extutils-makemaker=7.34=pl526_2
- perl-file-path=2.15=pl526_0
- perl-file-temp=0.2304=pl526_2
- perl-parent=0.236=pl526_1
- perl-threaded=5.22.0=13
- perl-xml-namespacesupport=1.12=pl526_0
- perl-xml-parser=2.44=pl526h3a4f0e9_6
- perl-xml-sax=1.00=pl526_0
- perl-xml-sax-base=1.09=pl526_0
- perl-xml-sax-expat=0.51=pl526_2
- perl-xml-simple=2.25=pl526_0
- perl-yaml=1.27=pl526_0
- prodigal=2.6.3=1
- prokka=1.12=3
- samtools=1.9=h8ee4bcc_1
- spades=3.12.0=1
- tbl2asn=25.3=0
- boost=1.60.0=py36_0
- bzip2=1.0.6=h14c3975_5
- ca-certificates=2018.03.07=0
- certifi=2018.10.15=py36_0
- curl=7.62.0=hbc83047_0
- expat=2.2.6=he6710b0_0
- fontconfig=2.13.0=h9420a91_0
- freetype=2.9.1=h8a8886c_1
- giflib=5.1.4=h14c3975_1
- icu=58.2=h9c2bf20_1
- jpeg=9b=h024ee3a_2
- libcurl=7.62.0=h20c2e04_0
- libedit=3.1.20170329=h6b74fdf_2
- libffi=3.2.1=hd88cf55_4
- libgcc=7.2.0=h69d50b8_2
- libgcc-ng=8.2.0=hdf63c60_1
- libgd=2.2.5=hceca4fd_3
- libpng=1.6.35=hbc83047_0
- libssh2=1.8.0=h1ba5d50_4
- libstdcxx-ng=8.2.0=hdf63c60_1
- libtiff=4.0.9=he85c1e1_2
- libuuid=1.0.3=h1bed415_2
- libwebp=1.0.0=h222930b_1
- libxml2=2.9.8=h26e45fe_1
- ncurses=6.1=hf484d3e_0
- openjdk=8.0.152=h46b5887_1
- openssl=1.1.1=h7b6447c_0
- perl=5.26.2=h14c3975_0
- pip=18.1=py36_0
- python=3.6.7=h0371630_0
- readline=7.0=h7b6447c_5
- setuptools=40.5.0=py36_0
- sqlite=3.25.2=h7b6447c_0
- tk=8.6.8=hbc83047_0
- wheel=0.32.2=py36_0
- xz=5.2.4=h14c3975_4
- zlib=1.2.11=ha838bed_2
Now you can run following command and most of the installation part is completed. Just the database for kraken is missing.
conda env create -f bio.yaml
I hope this gives you some help how to process your data. If you need some information just let me know ;)
[1]: https://anaconda.org/
[2]: http://emboss.sourceforge.net/
[3]: https://ccb.jhu.edu/software/kraken2/
[4]: http://www.vicbioinformatics.com/software.prokka.shtml
[5]: http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml
[6]: http://www.htslib.org/
[7]: https://samtools.github.io/bcftools/bcftools.html