<p>A de novo assembly with spades should work as well as a mapping with bowtie2.</p>
<p>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.</p>
<p>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.</p>
<p>I'd like to give you the advice to use some kind of distribution system like <a href="https://anaconda.org/" rel="nofollow noreferrer">Anaconda</a> which helps you to keep your system clean and make it easier to reproduce your steps.</p>
<p>Here are some possible steps for both ways. I am using an e.coli k12 sample, sra number SRR6940089 :</p>
<h1>De novo assembly</h1>
<h2>1. Quality Check ( you already did this )</h2>
<h2>2. Run Spades</h2>
<pre><code>spades.py -o k12 -m 8 -t 4 --pe1-1 SRR6940089_1.fastq --pe1-2 SRR6940089_2.fastq &amp;&gt; k12.log
</code></pre>
<p>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 <code>k12</code> which holds the assembly information.
To get some information about the contigs you can use the <a href="http://emboss.sourceforge.net/" rel="nofollow noreferrer">emboss</a> toolkit.</p>
<pre><code>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  
</code></pre>
<p>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.</p>
<h2>3. Remove low quality contigs</h2>
<p>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.</p>
<h2>4. Check for possible contamination</h2>
<p>So to check for possible contamination you can use <a href="https://ccb.jhu.edu/software/kraken2/" rel="nofollow noreferrer">kraken2</a>. But it needs some disk space. I am using here the minikraken database</p>
<pre><code>kraken2 --db minikraken --threads 4 --output ecoli.kraken --report ecoli_kraken.report k12/contigs.fasta &amp;&gt; kraken.log
</code></pre>
<p>The <code>ecoli_kraken.report</code> shows an overview of how the contigs could be classified to availible species from the database.</p>
<p>Based on this results you can also remove some contigs.</p>
<h2>5. Run Prokka</h2>
<p>Now after you removed unwanted contigs, you can run <a href="http://www.vicbioinformatics.com/software.prokka.shtml" rel="nofollow noreferrer">prokka</a> to get the annotation.</p>
<pre><code>prokka --cpus 4 --outdir prokka k12/contigs.fasta &amp;&gt; prokka.log
</code></pre>
<p>And you will find the information in the prokka folder</p>
<h1>Mapping</h1>
<p>The mapping process is a little bit easier. You can use <a href="http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml" rel="nofollow noreferrer">bowtie2</a> and <a href="http://www.htslib.org/" rel="nofollow noreferrer">samtools</a></p>
<pre><code>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 &gt;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 &gt;1 times
82.14% overall alignment rate
</code></pre>
<p>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.</p>
<pre><code>samtools sort ecoli.bam -o ecoli_sort.bam
samtools index ecoli_sort.bam
</code></pre>
<p>Now you can create a consensus sequence using <a href="https://samtools.github.io/bcftools/bcftools.html" rel="nofollow noreferrer">bcftools</a> and seqtk</p>
<pre><code>bcftools mpileup -f sequence.fasta ecoli_sorted.bam | bcftools call -c | vcfutils.pl vcf2fq | seqtk seq -aQ64 -q 20 -n N &gt; ecoli_consensus.fasta
</code></pre>
<p>This gives you now a sequence where unknown bases are marked as N.</p>
<p>As I wrote before, using anaconda would make the installation a lot easier. So you can write the following lines into the file <code>bio.yaml</code></p>
<pre><code>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
</code></pre>
<p>Now you can run following command and most of the installation part is completed. Just the database for kraken is missing.</p>
<pre><code>conda env create -f bio.yaml
</code></pre>
<p>I hope this gives you some help how to process your data. If you need some information just let me know ;)</p>

CC BY-SA 4.0
2018-11-14T21:50:58.357Z
2018-11-14T21:50:58.357Z
2020-06-18T08:30:41.743Z
Mr_Z
<p>I have been given the task of assembling a 'new' Ecoli genome and analysing the genes present etc.</p>

<p>The Ecoli is a new strain, and has been taken and run on a Nextseq 500 in high-output mode with 150bp paired end reads. The 'raw' files that I have is the forward and reverse reads.</p>

<p>I have initially QC checked the 'raw' files, and subsequently run them through trim galore and checked the QC after that.</p>

<p>For the next step, I now need to assemble my genome. I have been told that SPades will run a 'de novo' assembly for me, and then put that assembly into Prokka for Gene annotation.</p>

<p>Is this the best way to assemble the genome and annotate it? Or should I use another method? I am thinking that I should use a 'mapping' technique to assemble the genome using the Ecoli O157:H7 genome as a reference, but I have no idea how to do this. I would say that I am at an intermediate level with unix, but by no means am I a bioinformatician. Some help and guidance would be greatly appreciated!</p>

<p>Thanks,</p>

<p>Alex.</p>

CC BY-SA 4.0
2018-11-12T19:31:18.510Z
You can first do assembly of contigs with spades and then use a reference to get longer sequences, for example you can try this: https://github.com/fenderglass/Ragout
Peter Menzel
<ngs><assembly>
E.coli Sequencing & Analysis
# E.coli Sequencing & Analysis

(Asked by: username_0 on 2018-11-12)

I have been given the task of assembling a 'new' Ecoli genome and analysing the genes present etc.

The Ecoli is a new strain, and has been taken and run on a Nextseq 500 in high-output mode with 150bp paired end reads. The 'raw' files that I have is the forward and reverse reads.

I have initially QC checked the 'raw' files, and subsequently run them through trim galore and checked the QC after that.

For the next step, I now need to assemble my genome. I have been told that SPades will run a 'de novo' assembly for me, and then put that assembly into Prokka for Gene annotation.

Is this the best way to assemble the genome and annotate it? Or should I use another method? I am thinking that I should use a 'mapping' technique to assemble the genome using the Ecoli O157:H7 genome as a reference, but I have no idea how to do this. I would say that I am at an intermediate level with unix, but by no means am I a bioinformatician. Some help and guidance would be greatly appreciated!

Thanks,

Alex.

**Tags:** ngs, assembly
> **Comments:**
> **username_1**: You can first do assembly of contigs with spades and then use a reference to get longer sequences, for example you can try this: https://github.com/fenderglass/Ragout  

---

## 1 answers

### Answer by username_2 on 2018-11-14. Score: 1

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 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 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. 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 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 and samtools
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 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 ;)


---

