Answers.Body: <p>A relatively convenient method would be to create a .chain file, as used by tools like <code>liftOver</code> and <code>crossmap</code> (and probably remap from NCBI). For your use-case, you can simply use the <a href="http://genomewiki.ucsc.edu/index.php/DoSameSpeciesLiftOver.pl" rel="nofollow noreferrer">DoSameSpeciesLiftOver.pl</a> 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.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-11-08T14:28:39.067Z
Answers.Text: This post helped with the ssh error: https://stackoverflow.com/questions/17335728/connect-to-host-localhost-port-22-connection-refused#22255174
Answers.UserDisplayName: Biomagician
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-11-08T17:10:11.910Z
Answers.Text: Ah, it uses ssh to submit commands to itself as an ersatz cluster.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-08-28T19:48:11.213Z
Answers.Text: 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
Answers.UserDisplayName: Biomagician
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-08-28T20:26:03.727Z
Answers.Text: 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.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-08-31T18:17:10.087Z
Answers.Text: 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.
Answers.UserDisplayName: Biomagician
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-09-01T16:38:59.103Z
Answers.Text: I don’t know where ssh is coming from there either.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-08-14T08:29:56.280Z
Answers.LastActivityDate: 2018-08-14T08:29:56.280Z
Answers.UserDisplayName: Devon Ryan
Body: <p>I am working with two versions of the <em>C. elegans</em> 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 <code>show-coords</code> on the delta file and obtained a .coords file. It seems more human-friendly. The first few lines look like this:</p>
  
  <pre><code>/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
  </code></pre>
  
  <p>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:</p>
  
  <pre><code>   1   237139  |     4022   241181  |   237139   237160  |    99.99  | 15072434 15371751  |     1.57     1.54  | I  I
  </code></pre>
  
  <p>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?</p>
  
  <p>Can anybody guide me, please?</p>
  
  <p>I've also aligned the two genomes using Mauve:</p>
  
  <pre><code>progressiveMauve c_elegans.PRJNA13758.WS245.genomic.fa scaffolds.fa --output scaffolds.mauve &gt; scaffolds.txt
  </code></pre>
  
  <p>I've tried to import the .mauve file into Python using Biopython but I got an error:</p>
  
  <pre><code>alignment = SeqIO.read('scaffolds.mauve', "mauve")
  
  Traceback (most recent call last):
    File "code/mauve.py", line 4, in &lt;module&gt;
      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
  </code></pre>
  
  <p>Note that the alignment opened in the Mauve GUI looks fine. Unfortunately, the picture won't load here.</p>
  
  <p>I am working on Mac OSX.</p>
  
  <p>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.</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-08-28T14:01:09.977Z
Comments.Text: @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.
Comments.UserDisplayName: Biomagician
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-08-28T15:33:38.673Z
Comments.Text: Hmmmm, okay; that was not clear to me from your original post.  I see that one of the files you mention has ‘WS245’ in the name and I assumed that was a specific official build of the C. elegans genome, and I erroneously concluded that the other ‘version’ (you mention version 1 and version 2) was a different specific official genome build, and my comments were for that particular situation.  If you are using your own, homemade versions wont you eventually have to map the genes of interest back to one of the official genome builds (i.e., prior to publication)?
Comments.UserDisplayName: mdperry
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-08-28T19:51:47.033Z
Comments.Text: That's a good question. For the time being, I use my own build because I think it's better than WormBase's. But soon, they'll publish a better version and I'll use that one for the publication. However, I'll still have a situation of one official version + personal version. As far as I understand, REMAP only works with official versions.
Comments.UserDisplayName: Biomagician
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-08-09T08:08:01.333Z
Comments.Text: `SeqIO.read` is for a single record, try `SeqIO.parse`
Comments.UserDisplayName: Chris_Rands
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-08-10T13:12:03.610Z
Comments.Text: 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.
Comments.UserDisplayName: mdperry
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-08-11T09:14:47.320Z
Comments.Text: I don't see where I can feed my own file to the Remap service of NCBI: https://www.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."
Comments.UserDisplayName: Biomagician
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-08-16T18:27:56.127Z
Comments.Text: 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.
Comments.UserDisplayName: mdperry
Tags: <phylogenetics><python><mauve>
Title: Find a map/correspondence between two versions of a genome
ThreadText: # Find a map/correspondence between two versions of a genome
  
  (Asked by: username_0 on 2018-08-08)
  
  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.
  
  **Tags:** phylogenetics, python, mauve
  
  > **Comments:**
  > **username_2**: `SeqIO.read` is for a single record, try `SeqIO.parse`  
  > **username_1**: 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.  
  > **username_0**: I don't see where I can feed my own file to the Remap service of NCBI: https://www.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."  
  > **username_1**: 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.  
  > **username_0**: @username_1 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.  
  > **username_1**: Hmmmm, okay; that was not clear to me from your original post.  I see that one of the files you mention has ‘WS245’ in the name and I assumed that was a specific official build of the C. elegans genome, and I erroneously concluded that the other ‘version’ (you mention version 1 and version 2) was a different specific official genome build, and my comments were for that particular situation.  If you are using your own, homemade versions wont you eventually have to map the genes of interest back to one of the official genome builds (i.e., prior to publication)?  
  > **username_0**: That's a good question. For the time being, I use my own build because I think it's better than WormBase's. But soon, they'll publish a better version and I'll use that one for the publication. However, I'll still have a situation of one official version + personal version. As far as I understand, REMAP only works with official versions.  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_3 on 2018-08-14. Score: 1 | Accepted answer
  
  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](http://genomewiki.ucsc.edu/index.php/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.
  
  > **Comments:**
  > **username_0**: 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  
  > **username_3**: 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.  
  > **username_0**: 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.  
  > **username_3**: I don’t know where ssh is coming from there either.  
  > **username_0**: This post helped with the ssh error: https://stackoverflow.com/questions/17335728/connect-to-host-localhost-port-22-connection-refused#22255174  
  > **username_3**: Ah, it uses ssh to submit commands to itself as an ersatz cluster.  
  
  ---
  
  