Banana Slug Genomics

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== Mitochondrion ======

The mitochondrion was assembled by Kevin Karplus in the assemblies/slug/barcode-of-life/ directory. The reason for the strange name for the directory was that at first the attempt was just to recover the COX1 gene that is used for the [[http://www.boldsystems.org/|BOLD (barcode of life database)]] project to characterize eukaryotes by their mitochondrial sequences.  When it became clear that the whole mitochondrial genome was well covered in the Illumina data, the project switched to trying to reconstruct the full mitochondrial genome.

===== Method =====

There were many iterations and many different attempts at assembling the mitochondrial genome.  Most of the iterations consisted of taking some draft genome, mapping all the Illumina reads to it using bwa, and selecting out all reads that mapped and their paired ends (even if the pairs didn't map), and reassembling those reads.

The initial attempts used just some contigs from the previous year's (2010) attempt at a whole-genome assembly using SOAPdenovo.  Later, I found that the 2011 SOAPdenovo assembly had most of the genome in one contig (ending at two repeat regions).  Both SOAPdenovo and abyss were used to try to assemble the reads.  Generally, abyss got somewhat longer longest contigs, and the two assemblers agreed on the parts they had in common.
   * Started with a search of SOAPdenovo-assembly1/k31/soapSlug.scafSeq for scaffolds that matched examples from other mollusks.
   * Looked for 454 reads that extended or joined contigs in scaffold
   * Repeated (sometimes using more sensitive searches) until no more credible scaffolds from the SOAPdenovo-assembly1/k31/ assembly nor 454 reads were found.
   * The 454 coverage of the mitochondrion is so slight as to be nearly useless, so instead we can iterate:
        - find all Illumina reads that map to the mitochondrial draft, using BWA
        - assemble them using SOAPdenovo.
   * It looks like the Illumina reads have about 228x coverage of the mitochondrion, but coverage is patchy, and it seems to be difficult to close the circle (at least with SOAPdenovo).  
      
   * It turns out that a lot of the hard hand work and iterated searching to assemble the mitochondrion was not necessary, as the SOAPdenovo-assembly2/all/k63/illumina-454-all_63-mers.scafSeq assembly now has a 14960-long contig (not scaffold!) which is an almost-full-length mitochondrion, roughly as good as the best I'd managed to assemble from the SOAPdenovo-assembly1/ bits.  
   * Iterating mapping reads with BWA and assembling them with SOAPdenovo made some progress, but there was a gap that just wouldn't close.
   * Switching to abyss (version 1.2.7) for the assembly of the reads made a much larger contig (15535-long after pasting on a suggestion from one abyss assembly onto another).
   * Iterating search and abyss assembly does not lengthen the large contig.  Cleaning up and calling the consensus with bwa+samtools+bcftools doesn't change things much either.  There seems to be a large variation in coverage (from 20x to 2300x, with a median of 225x), so I suspect that there is a repeat region at the beginning of the current contig that may have 10 repeats in it.

Alternating finding new reads and assembling them made very slow progress, because the new reads only extended the assembled region by 50–100 bases.  Eventually, I wrote a new program (look-for-exit [needs its own page FIXME]) to manually extend the contigs and find exits from repeat regions, being more aggressive in extending the contig than the automatic assemblers.  I was eventually able to close the circle this way, and get a complete genome, though there is one repeat region with long repeats (about a dozen copies of a 615±1 long repeat) that I could not order, because the differences between repeats were far enough apart that I couldn't disambiguate the order with the [[bioinformatic_tools:bwa#determining_paired-end_insert_size|short fragment lengths]] of the data available.  I think I have all the variants of repeat, but in some cases I can't even tell which first half of the repeat goes with which second half.

At some point in the process, I rotated the genome to correspond to the closest previous mitochondrial genome: //Biomphalaria glabrata// strain M, a gastropod.

I used bwa with mpileup to make new consensus sequences, and iterated that a few times to get about as good an assembly as I can without some longer fragments or some PCR to determine the order.

===== Reads =====

One of the biggest tasks in assembling the mitochondrion was separating the mitochondrial reads from the full set of reads.
I needed a seed to start with, which I initially found by searching the earliest assembly for a contig that matched the COX1 protein from several different mollusks using tblastn, and to cox1 genes using blastn.  I got 3 contigs for this small gene in the low-quality assembly from SOAPdenovo_assembly1/k31, and had to do iterative searching to extend the set.  (I could have saved about 10 iterations by starting from the SOAPdenovo-assembly2/all/k63/ assembly instead, but it had not been completed when I started the search for mitochondrial reads.

The selection process used bwa to map reads that had been run through seqprep to the mitochondrial contigs found so far, keeping any reads that mapped or whose paired read mapped.  I sometimes used the quake-corrected set used for this year's assembly, and sometimes used the set that had been run through seqprep but not corrrected with quake.

In barcode-of-life/map-Illumina-raw-45/, which has a draft of the complete mitochondrial genome to map to, I mapped the reads that had not been run through quake (that is what the "raw" signifies in the directory name, though the reads had been run through seqprep, and so were not really the raw reads). I found 61,428 reads that are most likely mitochondrial.

>samtools flagstat merged.bam                                          
  61428 + 0 in total (QC-passed reads + QC-failed reads)                 
  0 + 0 duplicates                                                       
  60932 + 0 mapped (99.19%:nan%)                                         
  35560 + 0 paired in sequencing                                         
  17780 + 0 read1                                                        
  17780 + 0 read2                                                        
  31562 + 0 properly paired (88.76%:nan%)                                
  34568 + 0 with itself and mate mapped                                  
  496 + 0 singletons (1.39%:nan%)                                        
  0 + 0 with mate mapped to a different chr                              
  0 + 0 with mate mapped to a different chr (mapQ>=5)    

I'm wondering about the 3004 pairs that map both reads but are not "properly paired".  Are these indicating a misassembly, or is this just mismapping in the repeat region?                                  

I'm also curious about whether I could have found the mitochondrial reads more simply by extracting every read that had a high-count k-mer.  Looking at the jellyfish 19-mer counts from this set, I see counts as high as 5707 (for AATAAAGACTAATAAGATA) and 12,976 19-mers with multiplicity ≥100. But in the jellyfish counts for the whole set of reads, I see 7 low-complexity kmers that occur over 10 million times:               
  CCCTAACCCTAACCCTAAC     10752814                                       
  AACCCTAACCCTAACCCTA     10736174                                       
  CCTAACCCTAACCCTAACC     10715055                                       
  CTAACCCTAACCCTAACCC     10710935                                       
  AGGGTTAGGGTTAGGGTTA     10654682                                       
  ACCCTAACCCTAACCCTAA     10607845                                       
  AAAAAAAAAAAAAAAAAAA     10486112                                       

Even if entire reads are full of CCCTAA, that is still hundreds of thousands of reads, so separating by frequency would not have been good enough.                                                           
                                                                       
The most frequent 19-mer in the subset occurs 6035 times in the full set (so I may be missing 272 copies in the subset), but there are almost 209,000 more common 19-mers, so selecting by frequency would have gotten me mostly low-complexity junk, not mitochondrial sequence.      

After cleaning the mitochondrial reads with [[bioinformatic_tools:jellyfish|jellyfish]] and [[bioinformatic_tools:quake|quake]], 
  clean_19_dir/merged_1.fastq has 1253271 bases in 16885 reads.
  clean_19_dir/merged_2.fastq has 1252843 bases in 16885 reads.
  clean_19_dir/merged.fastq has 2860095 bases in 26498 reads. 
  Total: 5,366,209 bases in 60,268 reads. 

The estimated coverage of the genome from fitting a gamma distribution to the 19-mers is about 216x, giving a genome-size estimate of 24,844 bases, not far from the 23642 I currently have---perhaps there are a couple more repeats than I've been estimating?  But the accuracy of the genome size estimates is not likely to be that precise. 

===== Mitochondrial sequence =====

The first draft sequence is available as {{mitochondrion-draft1.fasta.gz|gzipped fasta file}}.  This corresponds to /campusdata/BME235/assemblies/slug/barcode-of-life/map-Illumina-raw-45/consensus-6 on campusrocks.

===== Annotation =====

I sent the sequence to [[http://dogma.ccbb.utexas.edu/]] for annotation, and it quickly created a very crude web interface that almost works.  Unfortunately, it seems to have missed some of the protein genes and it provides no way to download its annotation in a GENBANK format.  Since there are hundreds of tRNA genes (the repeat region is full of them), I'm not tempted to try screen-scraping.  A better way to annotate the mitochondrion needs to be found.
                                                                                                         

The repeat region starts at position 7037 in draft1, with CTGTAAGAGAATTATTTTAGTAATAAAATTTAATTTTAAGAAAAGAATTTTTCT
There are 12 full (615±1 long) repeats there, then a partial one (456 long) that ends with 
CTGTAAGAGAATTATTTTAGTAATAAAATTTAATTTTAAGAAAAGAATTTTTCT at 14867.

The closest sequenced mitochondrial genomes at NCBI are

  NC_010220.1     Biomphalaria tenagophila mitochondrion, complete genome
  >gb|EF433576.1| Biomphalaria tenagophila strain Taim-RS mitochondrion, complete genome
       	max     total   query   E-value max
	score   score   coverage        identity
	2679    4235    51%     0.0     71%

  NC_005439.1             Biomphalaria glabrata mitochondrion, complete genome
  >gb|AY380531.1| Biomphalaria glabrata strain 1742 mitochondrion, complete genome
 	max     total   query   E-value max
	score   score   coverage        identity
	2562    3934    52%     0.0     69%

So this is pretty far from the previously known genomes.
Banana Slug Genomics

User Tools

Site Tools

Page Tools
