text: user3668: 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.
  
  user3360: 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 ;)
  
original_text: 3668: <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>
   3360: <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>
  