hypothesis - high error rate in genes in flye assembly
Hypothesis: The genes in the flye assembly have errors:
- no effect on
prodigalresults because these genes are not necessarily partial - less hits in
diamondDB - low query/subject ratio - unique genes in SR vs. LR. vs hybrid
Designs
Relates to
- #73Analyses
Activity
changed milestone to %Analyses
added Analysis Investigation labels
barrnapIdentifying and comparing rRNA genes
- repo
- supported organisms: Bacteria, Archaea, Eukaryota, Metazoan Mitochondria
Edited by Valentina Galata-
find . -name "GDB_*.fna" | sort | while read f; do for k in mito bac arc euk; do echo "${f}" && barrnap --threads 10 --kingdom ${k} --outseq barrnap/$(basename -s '.fna' ${f}).${k}.fa < ${f} > barrnap/$(basename -s '.fna' ${f}).${k}.gff; done; donefind barrnap/ -type f -name "*.bac.gff" | sort | while read f; do echo -e "${f}\n$(grep -v "^#" ${f} | cut -f9 | cut -d";" -f1,2 | sed 's/;/\t/g' | sed 's/(partial)/\tpartial/' | sort | uniq -c | sort -k2)"; doneEdited by Valentina Galata
-
-
hmmsearch --cut_ga --cpu 10 --tblout bac.tblout -A bac.align Bacteria.hmm prot.faa > /dev/nullEdited by Valentina Galata
marked this issue as related to #69 (closed)
-
DIAMONDparameters to try:
-
--sensitive,--more-sensitive- w/o any sensitivity option: hits of >70% identity and short read alignment
-
--query-coverand--subject-cover, e.g.95or higher -
--id, e.g.95 -
--unal 1to output unaligned queries
Edited by Valentina Galata -
-
in an interactive session:
srun -p bigmem --time=1:00:0 -N1 -n1 -c5 --pty bash -iEdited by Valentina Galata
marked this issue as related to #73 (closed)
SMALLER GENES DETECTED IN FLYE - possible reasons (lest, I forget)
- the median size of small genes detected in flye is
~400 bp, and there are others even smaller - given the SR sequencing is typically 2x150 bp, it is plausible that we miss these “smaller” genes with the SR sequencing and subsequent assembly
- also the possible explanation for the disparity with the “UniProt/TrEmbl” database, i.e. the
0.5-07 ratioswhen mapping to TrEmbl - explanation: the TrEmbl database was in the past created using SR-sequencing, i.e. “short genes” may not have been recorded then. whereas ``flye seems to capture them because of the
possibly longersequences being fed into the sequencer itself.
- the median size of small genes detected in flye is
-
New findings: In the
flyeassembly, many genes with low query/subject length ratio (<= 0.5) appear in clusters. This is not the case of themegahitassembly.- filter
diamondhits be query/subject length ratio (<= 0.5) - sort by query ID (contains contig and gene number)
- find clusters of consecutive gene numbers (and same subject ID)
/scratch/users/vgalata/GDB/results/annotation/diamond/lr/flye/proteins.filtered.annot.tsv: looking for consecutive genes (clusters) where subject ID does not matter 75460 total hits 49632 (65.77%) queries are in clusters 14770 clusters w/ average length of 3.36 genes /scratch/users/vgalata/GDB/results/annotation/diamond/lr/flye/proteins.filtered.annot.tsv: looking for consecutive genes (clusters) where subject ID does matter 75460 total hits 24253 (32.14%) queries are in clusters 9810 clusters w/ average length of 2.47 genes/scratch/users/vgalata/GDB/results/annotation/diamond/sr/megahit/proteins.filtered.annot.tsv: looking for consecutive genes (clusters) where subject ID does not matter 20397 total hits 3016 (14.79%) queries are in clusters 1459 clusters w/ average length of 2.07 genes /scratch/users/vgalata/GDB/results/annotation/diamond/sr/megahit/proteins.filtered.annot.tsv: looking for consecutive genes (clusters) where subject ID does matter 20397 total hits 250 (1.23%) queries are in clusters 124 clusters w/ average length of 2.02 genes - filter
-
Quick check how many proteins in
flyehave low q/s length ratio and how many were not found inmegahit:cd /scratch/users/vgalata/GDB/results # intersections comm -13 <(tail -n +2 annotation/diamond/lr/flye/proteins.filtered.annot.tsv | cut -f 1 | sort) <(grep "^>" analysis/cdhit/sr_megahit__lr_flye__test.faa | sed "s/^>//" | cut -d" " -f1 | sort) | wc -l # 29721 comm -23 <(tail -n +2 annotation/diamond/lr/flye/proteins.filtered.annot.tsv | cut -f 1 | sort) <(grep "^>" analysis/cdhit/sr_megahit__lr_flye__test.faa | sed "s/^>//" | cut -d" " -f1 | sort) | wc -l # 40897 comm -12 <(tail -n +2 annotation/diamond/lr/flye/proteins.filtered.annot.tsv | cut -f 1 | sort) <(grep "^>" analysis/cdhit/sr_megahit__lr_flye__test.faa | sed "s/^>//" | cut -d" " -f1 | sort) | wc -l # 34563 # total counts grep "^>" analysis/cdhit/sr_megahit__lr_flye__test.faa | wc -l # 64284 wc -l annotation/diamond/lr/flye/proteins.filtered.annot.tsv # 75461 annotation/diamond/lr/flye/proteins.filtered.annot.tsv # - 1 for the headerEdited by Valentina Galata
mentioned in commit 14d7afdf
-
| From: | Mikhail Kolmogorov <fenderglass@gmail.com> | |----------|--------------------------------------------------------------------------------------| | To: | Cedric Christian LACZNY <cedric.laczny@uni.lu> | | Cc: | Susheel Bhanu BUSI <susheel.busi@uni.lu>, Valentina GALATA <valentina.galata@uni.lu> | | Subject: | Re: Question about discrepancy in assemblies using ONT and ONT + Illumina | | Date: | Wed, 19 Aug 2020 14:09:53 -0700 (19/08/20 23:09:53) | I think it is very likely that the effect is caused by the remaining indels - something that we have also seen in our analysis. My guess is that an indel can potentially split a gene into parts, thus artificially increasing the protein counts. You might try to check the predicted ORF length distribution might help to find out if this is indeed the case (e.g. I expect Flye to have more short proteins and less long proteins, as compared to the other assemblers in this case). Also, it makes sense to polish Flye assembly with Illumina and check if it changes the statistics. We used both Prodigal and MetaGeneMark (incorporated into QUAST) for the analysis, but both tools seem to be in agreement at large.
Please register or sign in to reply