This question is slightly related to this one: Improve a reference genome with sequencing data
Only in my case, I have a starting genome reference A, of a short phage genome, and reads from a related genome B, for which I haven't been able to find a reference.
I took the Illumina 2x75bp paired end reads from genome B and ran Velvet on it (I ran Ray, but it didn't complete). It gives me a fragmented collection of contigs. The contigs align 99-98% to reference genome A. I don't have long reads.
My question is: with the short reads and the contigs that I have, can I transform the reference genome A into a reference genome B?
PS: no, I don't have long reads, and there is no budget for them. The budget for long reads is zero and will remain zero. Even if we could sequence long reads for free, we would still not do it because of the lab work.
Tools that I have tried so far:
ococo:
$HOME/ococo/ococo -i contigs.bwamem.bam -c 1 -f $HOME/Genomes/T2_phage/starting.fasta -F consensus.fasta --vcf-cons - -V variants.vcf
It produced a consensus.fasta, and checking how much it changes from the reference, it changes. There are one or more changes in about 33% of the 100bp windows:
diff -y <(seqtk seq -l 100 $HOME/Genomes/T2_phage/starting.fasta) <(seqtk seq -l 100 consensus.fasta) | grep '\ \|\ ' | wc -l
1272
diff -y <(seqtk seq -l 100 $HOME/Genomes/T2_phage/starting.fasta) <(seqtk seq -l 100 consensus.fasta) | grep -v '\ \|\ ' | wc -l
419