Lost reads in pair 1 after the update to 1.502

[sam....@gmail.com]

HI Falk,

I updated from version 1.462 to version 1.5 and then 1.502 to get the fix for the combineSamples bug but now when I reran some data I find I have lost nearly all of the reads in pair 1. For example where 70-80% of the pair 1 reads were accepted now only 1-3% of the pair 1 reads are accepted however the number of pair 2 reads is exactly the same. I am using the same options apart from changing the dereplication value.

On top of this the combineSample bug is still not fixed, I have the same error as before.

Please help

Thank you very much, this is a wonderful pipeline.

Sam

[Falk Hildebrand]

Hey Sam,

can you please post the text document "demulti.log" from the log files, so I can see how reads were filtered? This is also a general first step when setting up the parameters in the sdm_opt*.txt files. So what I could think of, is that min read length is set too high, as I changed the behaviour of what part of the 5' DNA molecule will be cut, e.g. Primer are cut by default and this could effect your result.

Also the number of pair1 reads should never be higher than the number of pair2 reads, so this could be an error in the counting (the counting part was also changed a bit, but so far I had no problems). Can you please also post the LotuSrunlog and your sdm options file, used in your run? 

Which combine Samples bug do you mean? 

best,

Falk

[sam....@gmail.com]

HI Falk,

This is going to be a long post.

First the sdm_opt.txt file

#sdm options file to control sequence quality filtering, demultiplexing and preparation (can also be used without demultiplexing)

#* indicates alternative quality filtering options, saved in *.add.fna etc. files separately from initial quality filtered dataset

#sequence length refers to sequence length AFTER removal of Primers, Barcodes and trimming. this ensures that downstream analyis tools will have appropiate sequence information

#options with a star in front are lenient parameters for mid qual sequences (only used for estimating OTU abundance, not for OTU building itself).

minSeqLength 230

maxSeqLength 1000

minAvgQuality 27

*minSeqLength 230

*minAvgQuality 20

#truncate total Sequence length to X (length after Barcode, Adapter and Primer removals, set to -1 to deactivate)

TruncateSequenceLength 230

#Ambiguous bases in Sequence

maxAmbiguousNT 0

*maxAmbiguousNT 1

#sequence is discarded if a homonucleotide run in sequence is longer

maxHomonucleotide 8

#Filter whole sequence if one window of quality scores is below average

QualWindowWidth 20

QualWindowThreshhold 25

#Trim the end of a sequence if a window falls below quality threshhold. Useful for removing low qulaity trailing ends of sequence

TrimWindowWidth 20

TrimWindowThreshhold 25

#Probabilistic max number of accumulated sequencing errors. After this length, the rest of the sequence will be deleted. Complimentary to TrimWindowThreshhold. (-1) deactivates this option.

maxAccumulatedError 0.75

*maxAccumulatedError -1

#Binomial error model of expected errors per sequence (see https://github.com/fpusan/moira), to deactivate, set BinErrorModelAlpha to -1

BinErrorModelMaxExpError 2.5

BinErrorModelAlpha 0.005

#Max Barcode Errors 

maxBarcodeErrs 0

maxPrimerErrs 0

#keep Barcode / Primer Sequence in the output fasta file - in a normal 16S analysis this should be deactivated (0) for Barcode and de-activated (0) for primer

keepBarcodeSeq 0

keepPrimerSeq 0

#set fastqVersion to 1 if you use Sanger, Illumina 1.8+ or NCBI SRA files. Set fastqVersion to 2, if you use Illumina 1.3+ - 1.7+ or Solexa fastq files.  "auto" will look for typical characteristics of either of these and choose the quality offset score automatically.

fastqVersion auto

#if one or more files have a technical adapter still included (e.g. TCAG 454) this can be removed by setting this option

TechnicalAdapter

#delete X NTs (e.g. if the first 5 bases are known to have strange biases)

TrimStartNTs 0

#correct PE header format (1/2) this is to accomodate the illumina miSeq paired end annotations 2="@XXX 1:0:4" insteand of 1="@XXX/1". Note that the format will be automatically detected

PEheaderPairFmt 1

#sets if sequences without match to reverse primer (ReversePrimer) will be accepted (T=reject ; F=accept all); default=F

RejectSeqWithoutRevPrim F

#*RejectSeqWithoutRevPrim F

#sets if sequences without a forward (LinkerPrimerSequence) primer will be accepted (T=reject ; F=accept all); default=F

RejectSeqWithoutFwdPrim F

#*RejectSeqWithoutFwdPrim F

#this option should be "T" if your amplicons are possibly shorter than a single read in a paired end sequencing run (e.g. if the 16S amplicon length is 200bp in a 250x2 miSeq run, set this to "T"). This option increases runtime by 10%, if in doubt just set to "T". *Requires* LinkerPrimerSequence and ReversePrimer to be defined in mapping file.

AmpliconShortPE T

#options for difficulties during sequencing library construction

#checks if pair1 and pair2 were switched (ignore if single read data)

CheckForMixedPairs F

#checks if whole amplicon was reverse-transcribed sequenced (not switched, just reverse translated)

CheckForReversedSeqs F

second the demulti.log from version 1.462

sdm 1.26 beta

Input File:  several

Output File: /scratch/samantha/1154051.qserver.ctrl.ghpc.dk/demulti.1.fna

Statistics of high quality reads

Reads processed: 3,644,658; 3,644,658 (pair 1;pair 2)

Rejected: 987,439; 1,858,140

Accepted: 2,657,219; 1,786,518 (0; 0 were 5'-trimmed)

Singletons among these: 984,268; 113,567

Bad Reads recovered with dereplication: 406,211

Short amplicon mode.

Min/Avg/Max stats Pair 1

     - Seq Length : 230/230/230

     - Quality :   33/36.8/39

     - Median Seq Length : 0, Quality : 0

     - Accum. Error 0.113

Filtered due to:

  < min Seq length (230)  :                0; 98

       -after Quality trimming :           723,014; 1,806,335

  < avg Quality (27)  :                    4,153; 51,613

  < window (20 nt) avg. Quality (25)  :    4,153; 51,613

  > max Seq length (1000)  :               0; 0

  > (8) homo-nt run  :                     219; 98

  > (0) amb. Bases  :                      0; 0

  > (0.75) acc. errors :                   0; 0

  > (2.5) binomial est. errors :           264,081; 0

  -Barcode unidentified (max 0 errors) :   0

SampleID Barcode Instances

arch1 31,744

arch10 38,944

arch11 44,505

arch12 58,726

arch13 93,159

arch14 62,532

arch15 41,793

arch16 25,003

arch17 32,085

arch18 38,487

arch19 52,248

arch2 39,255

arch20 50,865

arch21 116,147

arch22 136,796

arch23 131,341

arch24 42,500

arch25 70,911

arch26 141,422

arch27 77,515

arch28 82,584

arch29 119,740

arch3 45,378

arch30 90,158

arch31 50,965

arch32 76,220

arch33 105,960

arch34 82,798

arch35 82,637

arch36 127,717

arch37 154,029

arch38 90,545

arch39 79,578

arch4 74,113

arch40 47,688

arch41 82,129

arch42 57,786

arch43 61,009

arch44 57,900

arch5 78,000

arch6 25,615

arch7 16,887

arch8 16,836

arch9 31,180

third the demulti.log from version 1.502

sdm 1.27 beta

Input File:  several

Output File: /scratch/samantha/1155475.qserver.ctrl.ghpc.dk/demulti.1.fna

Statistics of high quality reads

Reads processed: 3,644,658; 3,644,658 (pair 1;pair 2)

Rejected: 3,575,164; 1,858,140

Accepted: 69,494; 1,786,518 (0; 0 were end-trimmed)

Singletons among these: 24,952; 1,741,976

Bad Reads recovered with dereplication: 8,663

Short amplicon mode.

Min/Avg/Max stats Pair 1

     - Seq Length : 230/230/230

     - Quality :   34/36.6/38

     - Median Seq Length : 0, Quality : 0

     - Accum. Error 0.128

Trimmed due to:

  > 25 avg qual in 20 bp windows :         0; 0

  > (0.75) acc. errors, trimmed seqs :     0; 0

Rejected due to:

  < min Seq length (230)  :                3,521,373; 98

       -after Quality trimming :           43,785; 1,806,335

  < avg Quality (27)  :                    149; 51,613

  < window (20 nt) avg. Quality (25)  :    149; 51,613

  > max Seq length (1000)  :               0; 0

  > (8) homo-nt run  :                     14; 98

  > (0) amb. Bases  :                      0; 0

  > (2.5) binomial est. errors :           9,988; 0

Specific sequence searches:

  -With fwd Primer remaining (<= 0 mismatches) :           116,922; 3,643,774

  -Barcode unidentified (max 0 errors) :                   0

SampleID Barcode Instances

arch1 847

arch10 790

arch11 878

arch12 1,144

arch13 2,163

arch14 1,526

arch15 1,027

arch16 617

arch17 769

arch18 942

arch19 1,331

arch2 1,073

arch20 1,085

arch21 2,228

arch22 2,407

arch23 3,029

arch24 1,049

arch25 1,545

arch26 3,299

arch27 1,693

arch28 1,463

arch29 4,448

arch3 1,242

arch30 3,228

arch31 1,723

arch32 2,824

arch33 3,756

arch34 2,948

arch35 3,129

arch36 4,394

arch37 3,060

arch38 2,185

arch39 1,745

arch4 1,961

arch40 1,076

arch41 1,889

arch42 1,451

arch43 1,411

arch44 1,195

arch5 1,524

arch6 610

arch7 368

arch8 398

arch9 687

[sam....@gmail.com]

The combine samples bug, when I use the mapping file to combine samples it does not write a biom file. Instead it writes a biom.not file and says "Samples are being combined; biom format is not supported by LotuS in this case" in the LotuS_runlog.txt

Thanks 

Sam

[Falk Hildebrand]

Hey Sam,

this is clearly a problem of the read length, as suspected due to the primer removal, your seqs are too short. Change these parameters:

minSeqLength 230

*minSeqLength 230

TruncateSequenceLength 230

to something like

minSeqLength 200

*minSeqLength 180

TruncateSequenceLength 200

this is just guesswork, most likely 210 or 220 would already help a lot, but as you see in the sdm output, most seqs were filtered due to being too short:

Rejected due to:

  < min Seq length (230)  :                3,521,373; 98

       -after Quality trimming :           43,785; 1,806,335

hth,

Falk

[Falk Hildebrand]

Hey Sam, 

this combine samples bug was already fixed several versions ago. Can you make sure that you're running lotus 1.502 (1.50 was under some conditions printing a .biom with all OTUs being 0 abundant)?

It should say something like "combined samples, trying to merge" after the tax assignments.

best,
Falk

[sam....@gmail.com]

HI Falk,

Yes I can confirm I was running 1.502 and it still asys Samples are being combined; biom format is not supported by LotuS in this case

[sam....@gmail.com]

So what I understand, is the change is whether the primers are included or not, so if I subtract the length of my primers from the minSeqLength value I was using then it should be the same as previous lotuS ns. I ask this because I did quite a few runs on a subset of the data to determine which minSeqLength value was best for me.

thanks

Sam

[Falk Hildebrand]

Hey Sam,

I found the bug how this could happen, Can you please try with the attached lotus.pl (just replace the lotus.pl in your dir).

hope this fixes the bug for you,

Falk

[Falk Hildebrand]

Yes, substract the primer length from the previous used length threshholds. But beware that the results will still be different, as the qual in the first few bases (the primer usually) is better than in the end of the read. Alternatively you could leave the Primer on your OTUs (these are after all conserved sequences that are part of the rRNA), by renaming the FwdPrimer / RevPrimer column to something else in the mapping file.

hth,

Falk