0
$\begingroup$

I am working with two versions of the C. elegans genome. I am finding interesting regions (specifically, tRNA genes) in version 1 and then I would like to know if version 2 also has a tRNA gene in that region. For this, I need to know the correspondence between coordinates in genome 1 and coordinates in genome 2. I think a whole-genome alignment programme might help me, e. g. nucmer or mauve. I have run nucmer on the two genomes and got a delta file but I don't know how to interpret it. I have run show-coords on the delta file and obtained a .coords file. It seems more human-friendly. The first few lines look like this:

/Users/user/Documents/project/output/genome/celegans/ref/wormbase/seq/c_elegans.PRJNA13758.WS245.genomic.fa /Users/user/Documents/project/output/genome/celegans/strain/seq/bristol/scaffolds.fa
NUCMER

    [S1]     [E1]  |     [S2]     [E2]  |  [LEN 1]  [LEN 2]  |  [% IDY]  |  [LEN R]  [LEN Q]  |  [COV R]  [COV Q]  | [TAGS]
===============================================================================================================================
       1      432  |      498      929  |      432      432  |   100.00  | 15072434 15371751  |     0.00     0.00  | I  I
       1      432  |     1002     1433  |      432      432  |   100.00  | 15072434 15371751  |     0.00     0.00  | I  I

Does this mean that chromosome I starting at 1 and ending at 432 of version 1 corresponds to chromosome I starting at 498 and ending at 929 and that there is 100% identity between these two 432-long sequences? In this case, it's easy to find a correspondence between bases in version 1 and bases in version 2. However, some lines do not align sequences of the same length:

   1   237139  |     4022   241181  |   237139   237160  |    99.99  | 15072434 15371751  |     1.57     1.54  | I  I

Also, in this sequence position 1 of chromosome I is aligned again? So is it aligning to more than one position in version 2 of the genome? Also one sequence is 237139 bases long and the other 237160 bases long. How do I know the internal structure?

Can anybody guide me, please?

I've also aligned the two genomes using Mauve:

progressiveMauve c_elegans.PRJNA13758.WS245.genomic.fa scaffolds.fa --output scaffolds.mauve > scaffolds.txt

I've tried to import the .mauve file into Python using Biopython but I got an error:

alignment = SeqIO.read('scaffolds.mauve', "mauve")

Traceback (most recent call last):
  File "code/mauve.py", line 4, in <module>
    alignment = SeqIO.read('scaffolds.mauve', "mauve")
  File "/Users/user/anaconda/lib/python3.6/site-packages/Bio/SeqIO/__init__.py", line 712, in read
    raise ValueError("More than one record found in handle")
ValueError: More than one record found in handle

Note that the alignment opened in the Mauve GUI looks fine. Unfortunately, the picture won't load here.

I am working on Mac OSX.

I would appreciate your help in understanding the delta file, resolving the Python error or any other support in finding a function between coordinates in version 1 of the genome and coordinates in version 2 of the genome.

$\endgroup$
  • 1
    $\begingroup$ SeqIO.read is for a single record, try SeqIO.parse $\endgroup$ – Chris_Rands Aug 9 '18 at 8:08
  • $\begingroup$ NCBI has REMAP page/tool where you can feed it a file with coordinates of sequence features from one build version, and ask for those features mapped to a different build. I have only used it to go forward, not backword. Also, each version of the C. elegans genome used to have a liftover directory at the Sanger ftp site that showed the version changes, as I recall. Check out the WormBase site. $\endgroup$ – mdperry Aug 10 '18 at 13:12
  • $\begingroup$ I don't see where I can feed my own file to the Remap service of NCBI: ncbi.nlm.nih.gov/genome/tools/remap. Reply from NCBI to my request: "The remap service relies on assembly to assembly alignment data to do this remaping/conversion. These datasets came from NCBI's eukaryotic genome assembly workflow. For outside genomic assembly that are not deposited to NCBI, we have no access, so we will not be able to provide mapping between that and those available from the NCBI collection." $\endgroup$ – charlesdarwin Aug 11 '18 at 9:14
  • $\begingroup$ The file that you need is a bed file, or a GFF3 file, or similar, that contains the start and end coordinates for all of your tRNA gene regions from genome Version 1. You select and upload that file and tell REMAP the version/build using one drop down menu. Then in another drop down menu you select the genome version/build for Version 2. Click on the button and in a little while you will have some files to download showing the coordinates for those gene regions in Version 1, and how they map to Version 2. $\endgroup$ – mdperry Aug 16 '18 at 18:27
  • $\begingroup$ @mdperry the versions I am using are personal. I have reassembled the genomes using my own data so I cannot choose a version available on NCBI because I have not deposited my versions yet. $\endgroup$ – charlesdarwin Aug 28 '18 at 14:01
1
$\begingroup$

A relatively convenient method would be to create a .chain file, as used by tools like liftOver and crossmap (and probably remap from NCBI). For your use-case, you can simply use the DoSameSpeciesLiftOver.pl command as described in the link. The resulting chain files are pretty self-explanatory and can be used to map BED files of regions of interest.

$\endgroup$
  • $\begingroup$ I am following all the steps. I am stuck here: faToTwoBit ../../genomes/genbank/GCA_000004515.3_Glycine_max_v2.0_genomic.fna.gz GCA_000004515.3_Glycine_max_v2.0.2bit dyld: Library not loaded: /usr/local/opt/mysql/lib/libmysqlclient.20.dylib Referenced from: /usr/local/bin/faToTwoBit Reason: image not found $\endgroup$ – charlesdarwin Aug 28 '18 at 19:48
  • 1
    $\begingroup$ I suggest you install that with conda instead, then you'll get the mysqlclient dependency (presumably). BTW, I expect you need to gunzip the fasta file. $\endgroup$ – Devon Ryan Aug 28 '18 at 20:26
  • $\begingroup$ I have the following error: "ssh: connect to host localhost port 22: Connection refused" for the command "time (doSameSpeciesLiftOver.pl -verbose=2 -buildDir=pwd \ -ooc=pwd/${target}.ooc -fileServer=localhost -localTmp="/dev/shm" \ -bigClusterHub=localhost -dbHost=localhost -workhorse=localhost \ -target2Bit=pwd/${target}.2bit -targetSizes=pwd/${target}.chrom.sizes \ -query2Bit=../${query}/${query}.2bit \ -querySizes=../${query}/${query}.chrom.sizes ${target} ${query}) > do.log 2>&1". I don't understand why it's trying to ssh. $\endgroup$ – charlesdarwin Aug 31 '18 at 18:17
  • 1
    $\begingroup$ This post helped with the ssh error: stackoverflow.com/questions/17335728/… $\endgroup$ – charlesdarwin Nov 8 '18 at 14:28
  • 1
    $\begingroup$ Ah, it uses ssh to submit commands to itself as an ersatz cluster. $\endgroup$ – Devon Ryan Nov 8 '18 at 17:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.