<jupyter_start><jupyter_text>Prerequisites[MUMmer](http://mummer.sourceforge.net/)  MUMmer does not run on Windows systems unfortunately. As an alternative, you can use Mauve to align genomes, this has a slightly different output format than MUMmer. [Circos](http://circos.ca/) Programs will be run from the commandline, so you need to be able to open a terminal. Google `terminal` + the name of your operating system, and you'll fiind out how to do it. On a Mac it is in /Applications/Utilities/. In this document, command are written in unix form (Linux/MacOS). Most will work in widnows terminal too, just flip the slashes (change / into \). If not, check out a cheatsheet (e.g. [this one](http://ftp.kh.edu.tw/Linux/Redhat/en_6.2/doc/gsg/ch-doslinux.htm) for 'translations' )  Compare genomes of closely related species  Two relatively closely related pathogenic fungi, _Fusarium oxysporum_ (_Fo_) and _Fusarium graminearum_ (_Fg_) have a very different chromosome number. _Fg_ has 4 chromosomes, while _Fo_ has 15. They also differ in size: _Fg_ has a genome of 36.36 Mb (megabases: 1 million A,T,C or G) and _Fo_ has 57.91 Mb. We want to compare these two genomes and see whether for example one chromsoome in _Fg_ corresponds to two or more chromosomes in _Fo_.For this we need to download and align the genomes of _Fo_ and _Fg_. We will use Circos to visualize these alignments.    0. Project managementLet's first do a little project management, and make a directory structure in which we will save the data we download and the results of our analyses.For this, open a Terminal and create a working directory, something like:  `mkdir CircosPlot_FgFo`  `mkdir` is the command for creating a directory.    This notebook is run from this directory, that is I used `cd` to `c`hange `d`irectory to my working directory `CircosPlot_FgFo`, and used `jupyter notebook` to start this notebook.  Create a folder for data:  `mkdir data/`  We may add other data types to the circos plot (Circos is very suitable for that) so let's make a separate subdirectory for genomesequences:  `mkdir data/genomes/`     We also need a folder to store the results of our comparative analyses:  `mkdir results/`We need a place to store the results from genome comparisons:  `mkdir results/genome_alignments/`   We will visualize these results in a single plot, using Circos:  `mkdir circos/`Circos needs a number of input files, which we will store in a separate directory:  `mkdir circos/input/`Circos needs config files, files in which you state the configuration of your plot, i.e. how you plot will look.  `mkdir circos/input/config_files/`  Circos needs files describing the chromsomes you want to plot:  `mkdir circos/input/karyotypes/`  and the links:  `mkdir circos/input/links/`     Circos output we will store in:`mkdir circos/output/`      It is also good practice to redirect any output (like error messages, progress reports, etc. to a logfile. Let's save these in a separate directory as well:`mkdir docs/`We work from the commandline mostly, and it is vital you save the commands you've typed for future reference. This is like your labjournal. Take it seriously! It is also good practice to write some information in the project and the analyses you've done and where you stored what in a README file. For example a description of the directory structure, etc. as I have given above.  One way to combine all this is by using a notebook, like this one.       1. Download genomes of _Fusarium oxysporum_ and _Fusarium graminearum_Go to the [genbank website](https://www.ncbi.nlm.nih.gov/).  From `All Databases` (on the left), select `Assembly` and type `Fusarium oxysporum f. sp. lycopersici 4287` in the search bar [link-to-result](https://www.ncbi.nlm.nih.gov/assembly/?term=Fusarium+oxysporum+f.+sp.+lycopersici+4287).  You can choose from three different assemblies. Pick the one with `Assembly level: Chromosome` [ASM14995v2](https://www.ncbi.nlm.nih.gov/assembly/GCF_000149955.1).Click on the [link to the assembly](https://www.ncbi.nlm.nih.gov/nuccore/AAXH00000000.1/), scroll down, click on the chromosome sequences (accession starts with CM): [link](https://www.ncbi.nlm.nih.gov/nuccore?term=CM000589:CM000603[PACC]). From the menu on the right `Send to:` select `file`, format `FASTA`, sort by `Accession`. You will download a file called `sequence.fasta.txt`, rename this file (`mv` (`m`o`v`e) in the terminal) to `Fol4287.fasta`. Store it in `data/genomes/` (e.g. `mv Downloads/sequence.fasta.txt CircosPlot_FgFo/data/genomes/Fol4287.fasta`)  Check the contents of this file using e.g `grep '>' ~home/CircosPlot_FgFo/data/genomes/Fol4287.fasta`.  You should get:\>CM000589.1 Fusarium oxysporum f. sp. lycopersici 4287 chromosome 1, whole genome shotgun sequence  \>CM000590.1 Fusarium oxysporum f. sp. lycopersici 4287 chromosome 2, whole genome shotgun sequence  \>CM000591.1 Fusarium oxysporum f. sp. lycopersici 4287 chromosome 3, whole genome shotgun sequence  \>CM000592.1 Fusarium oxysporum f. sp. lycopersici 4287 chromosome 4, whole genome shotgun sequence  \>CM000593.1 Fusarium oxysporum f. sp. lycopersici 4287 chromosome 5, whole genome shotgun sequence  \>CM000594.1 Fusarium oxysporum f. sp. lycopersici 4287 chromosome 6, whole genome shotgun sequence  \>CM000595.1 Fusarium oxysporum f. sp. lycopersici 4287 chromosome 7, whole genome shotgun sequence  \>CM000596.1 Fusarium oxysporum f. sp. lycopersici 4287 chromosome 8, whole genome shotgun sequence  \>CM000597.1 Fusarium oxysporum f. sp. lycopersici 4287 chromosome 9, whole genome shotgun sequence  \>CM000598.1 Fusarium oxysporum f. sp. lycopersici 4287 chromosome 10, whole genome shotgun sequence  \>CM000599.1 Fusarium oxysporum f. sp. lycopersici 4287 chromosome 11, whole genome shotgun sequence  \>CM000600.1 Fusarium oxysporum f. sp. lycopersici 4287 chromosome 12, whole genome shotgun sequence  \>CM000601.1 Fusarium oxysporum f. sp. lycopersici 4287 chromosome 13, whole genome shotgun sequence  \>CM000602.2 Fusarium oxysporum f. sp. lycopersici 4287 chromosome 14, whole genome shotgun sequence  \>CM000603.1 Fusarium oxysporum f. sp. lycopersici 4287 chromosome 15, whole genome shotgun sequence   Note that chromosome 14 ends with `.2`. This is because the assembly was updated and has a second version for this chromosome.Now do the same for F. graminearum: (pick the one with accession GCF_000240135.3, sequenced and assembled by the Broad institute)(https://www.ncbi.nlm.nih.gov/assembly/GCF_000240135.3/)WGS project, select accession starting with CM, send to file, etc, etc. Rename the file to `Fg_PH1.fasta`. Store it in `CircosPlot_FgFo/data/genomes/` `grep '>' Fg_PH1.fasta`.  \>CM000574.1 Fusarium graminearum PH-1 chromosome 1, whole genome shotgun sequence  \>CM000575.1 Fusarium graminearum PH-1 chromosome 2, whole genome shotgun sequence  \>CM000576.1 Fusarium graminearum PH-1 chromosome 3, whole genome shotgun sequence  \>CM000577.1 Fusarium graminearum PH-1 chromosome 4, whole genome shotgun sequence     2. Align the genomes of _Fusarium oxysporum_ and _Fusarium graminearum_  Next, we want to compare these genomes, to see which genomic regions are homologous to eachother. For this, we will use [MUMmer](http://mummer.sourceforge.net/). MUMmer is a package that contains programs to align genomes, filter aligned regions and plot algnments in dotplots. MUMmer is well-documented.We will use promer to align _Fo_ and _Fg_. Promer first translates dna sequences into protein and then aligns those. It is more sensitive over longer evolutionary distances than nucmer (that aligns DNA sequences directly).  I installed MUMmer version 3.23 in my /Applications folder so I will run the genome alignment like this:`/Applications/MUMmer3.23/promer`  I will add prefix (using `-p`) to my output files, so they will be saved in the correct directory and they will start with this name. As I align Fol4287 with Fg_PH1 using promer with default settings, I choose:`-p results/genome_alignments/Fol4287__vs__Fg_PH1.promer_def`Then we need to tell promer which genomes to align:`~data/genomes/Fol4287.fasta` and `data/genomes/FG_PH1.fasta`.  and we'll redirect (`>&`: `>` is to redirect output that would normally be printed in the terminal (also known as `stdout`, I add `&` to also include errormessages) the output of this to a logfile, which we will store in `docs`:`docs/Fol4287__vs__Fg_PH1.promer_def.log`.  The entire command is:`/Applications/MUMmer3.23/promer -p results/genome_alignments/Fol4287__vs__Fg_PH1.promer_def data/genomes/Fol4287.fasta data/genomes/FG_PH1.fasta >& docs/Fol4287__vs__Fg_PH1.promer_def.log`.This takes a few minutes...When it is finished, you have a file called `Fol4287__vs__Fg_PH1.promer_def.delta` in the directory `results/genome_alignments/`. Open this file and have a look at it's contents.  It should look something like this:  `PROMER  \>CM000591.1 CM000574.1 5634991 11697295  173055 174734 11476307 11474640 200 129 2  412  1  1  1  0  5329501 5328974 2528 3052 65 45 14  `  We won't go further into any detail, but notice that this is not an intuitive format to read. Luckily MUMmer also has a more manageable format, and we can obtain this by running `show-coords`. We add the `-r` option, to sort the output by reference:  `path/to/MUMmer/show-coords -r /results/genome_alignments/Fol4287__vs__Fg_PH1.promer_def.delta > results/genome_alignments/Fol4287__vs__Fg_PH1.promer_def.coords`  >We will visualize this alignment using Circos, but MUMmer also has a method to quickly display alignments, called MUMmerplot. This needs gnuplot, a plotting program that you can install with `conda`. We won't go into further detail for now but I want show you what this plot looks like:  >   >NB I modified the .gp a little to get this exact plot, I made the points a little smaller and removed the mouse-over options as they raised an error...    >Take a look at this plot, so you a little what to expect.  You can see for example that _Fg_ is on the y-axis (this matches with how we called `nucmer`, we first put the fasta of _Fo_ and then of _Fg_: _Fo_ is considered the `reference` (`REF` on the x-axis) and _Fg_ the `query` (`QRY` on the y-axis).  You can see that chromosome 1 (CM000574.1) in _Fg_ corresponds to chromosome 1 (CM000589.1) and chromosome 7 (CM000595.1) in _Fo_.  We can also see that _Fo_ has a few chromosomes that do not resemble any chromsome in _Fg_: chromosome 3, (CM000591.1) 6 (CM000594.1), 14 (CM000602.2) and 15 (CM000603.1) have very few dots and no lines in this plot.  Now let's visualize this in Circos!   3. Use Circos to visualize homologous regionsAs you may have noticed, _F. graminearum_ has 4 chromosomes where _F. oxysporum_ has 15. This means that chromosomes have rearranged, some may merged in _Fg_, some may have split in _Fo_ or more complicated scenario's of breaking, merging and mixing may have occurred. We already get a nice idea from the dotplot of what is going on, let's visualzie in Circos for perhaps a more clear and more attractive view.Circos uses a config file (ends with `.conf`) that contains information on how a plot should look. First thing we should add to the .conf file is the karyotype: the number and size of chromosomes. This is stored in a karyotype file. 3.1 Karyotype fileThe karyotype file also contain information on the colors of the chromosomes:  >chr - ID LABEL START END COLORSo for example, chr 1 in _Fg_ is `CM000574.1`, and it has 11697295 basepairs and if we want plot these in orange we write:  `chr - CM000574.1 1 0 11697295 orange`    Let's create a such a file using Python. We will plot the chromosomes of different species in different colors in order to separate them.<jupyter_code># first karyotype file
import glob

# I run this notebook from the working directory hence I don't need to
# specify that path; instead of path/to/CircosPlot_FgFo/data/genomes/
# I can just say data/genomes/
indir   = 'data/genomes/' 
outdir  = 'circos/input/karyotypes/'
def karyotype_file_from_fasta(fasta_fname, out_fname, color = 'grey'):
    '''
    This function takes a fastafile (fasta_fname) and "translates" this 
    to a karyotype file (out_fname). Grey is the default color for chromosomes.
    '''
    outfile = open(out_fname, 'w')
    l      = 0 # counter for the sequence length
    cindex = 1 # counter for the chromosomes, note that this only works 
               # because the fastafile is ordered!
    cid    = ''
    for line in open(fasta_fname).readlines():
        if line[0] == '>': 
            '''
            This is the fastaheader. We need to store the name in this header 
            (we assign it to the variable cid), while we count the length of 
            the subsequent lines (the sequence of this chromosome), until we 
            encounter the next header
            '''
            if l > 0: 
                '''
                If l > 0 it means it is not the first sequence.
                We first have to process the previous sequence before we can move on 
                with this one
                '''    
                #process previous sequence
                outfile.write('chr - '+cid+' '+str(cindex)+' 0 '+' '+str(l)+' '+color+'\n')
                #make everything ready for this one
                l = 0
                cindex += 1
                
            '''If you encounter a header, store it in cid.
            '''
            cid = line[1:].strip().split()[0]
            print(cid, cindex)
                
        else:
            ''' 
            If it is not header it is part of the sequence, add the length of this line. 
            str.strip() takes of any newlines, spaces etc. so we don't count those.
            '''
            l += len(line.strip())
        
    ''' 
    When you've reached the end of the file, you still need to process the last header. 
    We do this here.
    '''
    if l > 0:
        outfile.write('chr - '+cid+' '+str(cindex)+' 0 '+' '+str(l)+' '+color+'\n')
            
    outfile.close()
    
colors = ['green', 'orange']
filenames = glob.glob(indir+'*.fasta')

if len(colors) != len(filenames): #check if your list of colors is long enough 
    print("I need more colors, I have ", len(colors), " color(s) for ", len(filenames), " files")
else:  
    for i, fasta_fname in enumerate(glob.glob(indir+'*.fasta')):
        print(fasta_fname.split('/')[-1].split('.fasta')[0],colors[i])
        out_fname = outdir + fasta_fname.split('/')[-1].split('.fasta')[0]+'.txt'
        karyotype_file_from_fasta(fasta_fname, out_fname, color = colors[i])<jupyter_output>Fol4287 green
CM000589.1 1
CM000590.1 2
CM000591.1 3
CM000592.1 4
CM000593.1 5
CM000594.1 6
CM000595.1 7
CM000596.1 8
CM000597.1 9
CM000598.1 10
CM000599.1 11
CM000600.1 12
CM000601.1 13
CM000602.2 14
CM000603.1 15
Fg_PH1 orange
CM000574.1 1
CM000575.1 2
CM000576.1 3
CM000577.1 4<jupyter_text>3.2 Config file 3.2.1 Karyotypes and labelsLet's see if we can plot these chromosomes in Circos. Circos has so many possibilities that we can not simply define all our wishes on the commandline. They're to many. Circos stores parameters settings in config files. There is one main config file, that can import/load/inlcude other config files. The minimum config file can be found [here](http://circos.ca/documentation/tutorials/quick_start/hello_world/).  (I added line numbers for convenience).      1:   circos.conf    2:      3:  karyotype = data/karyotype/karyotype.human.txt    4:      5:      6:      7:      8:  default = 0.005r    9:      10:      11: radius    = 0.9r      12: thickness = 20p      13: fill      = yes      14:      15:     16:        17:         18:  The remaining content is standard and required. It is imported       19:  from default files in the Circos distribution.      20:       21:  These should be present in every Circos configuration file and      22:  overridden as required. To see the content of these files,       23:  look in etc/ in the Circos distribution.      24:       25:       26:  Included from Circos distribution.      27: >      28:       29:       30:  RGB/HSV color definitions, color lists, location of fonts, fill patterns.      31:  Included from Circos distribution.      32: >      33:         34:  Debugging, I/O an dother system parameters      35:  Included from Circos distribution.      36: >       Open a text editor and copy-paste the contents of a minimum config file. Save it with a sensible filename such as Fo_vs_Fg.circos.conf.   First we change the first line, which contains the name of this file.Then we add the correct filename for the karyotype file (line 3). You can plot multiple karyotype files into a single plot by separating them with a ,.You can use radius (line 11) and thickness (line 12) to lot the chromosome in smaller circle or with thicker bands respectively. We also want to add labels to our chromosomes, so we know what we look at.To add labels, you have to add `show_labels = yes` to the configfile and you need to specify where (`label_radius   = dims(ideogram,radius_outer) + 10p`) and how (`label_font = default`, `label_size = 24p`, `label_parallel = yes`) to show the labels. After this, the file becomes:    1:   Fo_vs_Fg.circos.conf.    2:      3:  karyotype = circos/input/karyotypes/Fg_PH1.karyotype.txt, circos/input/karyotypes/Fol4287.karyotype.txt    4:      5:      6:      7:      8:  default = 0.005r    9:      10:      11: radius    = 0.9r      12: thickness = 20p      13: fill      = yes      14:     15: show_label     = yes    16: label_radius   = dims(ideogram,radius_outer) + 10p    17: label_font     = default    18: label_size     = 24p    19: label_parallel = yes    20:     21:        22:         23:  The remaining content is standard and required. It is imported       24:  from default files in the Circos distribution.      25:       26:  These should be present in every Circos configuration file and      27:  overridden as required. To see the content of these files,       28:  look in etc/ in the Circos distribution.      29:       30:       31:  Included from Circos distribution.      32: >      33:       34:       35:  RGB/HSV color definitions, color lists, location of fonts, fill patterns.      36:  Included from Circos distribution.      37: >      38:         39:  Debugging, I/O an dother system parameters      40:  Included from Circos distribution.      41: >  Let's try!   We tell Circos to use this config file with `-conf your/configfile`.I have installed circos is the `/Applications` folder (MacOS), change this to your Circos installation folder.<jupyter_code>%%bash
/Applications/circos-0.69-6/bin/circos -conf circos/input/config_files/Fo_vs_Fg.circos.conf<jupyter_output>debuggroup summary 0.17s welcome to circos v0.69-6 31 July 2017 on Perl 5.018002
debuggroup summary 0.18s current working directory /Users/like/Dropbox/00.Projects/CircosPlot_FgFo
debuggroup summary 0.18s command /Applications/circos-0.69-6/bin/circos -conf circos/input/config_files/Fo_vs_Fg.circos.conf
debuggroup summary 0.18s loading configuration from file circos/input/config_files/Fo_vs_Fg.circos.conf
debuggroup summary 0.18s found conf file circos/input/config_files/Fo_vs_Fg.circos.conf
debuggroup summary 0.33s debug will appear for these features: output,summary
debuggroup summary 0.33s bitmap output image ./circos.png
debuggroup summary 0.33s SVG output image ./circos.svg
debuggroup summary 0.33s parsing karyotype and organizing ideograms
debuggroup summary 0.34s karyotype has 19 chromosomes of total size 94,271,448
debuggroup summary 0.34s applying global and local scaling
debuggroup summary 0.36s allocating image, colors and brushes
debuggroup summary 7.47s drawing 19 ideogram[...]<jupyter_text>Everything seemed to have worked well, now where is the plot?From the last lines of the output:      debuggroup output 7.53s created PNG image ./circos.png (106 kb)      debuggroup output 7.53s created SVG image ./circos.svg (9 kb)     We can see that Circos saved the figures in the folder it is run from. Let's have a look: We can change this either in the configfile (see an other notebook Circos_output_config_files on how to do this and understand better how Circos uses config files) or on the commandline. Let's do that:<jupyter_code>%%bash
/Applications/circos-0.69-6/bin/circos -conf circos/input/config_files/Fo_vs_Fg.circos.conf -outputdir circos/output -outputfile circos_karyotype_with_labels.png<jupyter_output>debuggroup summary 0.24s welcome to circos v0.69-6 31 July 2017 on Perl 5.018002
debuggroup summary 0.24s current working directory /Users/like/Dropbox/00.Projects/CircosPlot_FgFo
debuggroup summary 0.24s command /Applications/circos-0.69-6/bin/circos -conf circos/input/config_files/Fo_vs_Fg.circos.conf -outputdir circos/output -outputfile circos_karyotype_with_labels.png
debuggroup summary 0.24s loading configuration from file circos/input/config_files/Fo_vs_Fg.circos.conf
debuggroup summary 0.24s found conf file circos/input/config_files/Fo_vs_Fg.circos.conf
debuggroup summary 0.39s debug will appear for these features: output,summary
debuggroup summary 0.39s bitmap output image circos/output/circos_karyotype_with_labels.png
debuggroup summary 0.39s SVG output image circos/output/circos_karyotype_with_labels.svg
debuggroup summary 0.39s parsing karyotype and organizing ideograms
debuggroup summary 0.39s karyotype has 19 chromosomes of total size 94,271,448
debuggroup summary 0.40s ap[...]<jupyter_text>3.3 Links 3.3.1 Write links file from MUMmer coords fileWe now have plotted the chromosomes of _Fo_ and _Fg_, but we want to see whther and how different chromosomes in these two species correspond to eachother. Circos is desiged to show links between ordered sets, so let's tart using it what is was designed to do.    Link files have the [following format](http://circos.ca/documentation/tutorials/links/basic_links/lesson):  >chrA from to chrB from toWe can get this information from the coords file:<jupyter_code>lines   = open('results/genome_alignments/Fol4287__vs__Fg_PH1.promer_def.coords').readlines()[5:]
print(lines[0])
columns = lines[0].split()
print(columns)<jupyter_output>9030     9248  |  7107985  7108206  |      219      222  |    63.51    68.92     0.68  |  3  1  CM000589.1	CM000574.1

['9030', '9248', '|', '7107985', '7108206', '|', '219', '222', '|', '63.51', '68.92', '0.68', '|', '3', '1', 'CM000589.1', 'CM000574.1']<jupyter_text>Now let's write the 'translation' from this format to the circos link format to a function<jupyter_code>def coords_to_links(coords_fname, links_fname):
    linksfile = open(links_fname, 'w')
    lines     = open(coords_fname).readlines()[5:] #skip header
    for line in lines:
        columns = line.strip().split()
        
        if len(columns) != 17:
            print('Unexpected fileformat, can not continue!')
            break
            
        chrA  = columns[15]
        chrB  = columns[16]
        fromA = columns[0]
        toA   = columns[1]
        fromB = columns[3]
        toB   = columns[4]
        
        linksfile.write(chrA+' '+fromA+' '+toA+' '+chrB+' '+fromB+' '+toB+'\n')
    linksfile.close()
    print(coords_fname,": Done")
    
coords_fname ='results/genome_alignments/Fol4287__vs__Fg_PH1.promer_def.coords'
links_fname  = 'circos/input/links/links_Fo_vs_Fg.txt'
coords_to_links(coords_fname, links_fname)<jupyter_output>results/genome_alignments/Fol4287__vs__Fg_PH1.promer_def.coords : Done<jupyter_text>Let's have a look:<jupyter_code>%%bash
less circos/input/links/links_Fo_vs_Fg.txt | head -10<jupyter_output><empty_output><jupyter_text>3.3.2 Include links in config fileLet's add this to our conf file (save it under a new name):     Fo_vs_Fg.circos.incl_links.conf        karyotype = karyotypes/Fg_PH1.txt,karyotypes/Fol4287.txt                default = 0.005r        position    radius    	= 0.80r    thickness 	= 50p    fill     	= yes    label    show_label     = yes    label_font     = default    label_radius   = dims(ideogram,radius_outer) + 10p    label_size     = 24p    label_parallel = yes        **Here we can add a `links` section:**                file          = input/links/links_Fo_vs_Fg.txt    radius        = 0.95r    color         = black_a4     Curves look best when this value is small (e.g. 0.1r or 0r)    bezier_radius = 0.1r    thickness     = 1                 The remaining content is standard and required. It is imported      from default files in the Circos distribution.         These should be present in every Circos configuration file and     overridden as required. To see the content of these files,      look in etc/ in the Circos distribution.         Included from Circos distribution.    >         RGB/HSV color definitions, color lists, location of fonts, fill patterns.     Included from Circos distribution.    >     Debugging, I/O an dother system parameters     Included from Circos distribution.    ><jupyter_code>%%bash
/Applications/circos-0.69-6/bin/circos -conf circos/input/config_files/Fo_vs_Fg.circos.incl_links.conf  -outputdir circos/output -outputfile circos_karyotype_with_labels_and_links.png<jupyter_output>debuggroup summary 0.14s welcome to circos v0.69-6 31 July 2017 on Perl 5.018002
debuggroup summary 0.14s current working directory /Users/like/Dropbox/00.Projects/CircosPlot_FgFo
debuggroup summary 0.14s command /Applications/circos-0.69-6/bin/circos -conf circos/input/config_files/Fo_vs_Fg.circos.incl_links.conf -outputdir circos/output -outputfile circos_karyotype_with_labels_and_links.png
debuggroup summary 0.14s loading configuration from file circos/input/config_files/Fo_vs_Fg.circos.incl_links.conf
debuggroup summary 0.14s found conf file circos/input/config_files/Fo_vs_Fg.circos.incl_links.conf
debuggroup summary 0.37s debug will appear for these features: output,summary
debuggroup summary 0.37s bitmap output image circos/output/circos_karyotype_with_labels_and_links.png
debuggroup summary 0.37s SVG output image circos/output/circos_karyotype_with_labels_and_links.svg
debuggroup summary 0.38s parsing karyotype and organizing ideograms
debuggroup summary 0.38s karyotype has 19 c[...]<jupyter_text>We get an error that we're trying to display to many links. We can either reset the maximum number of links, or we can filter the links we plot. Let's do the latter:<jupyter_code>def coords_to_links(coords_fname, links_fname, min_length = 0):
    linksfile = open(links_fname, 'w')
    lines     = open(coords_fname).readlines()[5:] #skip header
    for line in lines:
        columns = line.strip().split()
        
        if len(columns) != 17:
            print('Unexpected fileformat, can not continue!')
            break
            
        chrA  = columns[15]
        chrB  = columns[16]
        fromA = columns[0]
        toA   = columns[1]
        fromB = columns[3]
        toB   = columns[4]
        
        lenA  = int(columns[6])
        lenB  = int(columns[7])
        if lenA >= min_length and lenB >= min_length: 
            linksfile.write(chrA+' '+fromA+' '+toA+' '+chrB+' '+fromB+' '+toB+'\n')
    linksfile.close()

coords_fname = 'results/genome_alignments/Fol4287__vs__Fg_PH1.promer_def.coords'
min_length   = 1000
links_fname  = 'circos/input/links/links_Fo_vs_Fg.min_length'+str(min_length)+'bp.txt'
coords_to_links(coords_fname, links_fname, min_length = min_length)<jupyter_output><empty_output><jupyter_text>Here we filter links based on their length, only including alignments that span 1000 bp or more in both genomes. ALternatively, you can also filter based on sequence similarity or both.Let's change our config file to display the filtered links only (again saving the file under a new name):     Fo_vs_Fg.circos.incl_links.filtered.conf        karyotype = input/karyotypes/Fg_PH1.txt, input/karyotypes/Fol4287.txt                default = 0.005r        position    radius    	= 0.80r    thickness 	= 50p    fill     	= yes    label    show_label     = yes    label_font     = default    label_radius   = dims(ideogram,radius_outer) + 10p    label_size     = 24p    label_parallel = yes                    **Change the file in the `links` section:**    file          = input/links/links_Fo_vs_Fg.min_length1000bp.txt    radius        = 0.95r    color         = black_a4     Curves look best when this value is small (e.g. 0.1r or 0r)    bezier_radius = 0.1r    thickness     = 1                 The remaining content is standard and required. It is imported      from default files in the Circos distribution.         These should be present in every Circos configuration file and     overridden as required. To see the content of these files,      look in etc/ in the Circos distribution.         Included from Circos distribution.    dir  = /path/to/your/output/directory    file = yourimage.png    >         RGB/HSV color definitions, color lists, location of fonts, fill patterns.     Included from Circos distribution.    >     Debugging, I/O an dother system parameters     Included from Circos distribution.    >    Let's try and run this (this may take a while, be a little patient)<jupyter_code>%%bash
/Applications/circos-0.69-6/bin/circos -conf circos/input/config_files/Fo_vs_Fg.circos.incl_links.filtered.conf  -outputdir circos/output -outputfile circos_karyotype_with_labels_and_links.png<jupyter_output>debuggroup summary 0.15s welcome to circos v0.69-6 31 July 2017 on Perl 5.018002
debuggroup summary 0.15s current working directory /Users/like/Dropbox/00.Projects/CircosPlot_FgFo
debuggroup summary 0.15s command /Applications/circos-0.69-6/bin/circos -conf circos/input/config_files/Fo_vs_Fg.circos.incl_links.filtered.conf -outputdir circos/output -outputfile circos_karyotype_with_labels_and_links.png
debuggroup summary 0.15s loading configuration from file circos/input/config_files/Fo_vs_Fg.circos.incl_links.filtered.conf
debuggroup summary 0.15s found conf file circos/input/config_files/Fo_vs_Fg.circos.incl_links.filtered.conf
debuggroup summary 0.31s debug will appear for these features: output,summary
debuggroup summary 0.32s bitmap output image circos/output/circos_karyotype_with_labels_and_links.png
debuggroup summary 0.32s SVG output image circos/output/circos_karyotype_with_labels_and_links.svg
debuggroup summary 0.32s parsing karyotype and organizing ideograms
debuggroup summa[...]<jupyter_text><jupyter_code>import glob

def coords_to_links(coords_fname, links_fname, chr2color = None, min_length = 0):
    linksfile = open(links_fname, 'w')
    lines     = open(coords_fname).readlines()[5:] #skip header
    for line in lines:
        columns = line.strip().split()
        
        if len(columns) != 17:
            print('Unexpected fileformat, can not continue!')
            break
            
        chrA  = columns[15]
        chrB  = columns[16]
        fromA = columns[0]
        toA   = columns[1]
        fromB = columns[3]
        toB   = columns[4]
        
        lenA  = int(columns[6])
        lenB  = int(columns[7])
        if lenA >= min_length and lenB >= min_length: 
            out = chrA+' '+fromA+' '+toA+' '+chrB+' '+fromB+' '+toB
            
            if chr2color != None:
                if chrA in chr2color:
                    linksfile.write(out+' color='+chr2color[chrA]+'\n')
                elif chrB in chr2color:
                    linksfile.write(out+' color='+chr2color[chrB]+'\n')
            else:
                linksfile.write(out+'\n')
    linksfile.close()


def karyotype_file_from_fasta(fasta_fname, out_fname, chr2color = None):

    outfile = open(out_fname, 'w')
    l      = 0 #counter for the sequence length
    cindex = 1 #counter for the chromosomes
    cid    = ''
    color  = 'grey'
    for line in open(fasta_fname).readlines():
        if line[0] == '>':
            if l > 0:
                outfile.write('chr - '+cid+' '+str(cindex)+' 0 '+' '+str(l)+' '+color+'\n')
                l = 0
                cindex += 1
            cid = line[1:].strip().split()[0]
            if chr2color != None and cid in chr2color:
                color = chr2color[cid]
            print(cid, cindex)
                
        else:
            l += len(line.strip())
        
    #last one:
    if l > 0:
        outfile.write('chr - '+cid+' '+str(cindex)+' 0 '+' '+str(l)+' '+color+'\n')
            
    outfile.close()
    
# write karyotype files from fasta
chr2color = {}
chr2color['CM000574.1'] = 'blue'
chr2color['CM000575.1'] = 'yellow'
chr2color['CM000576.1'] = 'green'
chr2color['CM000577.1'] = 'red'

for i, fasta_fname in enumerate(glob.glob('data/genomes/*.fasta')):
    
    out_fname = 'circos/input/karyotypes/' + fasta_fname.split('/')[-1].split('.fasta')[0]+'.color-coded.txt'
    karyotype_file_from_fasta(fasta_fname, out_fname, chr2color = chr2color)

coords_fname = 'results/genome_alignments/Fol4287__vs__Fg_PH1.promer_def.coords'
links_fname  = 'circos/input/links/links_Fo_vs_Fg.min_length1000bp.color-coded.txt'
coords_to_links(coords_fname, links_fname, chr2color = chr2color, min_length = 1000)<jupyter_output>CM000589.1 1
CM000590.1 2
CM000591.1 3
CM000592.1 4
CM000593.1 5
CM000594.1 6
CM000595.1 7
CM000596.1 8
CM000597.1 9
CM000598.1 10
CM000599.1 11
CM000600.1 12
CM000601.1 13
CM000602.2 14
CM000603.1 15
CM000574.1 1
CM000575.1 2
CM000576.1 3
CM000577.1 4<jupyter_text>Now adjust the config file, to incorporate these karyotype and links files:         Fo_vs_Fg.circos.incl_colored_links.filtered.conf        karyotype = input/karyotypes/Fg_PH1.color-coded.txt, input/karyotypes/Fol4287.color-coded.txt                default = 0.005r        position    radius    	= 0.80r    thickness 	= 50p    fill     	= yes    label    show_label     = yes    label_font     = default    label_radius   = dims(ideogram,radius_outer) + 10p    label_size     = 24p    label_parallel = yes                    **Change the file in the `links` section:**    file          = input/links/links_Fo_vs_Fg.min_length1000bp.color-coded.txt    radius        = 0.95r    color         = black_a4     Curves look best when this value is small (e.g. 0.1r or 0r)    bezier_radius = 0.1r    thickness     = 1                 The remaining content is standard and required. It is imported      from default files in the Circos distribution.         These should be present in every Circos configuration file and     overridden as required. To see the content of these files,      look in etc/ in the Circos distribution.         Included from Circos distribution.    dir  = /path/to/your/output/directory    file = yourimage.png    >         RGB/HSV color definitions, color lists, location of fonts, fill patterns.     Included from Circos distribution.    >     Debugging, I/O an dother system parameters     Included from Circos distribution.    >    Let's try and run this (this may take a while, be a little patient)<jupyter_code>%%bash
/Applications/circos-0.69-6/bin/circos -conf circos/input/config_files/Fo_vs_Fg.circos.incl_colored_links.filtered.conf  -outputdir circos/output -outputfile circos_karyotype_with_labels_and_colored_links.png<jupyter_output>debuggroup summary 0.12s welcome to circos v0.69-6 31 July 2017 on Perl 5.018002
debuggroup summary 0.13s current working directory /Users/like/Dropbox/00.Projects/CircosPlot_FgFo
debuggroup summary 0.13s command /Applications/circos-0.69-6/bin/circos -conf circos/input/config_files/Fo_vs_Fg.circos.incl_colored_links.filtered.conf -outputdir circos/output -outputfile circos_karyotype_with_labels_and_colored_links.png
debuggroup summary 0.13s loading configuration from file circos/input/config_files/Fo_vs_Fg.circos.incl_colored_links.filtered.conf
debuggroup summary 0.13s found conf file circos/input/config_files/Fo_vs_Fg.circos.incl_colored_links.filtered.conf
debuggroup summary 0.27s debug will appear for these features: output,summary
debuggroup summary 0.28s bitmap output image circos/output/circos_karyotype_with_labels_and_colored_links.png
debuggroup summary 0.28s SVG output image circos/output/circos_karyotype_with_labels_and_colored_links.svg
debuggroup summary 0.28s parsing kar[...]<jupyter_text>This is getting somewhere! We can see that chromosom1 in _Fg_ corresponds roughly to chr 7, half of chr 4 and chr 1 in _Fo_, but e.g. for chr 02 in _Fg_ it is tricky to see which part belongs to which chromosome in _Fo_. We can further improve our visualization by color-coding the chromosomes in _Fo_, rather than the chromosomes in _Fg_ and make 4 plots: one per chromosome.Circos implements [brewer palettes for qualitative data] (http://colorbrewer2.org/type=qualitative&scheme=Set3&n=12) for up to 12 categories. _Fo_ has 15 chromosomes, but chromosome 3, 6, 14 and 15 have very little regions that are homologous to regions in _Fg_ (are lineage-specific), hence we use a single color (e.g. black) for those. If we use `set3` of 11 colors (15 - 4 lineage-specific chromosomes) corresponding to categories (rather e.g. sequential data):`qual`, we can refer to them with `set3-11-qual-N`, where N is the number 1 to 11.<jupyter_code>import glob, os

# Write karyotype files from fasta
chr2color = {}
chr2color['CM000589.1'] = 'set3-11-qual-1'
chr2color['CM000590.1'] = 'set3-11-qual-2'
chr2color['CM000591.1'] = 'light-grey' #chr 3: lineage-specific
chr2color['CM000592.1'] = 'set3-11-qual-3'
chr2color['CM000593.1'] = 'set3-11-qual-4'
chr2color['CM000594.1'] = 'light-grey' #chr 6: lineage-specific
chr2color['CM000595.1'] = 'set3-11-qual-5'
chr2color['CM000596.1'] = 'set3-11-qual-6'
chr2color['CM000597.1'] = 'set3-11-qual-7'
chr2color['CM000598.1'] = 'set3-11-qual-8'
chr2color['CM000599.1'] = 'set3-11-qual-9'
chr2color['CM000600.1'] = 'set3-11-qual-10'
chr2color['CM000601.1'] = 'set3-11-qual-11'
chr2color['CM000602.2'] = 'light-grey' #chr 14: lineage-specific
chr2color['CM000603.1'] = 'light-grey' #chr 15: lineage-specific


# Write the karyotype file for Fol4287
Fol_karyo_fname = 'circos/input/karyotypes/Fol4287_chr_color-coded.txt'
fasta_fname     = 'data/genomes/Fol4287.fasta'
karyotype_file_from_fasta(fasta_fname, Fol_karyo_fname, chr2color = chr2color)

#write the links file, now color them according to the Fol4287 chromosomes instead of the Fg_PH1 chromosomes
coords_fname = 'results/genome_alignments/Fol4287__vs__Fg_PH1.promer_def.coords'
links_fname  = 'circos/input/links/links_Fo_vs_Fg.min_length1000bp.color-coded_Fol.txt'
coords_to_links(coords_fname, links_fname, chr2color = chr2color, min_length = 1000)

# We need to make 4 plots, one for each Fg chromosome.
# There are a lot of different ways to do this. 
#
# I decided here to write a separate file for each chromosome in F. graminearum, 
# splitting up the existing karyotype file.
#
# Moreover, as I am looping over the chromosomes in the karyotypefile, 
# I also write the config file.
#
# First open a 'template' config file and store it's lines, so we can write them to 
# the 4 new chromosome-specific config files:

config_lines = open('circos/input/config_files/Fo_vs_Fg.circos.incl_colored_links.filtered.conf').readlines()

'''
You can run the code below, to print the line number and line of the config file to make sure 
you're changing the correct line:

for i,line in enumerate(config_lines):
    print(i,line.strip())
    
'''

# loop over the Fg chromosomes:
karyotype_file_Fg = 'circos/input/karyotypes/Fg_PH1.txt'
for i, line in enumerate(open(karyotype_file_Fg).readlines()):
    
    # write the karyotype file
    chrom_spec_fname = 'circos/input/karyotypes/Fg_chr'+str(i+1)+'.txt'
    Fg_chrom_file = open(chrom_spec_fname, 'w')
    Fg_chrom_file.write(line.replace('orange', 'grey')) # change color of the chromosomes 
                                                        # (one of the Fo chrs is orange)
    Fg_chrom_file.close()
    
    # Make the config file:
    config_fname         = 'circos/input/config_files/Fo_vs_Fg_chr'+str(i+1)\
                                +'.circos.incl_colored_links.filtered.conf'
    # change the line with the name of the config file
    config_lines[0] = '# '+config_fname+'\n'
    # change the line listing the karyotype files, to use the chr-specific karyotype file for Fg
    config_lines[2] = 'karyotype = '+chrom_spec_fname+','+Fol_karyo_fname+'\n'
    # change the line with the file to use for the links
    config_lines[31] = 'file\t= '+links_fname+'\n'
    
    # Write the config file
    Fg_chrom_config_file = open(config_fname, 'w')
    for line in config_lines:
        Fg_chrom_config_file.write(line)
    Fg_chrom_config_file.close()
    
    # Finally, run circos using this config file:
    circos_command = '/Applications/circos-0.69-6/bin/circos' # path to Circos
    circos_command +=' -conf '+config_fname                   # the config file we want circos to use
    circos_command +=' -outputdir circos/output'              # where to save the output
    # filename of the generated plot
    circos_command +=' -outputfile circos_Fg_chr'+str(i+1)+'.vs.Fol4287_with_labels_and_colored_links.png' 

    print(circos_command)
    os.system(circos_command)<jupyter_output>CM000589.1 1
CM000590.1 2
CM000591.1 3
CM000592.1 4
CM000593.1 5
CM000594.1 6
CM000595.1 7
CM000596.1 8
CM000597.1 9
CM000598.1 10
CM000599.1 11
CM000600.1 12
CM000601.1 13
CM000602.2 14
CM000603.1 15
/Applications/circos-0.69-6/bin/circos -conf circos/input/config_files/Fo_vs_Fg_chr1.circos.incl_colored_links.filtered.conf -outputdir circos/output -outputfile circos_Fg_chr1.vs.Fol4287_with_labels_and_colored_links.png
/Applications/circos-0.69-6/bin/circos -conf circos/input/config_files/Fo_vs_Fg_chr2.circos.incl_colored_links.filtered.conf -outputdir circos/output -outputfile circos_Fg_chr2.vs.Fol4287_with_labels_and_colored_links.png
/Applications/circos-0.69-6/bin/circos -conf circos/input/config_files/Fo_vs_Fg_chr3.circos.incl_colored_links.filtered.conf -outputdir circos/output -outputfile circos_Fg_chr3.vs.Fol4287_with_labels_and_colored_links.png
/Applications/circos-0.69-6/bin/circos -conf circos/input/config_files/Fo_vs_Fg_chr4.circos.incl_colored_links.filtered.conf -outputd[...]<jupyter_text>Note on versioning  During this workshop/tutorial, we've been adjusting our code and circos input and outfiles, constantly improving them.This resulted in the following output:<jupyter_code>%%bash
echo '     circos/input/config_files : '
ls circos/input/config_files
echo '     circos/input/karyotypes : '
ls circos/input/karyotypes
echo '     circos/input/links : '
ls circos/input/links
echo '     circos/output :'
ls circos/output<jupyter_output>circos/input/config_files : 
Fo_vs_Fg.circos.conf
Fo_vs_Fg.circos.incl_colored_links.filtered.conf
Fo_vs_Fg.circos.incl_links.conf
Fo_vs_Fg.circos.incl_links.filtered.conf
Fo_vs_Fg.circos.specify_output.conf
Fo_vs_Fg_chr1.circos.incl_colored_links.filtered.conf
Fo_vs_Fg_chr2.circos.incl_colored_links.filtered.conf
Fo_vs_Fg_chr3.circos.incl_colored_links.filtered.conf
Fo_vs_Fg_chr4.circos.incl_colored_links.filtered.conf
image.Fo_vs_Fg.conf
image.generic.Fo_vs_Fg.conf
     circos/input/karyotypes : 
Fg_PH1.color-coded.txt
Fg_PH1.txt
Fg_chr0.txt
Fg_chr1.txt
Fg_chr2.txt
Fg_chr3.txt
Fg_chr4.txt
Fol4287.color-coded.txt
Fol4287.txt
Fol4287_chr_color-coded.txt
     circos/input/links : 
links_Fo_vs_Fg.min_length1000bp.color-coded.txt
links_Fo_vs_Fg.min_length1000bp.color-coded_Fol.txt
links_Fo_vs_Fg.min_length1000bp.txt
links_Fo_vs_Fg.txt
     circos/output :
circos_Fg_chr1.vs.Fol4287_with_labels_and_colored_links.png
circos_Fg_chr1.vs.Fol4287_with_labels_and_colored_links.svg
circos_Fg_[...]
{'fraction_non_alphanumeric': 0.08493276634869555, 'fraction_numerical': 0.051304447764624755, 'mean_word_length': 3.530931056029994, 'pattern_counts': {'":': 1, '<': 42, '<?xml version=': 0, '>': 100, 'https://': 6, 'lorem ipsum': 0, 'www.': 6, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '408521', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 15754, 'n_tokens_neox': 13822, 'n_words': 4623, 'path': 'lessons/20180403_Comparative_genomics_with_Circos_Like/Circos_demo_03-04-2018.ipynb', 'repo_name': 'JoachimGoedhart/studyGroup', 'size': '43505'}