Paired-end support

[mictadlo]

Hi,
Does SGA-ICE support paired-end files and if yes then is it enough just to do SGA-ICE.py /path/to/fastq/data/ -t 8?

Thank you in advance.

Best wishes,

Michal

[juroscito]

Hi Michal,
SGA-ICE does not distinguish between paired or single-end reads. Each fastq file will be corrected separately in all k-mer/overlap correction runs.
You should keep in mind that the fastq files will be converted to .fasta after the overlap correction (which, by default, is on) because the quality values in the .fastq file produced by 'sga correct' not always match the length of the reads if insertions/deletions were corrected (newer versions of sga may solve this issue).

And yes, calling SGA-ICE like that will work just fine. You'll get a runMe.sh file with all the commands that will run.

Let us know if you have further questions.

Best,
Juliana

[mictadlo]

Hi Julian,
SGA-ICE has finished to run and created fastq and fasta files. How could I check whether the quality values in the .fastq file produced by 'sga correct' did not match the length of the reads if insertions/deletions were corrected?

Thank you in advance.

Mic

[juroscito]

Hi Michal,

The quality values are not changed throughout the correction process. This means that the length of the quality val. line in all fastq files are always the same (equal to the original sequence length), irrespective of the sequence being corrected for insertions or deletions.

Does that answer your question?
Best,
Juliana

[mictadlo]

Hi Juliana,
I still not quite sure. SGA-ICE has created 4 output files:

• R1-rmdup-rmconta.final.ecOv.fastq
• R1-rmdup-rmconta.final.ecOv.fastq.fasta
• R2-rmdup-rmconta.final.ecOv.fastq
• R2-rmdup-rmconta.final.ecOv.fastq.fasta

I was not sure and I took the FASTA files for SGA assembly. Below are steps which shows how I try to run SGA to get an assembly. Should I have used the FASTQ files rather than FASTA file for SGA assembly and does the below steps make sense?

#! /bin/bash -x

cd inputFq/ec/

IN1=R1-rmdup-rmconta.final.ecOv.fastq.fasta
IN2=R2-rmdup-rmconta.final.ecOv.fastq.fasta

#
# Parameters
#


# The number of threads to use
CPU=10

# The minimum length of contigs to include in a scaffold
MIN_CONTIG_LENGTH=200


#
# Preprocessing
#

# Preprocess the data to remove ambiguous baseballs
sga preprocess --pe-mode 1 -o reads.pp.ec.fasta $IN1 $IN2

#
# Primary (contig) assembly
#

# Index the corrected data.
sga index -a ropebwt -t $CPU reads.pp.ec.fasta

# Remove exact-match duplicates and reads with low-frequency k-mers
sga filter --homopolymer-check --low-complexity-check -t $CPU reads.pp.ec.fasta

# Compute the structure of the string graph
sga overlap -t $CPU reads.pp.ec.filter.pass.fa

# Perform the config assembly
sga assemble -o primary reads.pp.ec.filter.pass.asqg.gz

#
# Scaffolding
#

PRIMARY_CONTIGS=primary-contigs.fa
PRIMARY_GRAPH=primary-graph.asqg.gz

# Align the reads to the contains
bwa index $PRIMARY_CONTIGS
bwa aln -t $CPU $PRIMARY_CONTIGS $IN1 > $IN1.sai
bwa aln -t $CPU $PRIMARY_CONTIGS $IN2 > $IN2.sai
bwa sampe $PRIMARY_CONTIGS $IN1.sai $IN2.sai $IN1 $IN2 | $SAMTOOLS_BIN view -Sb - > libPE.bam

# Convert the BAM file into a set of contig-contig distance estimates
sga-bam2de.pl -n $MIN_PAIRS -m $MIN_CONTIG_LENGTH --prefix libPE libPE.bam

# Compute copy number estimates of the contains
sga-astat.py -m $MIN_CONTIG_LENGTH libPE.bam > libPE.astat

# Build the scaffolds
sga scaffold -m $MIN_CONTIG_LENGTH -a libPE.astat -o scaffolds.scaf --pe libPE.de $PRIMARY_CONTIGS

# Convert the scaffolds to FASTA format
sga scaffold2fasta --use-overlap --write-unplaced -m $MIN_CONTIG_LENGTH -a $PRIMARY_GRAPH -o sga-scaffolds.fa scaffolds.scaf

Thank you in advance.

Mic

[juroscito]

Hi Mic,
Yes, your SGA commands for assembly seem fine. You don't need the " sga preprocess" step, since it has already been done with SGA-ICE. After SGA-ICE overlap correction, you index the *fastq.fasta with sga index and give the fasta files to sga filter and overlap.

As I said before, the quality value string remains untouched throughout the correction process. Therefore, after correcting for indels with "sga correct --algorithm overlap" , you will see mismatches between sequence length and qual.value string length in the *ecOv.fastq files. That's why we included the fastq>fasta conversion.
SGA assembler can handle both fasta and fastq, but I don't think it can deal with this mismatch in seq length and quality value length if you provide the ecOv.fastq for assembly (I'm not sure, though). Bottom line is, you'd better use the ecOv.fastq.fasta for the downstream steps of SGA, as you've done in your commands.

Keep in mind that other assemblers might need fastq as input, so if you ever decide to try something else than SGA you will need to fake quality values for the ecOv.fastq.fasta instead of using the ecOv.fastq.

Hope that helps.
Best,
Juliana

[mictadlo]

Hi Juliana,
How is it possible to merge two R1*.final.ecOv.fastq.fasta and R2*.final.ecOv.fastq.fastawithout sga preprocess?

Or do I have to run twice index, filter, overlap and assemble for R1 and R2 reads separately?

Thank you in advance,

Best wishes,

Michal

[juroscito]

Hi Mic,
you can cat the files into a single .fasta, as SGA-ICE does after every correction step, then index it. In the runMe.sh, you'll see that after each kmer correction we print the following commands:

cat _.ec.k$kmer.f_q > all.ec.k$kmer.fastq
sga index -a ropebwt -t 4 all.ec.k$kmer.fastq

You could just do the same, but with the *fastq.fasta

Cheers

[mictadlo]

HI Juliana,
Thank you. Are the following commands correct:

# merging
cat R*-rmdup-rmconta.final.ecOv.fastq.fasta > reads.pp.ec.fasta

# Index the corrected data.
sga index -a ropebwt -t $CPU reads.pp.ec.fasta

# Remove exact-match duplicates and reads with low-frequency k-mers
sga filter --homopolymer-check --low-complexity-check -t $CPU reads.pp.ec.fasta

# Compute the structure of the string graph
sga overlap -t $CPU reads.pp.ec.filter.pass.fa

# Perform the config assembly
sga assemble -o primary reads.pp.ec.filter.pass.asqg.gz

Thank you in advance,

Best wishes,

Mic

[juroscito]

Hi Mic,
Yes, the commands are ok.
Hope you get a nice assembly out of it, we would be happy to know the results.

Cheers,
Juliana

[mictadlo]

Thank you, I will do.