2016-09-19 15:04 GMT+02:00 Bastien Chevreux <bach@xxxxxxxxxxxx>:
On 19 Sep 2016, at 8:16 , Sven Klages <sir.svencelot@xxxxxxxxx> wrote:
Currently no idea of how diverse this repository is.
We have ~600bp PCR fragments containing 3 variable regions of roughly
30,30 and 60 bases.
Bonus question: how far away are these regions from each other? Then
again, it does not really matter as there is a good chance for them to be
“too far” for 100bp reads.
It’s probably too late, but I would have envisaged 250bp MiSeqs with a
library size of 400bp for sequencing that. With amplicons of size 600 I
might be even tempted to try not to fragment that.
Totally agreed. But I have PE125 data and need to proceed, .. somehow :-)
For each pool we have 7.5mio reads (3.8mio pairs) and 6.8mio reads (3.4mio
pairs).
I suspect you will have another problem waiting for you: sequencing
errors. Let’s assume a 1% error probability for Illuminas. In practice,
some regions will have more errors and others less than that, but it’s a
good approximation to show things.
So, for every 100x coverage, each position will have 1 erroneous base in
on one the reads covering it. With 300x, you have on average one of each of
the other bases as error. For bases to be valid, MIRA wants minimum 3x
forward and 3x reverse occurrence, so 300x*6 = 1800x. With that coverage,
you have on average every consensus base covered by “valid” alternative
bases … and MIRA will see that and create different contigs based on that.
As I said, in practice some regions need way less coverage for this to
happen. And I am not sure atm I made *all* the algorithms available to
parametrisation to change this 3/3 rule of MIRA.
What I would try if I had to do this project:
1. preprocess only of data, maybe each lib separately:
job=genome,accurate,de-novo
parameters= -GE:ppo=yes
2. Fetch preprocessed data, back to (clipped) FASTQ:
miraconvert -F -C …_assembly/…_d_chkpt/readpool.maf cleaned_data.fastq
3. Merge pairs on cleaned data. Conservative merging: 30bp overlap with no
errors. I use FLASH, like this: ...
4. Use only merged pairs (“.extended” file from FLASH), filtered for a
length which is, say, 10 to 20 bp larger than the longest non-variable
stretch. let’s say that this is 300bp … and this is the point where this
strategy fails with 100bp PE on your data :-(
miraconvert -X 300 …
That last step would result in 0 reads with 100bp PE.
5. Assemble only the surviving merged reads, switching of Illumina adaptor
clipping.
Every other strategy I can think of will create fragments and not fullThanks for your comment, .. I will think about a strategy ...
amplicons and that is somehow depressing.
