In [ ]:
from random import gauss
def illumina_qualgen (m_start, m_mid, m_end, sd_start, sd_mid, sd_end, length):
start_border = int(0.2*length)
mid_border = int(0.66*length)
end_border = length
mean_qual = []
sd_qual =[]
for i in range (start_border) :
mean_qual.append (int(m_start + (float (m_mid - m_start))/start_border*i))
sd_qual.append (int(sd_start + (float (sd_mid - sd_start))/start_border*i))
for i in range (mid_border - start_border) :
mean_qual.append(m_mid)
sd_qual.append(sd_mid)
for i in range(end_border - mid_border):
mean_qual.append (int(m_mid + (float (m_end - m_mid)/ (end_border - mid_border)*i)))
sd_qual.append (int(sd_mid + (float (sd_end - sd_mid)/ (end_border - mid_border)*i)))
return (mean_qual, sd_qual)
def random_qual_string (mean_qual, sd_qual, length):
return [valid_qual(mean_qual[i], sd_qual[i]) for i in range (length)]
def valid_qual(mean, sd):
while True:
qual = int(gauss(mean, sd))
if 0 <= qual <= 40:
return qual
def illumina_qualgen2 (len, Q0m, Q1m, Q2m, Q3m, Q4m, Q0r, Q1r, Q2r, Q3r, Q4r, Q0f, Q1f, Q2f, Q3f, Q4f):
"""Create a quality template for mean value and sd as """
l = [[Q0m, Q0m - Q0r/2, Q0m + Q0r/2, int(Q0f * len)],
[Q1m, Q1m - Q1r/2, Q1m + Q1r/2, int(Q1f * len)],
[Q2m, Q2m - Q2r/2, Q2m + Q2r/2, int(Q2f * len)],
[Q3m, Q3m - Q3r/2, Q3m + Q3r/2, int(Q3f * len)],
[Q4m, Q4m - Q4r/2, Q4m + Q4r/2, int(Q4f * len)]]
qual_pattern = []
for area in range (4):
start_mean = l [area] [0]
delta_mean = int(l [area + 1] [0] - l [area] [0])
start_min = l [area] [1]
delta_min = int(l [area + 1] [1] - l [area] [1])
start_max = l [area] [2]
delta_max = int(l [area + 1] [2] - l [area] [2])
len_area = l [area + 1] [3]
for pos in range(len_area):
mean_qual = int( start_mean + (float (delta_mean)/ (len_area)*pos))
min_qual = int( start_min + (float (delta_min) / (len_area)*pos))
max_qual = int( start_max + (float (delta_max) / (len_area)*pos))
qual_pattern.append([mean_qual, min_qual, max_qual])
return qual_pattern
Implementation d'une classe accessoire permettant de générer des listes de score de qualité. Un pattern modèle est déterminé à l'instanciation de la classe (mean, min et max pour chaque position) à partir de paramètres donnés par l'utilisateur pour le départ, le premier quart, le milieu, le troisième quart et la fin. Pour chaque zone il faut fournir la moyenne et la sd
A partir de ce pattern, l'objet peut être appelé pour générer des listes de score de qualité suivant suivant le pattern (facteur d'attraction à la moyenne) mais avec une une part stochastique suivant une distribution Gaussienne centrée sur de 0 en variant dans l'intervalle min, max autorisé
Ci dessous voici le script utilisé pour tester le programme.
In [ ]:
from QualGenerator import QualGenerator
q = QualGenerator(100, 30,35,37,35,20, 10,4,2,6,20, 0,15,25,25,35)
with (open("test.csv", "w")) as file:
for i in range (10):
s =""
for i in q.random_qual_string():
s += "{};".format(i)
file.write(s+"\n")
Quelques modifications ont été apporté par la suite pour ne plus travailler en intervalle de qualité mais plutôt en sd directement. D'autre part, j'ai simplifié considérablement les paramètres à passer pour instancier la classe. Il s'agit désormais d'indiquer si on souhaite une qualité : "very good", "good", "medium", "bad", "very bad". Les paramètres de chaque catégorie de qualité sont inclus en dur dans la classe d'après des approximations expérimentales
In [ ]:
qual_param_dict = {
"very_good" : [[33,38,39,38,36], [2,1,1,1,5], [0,a1,a2,a3,a4]],
"good" : [[32,37,38,37,32], [2,2,1,3,10], [0,a1,a2,a3,a4]],
"medium" : [[30,36,37,32,25], [4,3,2,10,10], [0,a1,a2,a3,a4]],
"bad" : [[25,30,32,28,20], [4,3,2,10,10], [0,a1,a2,a3,a4]],
"very_bad" : [[15,20,20,15,5], [4,3,2,10,10], [0,a1,a2,a3,a4]]}
Pour chaque entrée du dictionnaire la première liste correspond à la qualité moyenne, la seconde à la sd, et la troisième à la taille en pb de la région (recalculé pour chaque longueur donnée par l'utilisateur
Il reste potentiellement encore un peu de travail sur cette classe :
Reformatage du code pour se plier aux recommandations du "PEP8 style guide for python code"
Script ci dessous utilisé pour tester plus en detail QualGenerator
In [ ]:
from QualGenerator import QualGenerator
from matplotlib import pyplot
def qtester(length, quality, iteration, title):
q = QualGenerator(length, quality)
fig = pyplot.figure(figsize=(length/2, 10), dpi=100)
pyplot.title(title)
pyplot.ylabel('PHRED Quality')
pyplot.xlabel('Position')
position = [i+1 for i in range (length)]
for i in range (iteration):
q_list = q.random_qual_string()
pyplot.plot(position, q_list)
fig.savefig(title+'.png')
Exemple avec qtester(75, "good", 10, "len=75_qual=good_divider=2.png
In [ ]:
from ReferenceGenome import ReferenceGenome
from SlicePicker import SlicePicker
from QualGenerator import QualGenerator
def FastqGenerator (source, file, quality, length, mutfreq):
g = ReferenceGenome (source, file)
print(repr(g))
s = SlicePicker ()
print(repr(s))
q = QualGenerator (length, quality)
print(repr(q))
read = s.pick_single(g,length,1,1,mutfreq)
read.letter_annotations["phred_quality"] = q.random_qual_string()
read.id = "{}|{}|loc_{}_{}_{}".format(
read.annotations["source"],
read.annotations["refseq"],
read.annotations["location"][0],
read.annotations["location"][1],
read.annotations["orientation"])
print(read)
print(read.format("qual"))
print(read.format("fastq-illumina"))
return read
Pendant l'ecriture du module, il m'a semblé que la direction que j'avais prise auparavant pour SlicePicker (à savoir une classe portant les methodes pour single et pair read picking) était mauvaise. Il y a en effet beaucoup de paramètres à passer par FastqGenerator.
Séparation de la classe en 2 sous classes SlicePickerSingle et SlicePickerPair chacune spécialisée dans les read single ou paired. Le code est mieux factorisé et lors de l'instanciation les classes auront directement plusieurs variables stockées. Pour le moment j'ai un problème à l'instanciation via init de la super classe. A revoir
In [ ]:
# Running this code will raise `TypeError: must be type, not classobj`
class class_A:
def __init__( self ):
print 'Class A'
class class_B( class_A ):
def __init__( self ):
print 'Class B'
super( class_B, self ).__init__()
class_B()
# Running this code will no longer raise `TypeError: must be type, not classobj`
class class_A( object ): # < Note the inheritance of `object`.
def __init__( self ):
print 'Class A'
class class_B( class_A ):
def __init__( self ):
print 'Class B'
super( class_B, self ).__init__()
class_B()
Aprés un peu de débogage les classes SlicePicker, SlicePickerSingle et SlicePickerPair fonctionnent toutes 3. Voir code de la fonction draft de FastqGenerator pour les paires
In [ ]:
from ReferenceGenome import ReferenceGenome
from SlicePicker import SlicePickerPair
from QualGenerator import QualGenerator
def FastqGenerator (source, file, quality, read_len, sonic_min, sonic_max, sonic_mode, sonic_certainty, repeats, ambiguous, mut_freq):
g = ReferenceGenome (source, file)
print(repr(g))
s = SlicePickerPair(read_len, sonic_min, sonic_max, sonic_mode, sonic_certainty, repeats, ambiguous, mut_freq)
print(repr(s))
q = QualGenerator (read_len, quality)
print(repr(q))
read1, read2 = s.pick_slice(g)
read1.letter_annotations["phred_quality"] = q.random_qual_string()
read2.letter_annotations["phred_quality"] = q.random_qual_string()
read1.id = "R1"
read2.id = "R2"
print(read1)
print(read2)
print(read1.format("fastq-illumina"))
print(read2.format("fastq-illumina"))
return (read1, read2)
TO DO pour ce package
Reformatage du code pour être conforme à la norme PEP8
Test plus avancé des bogues possibles
Le nommage des id des reads sera fait désormais au fur et à mesure et pas uniquement dans SlicePicker. Nommage partielle dans SlicePickerPair (pas besoin dans SlicePickerSingle)
Suite formatage PEP8
Formatage PEP8
Ajout de commentaires dans la code
Simplification de certaines fonctions et début du nommage des id des SeqReccord.
PEP8 formatting and modifications of variables
A noter qu'il existe déjà une classe de base dans python dédié au parsing des fichiers de configuration.
Les 2 classes que j'ai codé semblent cependant plus adapté à mon programme, mais il faudrait peut être voir si ce n'est pas adaptable qd même.
In [ ]:
from ReferenceGenome import ReferenceGenome
from ReferenceJunctions import ReferenceJunctions
from SlicePicker import SlicePickerSingle
from SlicePicker import SlicePickerPair
bacteria = ReferenceGenome ("datasets/Bacterial_backbone.fa")
helper = ReferenceGenome ("datasets/Helper_plasmid.fa")
junction = ReferenceJunctions(50, 200, 10, bacteria, helper, 1, 1)
slicer1 = SlicePickerSingle(200,1,1,0.1)
slicer2 = SlicePickerPair(200, 300, 700, 350, 10, 1,1,0.1)
s1b = slicer1.pick_slice (bacteria)
print "S1B\n{}\n".format(s1b)
s1j = slicer1.pick_slice (junction)
print "S1J\n{}\n".format(s1j)
s2b = slicer2.pick_slice (bacteria)
print "S2B\nR1\n{}\nR2\n{}".format(s2b[0],s2b[1])
s2j = slicer2.pick_slice (junction)
print "S2J\nR1\n{}\nR2\n{}".format(s2j[0],s2j[1])
Retour de paires sans séquence parfois pour tirage dans la référence jonction. Cela venait du fait que la junction source était trop petite en cas de grande taille de sonication. Taille minimale des jonctions = 2 taille read si single end, mais = 2 sonic_max si paired ends... Attention a faire attention dans le futur.
Script de test pour voir si plus de read vide générés après n tentatives.
In [ ]:
def test(n):
for i in range (n):
s2 = slicer2.pick_slice (junction)
if len(s2 [0]) == 0:
print "Exit at iteration {}\n".format(i)
print("S2\nR1\n{}\nR2\n{}".format(s2[0],s2[1]))
return
print "All valid"
return
Comme pour SlicePicker auparavant, j'ai splitté le code en 2 FastqGeneratorSingle pour les read Single et FastqGeneratorPair pour les read paired. Les classes sont fonctionnelles
Les objets seq record possèdent désormais un compteur qui s'incrément à chaque fois qu'un read est pioché dans la référence.
Etant donnée que l'écriture du rapport va être plus simple maintenant, FastqGenerator va se charger de l'ecriture de n read dans le(s) fichiers de sortie. Apport des modifications nécessaires aux classes en question
Premier est grandeur nature pour les read single
In [ ]:
from ReferenceGenome import ReferenceGenome
from SlicePicker import SlicePickerSingle
from QualGenerator import QualGenerator
from FastqGenerator import FastqGeneratorSingle
helper = ReferenceGenome ("datasets/Helper_plasmid.fa")
slicer1 = SlicePickerSingle(200,1,1,0.1)
qualgen = QualGenerator (200, "good")
fast1 = FastqGeneratorSingle (slicer1, qualgen, "fastq-sanger")
fast1.write_fastq(helper, 1000, "test")
for name, rec in helper.d.items():
print "id = {}, nb samp = {}, length = {}".format(name, rec.annotations["nb_samp"], len(rec))
Affichage de sortie
Finalement j'ai remis le "name" de la source qui est en fait important pour nommer les fichiers exportés J'ai ajouté une méthode dans chacun des classe Reference pour écrire un rapport des références avec le nombre d'échantillonnage dans chacune. La fonction est différente dans ReferenceGenome et ReferenceJunctions.
Suppression de l'explication de _random_slice du code de ReferenceJunctions
Example of the strategy with size = 20 and min chimeric = 6
Junction -----------------------------|-----------------------------
Chimeric bases /////////////
Start area ///////////////ooooooo////////////////////////////////////
End area ////////////////////////////////////ooooooo////////////////Ecriture des fastq directement dans une archive gz Programme fonctionnel sans utilisation de la classe de configuration, mais lent.
Fin de reformatage de IsisConf. Desormais, la classe est plus simple et affiche seulement un message d'erreur si une option est mal paramétrée. La classe ConfFileParser a été supprimé et est maintenant remplacée par la classe de la librairy standard python ConfigParser. La gestion des exception a été ammelioré par la création d'une classe d'exception personalisée IsisConfException dont la fonction est de centralisé toute les exception possible qui peuvent intervenir pendant le parsing des options.
Reformatage en profondeur du fichier de configuration
#########################################################################################
# #
# Isis.py configuration file #
# Values can by customized with users values, but the file template must #
# remain unchanged #
# * = Changes of value are not recommended for proper programm execution #
# #
#########################################################################################
[General]
### General parameters
# Overall quantity of reads (min = 1, INTEGER)
read_num : 100000
# Read length (min = 1, INTEGER)
read_len : 150
# Frequency of mutations in final reads (min = 0, max = 1, FLOAT)
mut_freq : 0.01
# Allow sampling in repeat regions identified by lowercase characters (BOOLEAN)
repeats : True
# Allow sampling in ambigous IUPAC DNA bases (BOOLEAN)
ambiguous : False
[Frequency]
### Relative frequencies of DNA source in fastq (sum of frequencies should be equal to 1)
# Host genome frequency (min = 0, max = 1, FLOAT)
freq_host : 0.49
# Virus genome frequency (min = 0, max = 1, FLOAT)
freq_virus : 0.49
# True junctions frequency (min = 0, max = 1, FLOAT)
freq_tjun : 0.01
# False junctions frequency (min = 0, max = 1, FLOAT)
freq_fjun : 0.01
[Junction]
### Parameters specific to junctions between host and viral DNA
# Minimal number bases of bases from each references in reads (min = 0, max = read size/2 for se mode, max = read_len for pe mode, INTEGER)
min_chimeric : 25
# Mean number of sampling in true junctions (min = 0, max = freq_tj * read_num, FLOAT)
samp_tjun : 20
# Mean number of sampling in false junction (min = 0, max = freq_tj * read_num, FLOAT) *
samp_fjun : 0.1
[Sonication]
## Parameters of fragment sonication distribution to set up if paired end mode
# Minimal sonication size (min = read_len + min_chimeric, max = size of reference sequences, INTEGER)
sonic_min : 200
# Minimal sonication size (min = sonic_min, max = size of reference sequences, INTEGER)
sonic_mode : 400
# Minimal sonication size (min = sonic_min, max = size of reference sequences, INTEGER)
sonic_max : 1000
# Certainty of the sonication smire (min = 5 (wide peak), max = 50 (thin peak), INTEGER) *
sonic_certainty : 8
[Quality]
# Quality score scale of fastq output file (fastq-sanger OR fastq-solexa OR fastq-illumina) *
qual_scale : fastq-sanger
# Quality range for mean Phred scores (very-good OR good OR medium OR bad OR very-bad) *
qual_range : goodEcriture d'un petit script de test de la classe
In [ ]:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from IsisConf import IsisConf, IsisConfException
def main ():
"""main function"""
# Instantiate the configuration file parser/verifier
try:
config = IsisConf()
except IsisConfException as E:
print (E)
exit (0)
print (repr (config))
if __name__ == '__main__':
main()
Resultats execution
./Isis.py -V Bacterial_backbone.fa -H Helper_plasmid.fa -C Conf.txt -o Test -p
Parsing and verification of config file with IsisConf
<Instance of IsisConf>
ambiguous : False
basename : Test
conf_file : Conf.txt
host_genome : Helper_plasmid.fa
min_chimeric : 25
mut_freq : 0.01
nread_fjun : 1000
nread_host : 49000
nread_tjun : 1000
nread_virus : 49000
pair : True
qual_range : good
qual_scale : fastq-sanger
read_len : 150
read_num : 100000
repeats : True
sonic_certainty : 8
sonic_max : 1000
sonic_min : 200
sonic_mode : 400
uniq_fjun : 10000
uniq_tjun : 50
virus_genome : Bacterial_backbone.faModification du main pour ajouter une recupération initiale de tout les paramètres utiles dans dans variables locales à partir de IsisConf
Un dossier contenant les fichiers générés par un run de test avec l'ADN genomique mm1O et le génome AAV et disponible ICI. Le temps d'execution pour les 100 000 paires de read demandé est raisonable (5 minutes) mais pourrait encore être ammelioré considérablement.
En ce qui concerne la sortie permettant de verifier l'alignement sur chaque jonctions, le mieux serait de produire un fichier fasta du dictionnaire de jonctions (à implementer dans la classe Referencejunction) et un fichier SAM regroupant les reads samplés dans la jonctions en questions à partir de ça il sera simple de faire une visualisation sous igv.
En reflechissant bien il serait plus pertinant de generer un graph de couverture directement dans le programme. Il me semblerait même être bien d'ajouter une option dans le fichier de configuration qui autorise on non la génération des graphiques de distribution de taille des fragments et de couverure autour des jonctions
La localisation des reads est maintenant donnée de façon precise pour les read issu de pair (au lieu de la position du fragment dont ils sont tirés.
La localisation est desormais toujours donnée en 5' vers 3' et l'orientation est indiquée par + ou -
Un programme tirant parti de matplotlib trace un graph de couverture à partir d'une liste indiquant la couverture à chaque position = jun_cov_graph
Les programmes write_fastq_single et write fastq_pair ont été modifié. Ils accepte maintenant une liste de source contenant : L'objet source, le nb de read a sampler et un booléen indiquant si il faut ajouter les paires au graph de couverture (True pour junctions uniquement).
Implementation de 2 options supplémentaires dans le fichier de configuration:
Si min_chimeric est faible une grande latitude est laissé dans le choix de la position au dessus de la jonction. Au contraire un min_chimeric au max (read_len / 2) autorise une seule positionjuste sur le milieu de la jonction. # min_chimeric indique le minimal de base sur l'une des 2 references.#
On observe le même phenomene, mais il y a en plus une variabilité de taille maximale causée par la distribution beta de la taille des fragments. Dans ce mode, on peut monter à max = read_len, dans ce cas, il y un read d'un coté de la jonction, et le mate toujours de l'autre coté.
To simplify id naming will be done just before fastq writting from the annotation dictionnary
Exemple of origin reference sequence retrieval t
In [ ]:
s.annotations["source"] = self
In [34]:
from ReferenceGenome import ReferenceGenome as RefGen
a = RefGen("test", "../test/data/Bacterial_backbone.fa")
b = a.get_slice(100)
print(b.annotations)
str(b.annotations["source"].d[b.annotations["refseq"]].seq[b.annotations["location"][0]:b.annotations["location"][1]])
Out[34]:
In [1]:
# VERSION SINGLE END
from ReferenceGenome import ReferenceGenome as RefGen
a = RefGen("test", "../test/data/Bacterial_backbone.fa")
from SlicePicker import SlicePickerSingle
b = SlicePickerSingle(100,1,1,0)
c = b.pick_slice(a)
print (c)
# Sequence stored in the object
print(str(c.seq))
# Sequence retrieved through annotation dictionnary
seq = c.annotations["source"].d[c.annotations["refseq"]].seq[c.annotations["location"][0]:c.annotations["location"][1]]
if c.annotations["orientation"] =='+':
print (str(seq))
else:
print (str(seq.reverse_complement()))
In [1]:
# VERSION PAIR END
from ReferenceGenome import ReferenceGenome as RefGen
a = RefGen("test", "../test/data/Bacterial_backbone.fa")
from SlicePicker import SlicePickerPair
b = SlicePickerPair(100,150,200,300,10,1,1,0)
c,d = b.pick_slice(a)
print (c)
print (d)
# Sequence stored in the object
print(str(c.seq))
# Sequence retrieved through annotation dictionnary
seq = c.annotations["source"].d[c.annotations["refseq"]].seq[c.annotations["location"][0]:c.annotations["location"][1]]
if c.annotations["orientation"] =='+':
print (str(seq))
else:
print (str(seq.reverse_complement()))
# Sequence stored in the object
print(str(d.seq))
# Sequence retrieved through annotation dictionnary
seq = d.annotations["source"].d[d.annotations["refseq"]].seq[d.annotations["location"][0]:d.annotations["location"][1]]
if d.annotations["orientation"] =='+':
print (str(seq))
else:
print (str(seq.reverse_complement()))
Possibility to recover the seq directly thought the dictionnary of annotations. Good but much heavier + do not take into acount possible mutations.
In [ ]:
def origin_coord (self, refseq, start, end):
""" Return a string describing the the origin of a sequence from a
junction of the dictionnary
"""
# If the give coord overlap only the left reference of a junction
if start <= end <= self.half_len:
return ("{}-{}={}".format(
1, end-start,
self.coord_to_str(self.d[refseq].annotations, "ref1", start, end)))
# If the give coord overlap only the right reference of a junction
elif self.half_len <= start <= end:
return ("{}-{}={}".format(
1,end-start,
self.coord_to_str(self.d[refseq].annotations, "ref2", start-self.half_len, end-self.half_len)))
# If the give coord overlap both references of a junction
else :
return ("{}-{}={}|{}-{}={}".format(
1,self.half_len-start,
self.coord_to_str(self.d[refseq].annotations, "ref1", start, self.half_len-1),
self.half_len-start+1, end-start,
self.coord_to_str(self.d[refseq].annotations, "ref2", self.half_len, end)))
def coord_to_str (self, d, ref, start, end, ):
"""
"""
ref_id = d[ref+"_refseq"].id
if d[ref+"_orientation"] == '+':
ref_start = d[ref+"_location"][0] + start
ref_end = d[ref+"_location"][0] + end
else:
ref_start = d[ref+"_location"][1] - end
ref_end = d[ref+"_location"][1] - start
return ("{}:{}-{}".format(ref_id, ref_start, ref_end))
In [1]:
from ReferenceGenome import ReferenceGenome as RefGen
from ReferenceJunctions import ReferenceJunctions as RefJun
a = RefGen("Bact", "../test/data/Bacterial_backbone.fa")
b = RefGen("Help", "../test/data/Helper_plasmid.fa")
c = RefJun("Tjun", 50, 200, 10, a, b, 1, 1)
In [2]:
refseq = c.d["Junction_05"]
print ("\nJUNCTION SEQUENCE")
print (str(refseq.seq[180:219]))
print ("\nJUNCTION CARACT")
for ref in ["ref1", "ref2"]:
print ("{} : {} - {} ({})".format(
refseq.annotations[ref+"_refseq"].id,
refseq.annotations[ref+"_location"][0],
refseq.annotations[ref+"_location"][1],
refseq.annotations[ref+"_orientation"]))
print ("\nCOORDINATES")
# left ref only
print(c.origin_coord(refseq, 180, 199))
# right ref only
print(c.origin_coord(refseq, 200, 219))
# Overlap both ref
print(c.origin_coord(refseq, 180, 219))
@virus|018|0|All=Bacterial_backbone:810-960
gatgcgcgagctagctagctacgatcgtagctagtcatgcgtagcttagccgatcgccgataatatcgtacgtatgctacgagcagtctagctacgcctagagtgctacgtatcatctattgcgcatatatagccgaatcagatagttat
+
B@?@A?A@BFFEEEEGFFFFFFFFFEEECDEEEEEDDEEEEEEEEFFFGFFEDEEEFFFFFEEEEEEDFDEFGHGGCCGFEEDDECDFHE@EGFEEEBFC?EFC:>;>E??;=F?CIHFF<?EGFDCB@AD=C>9BE8<>>?;>?GHF:;
...
@host|063|0|All=Phage_lambda_genome_J02459.1:21013-21163
AACCTGAGTCAGTTCAGTCAGGCTGGCGGCATCATTTTCCGCAAGATACGGTAATTTATTTTTCACCGTGGAAAGCCCTGCCAGCGCCGTCAGTGTCGCATTCTTCGGTTGTTTACCCGCAAGCGCGTTAGTCATGGTGGTAGCAAAATC
+
ACCBCA@B?ABBDCDGGGDFGFFFFFEFFFFFGFFFFFGFFDEEEEEEGGEBDFEEDFFGFFFGFGEFFIHIFEEEFECDECCDBFFHHGFIHFEACFGD@A<@BDFIACGCACCABA@>B>@AFDDA?C?6A=?C@<?@?<>E=DBAAC
...
@True_Junction|112|1|0-150=Helper_plasmid:19657-19808
TACGTCACTTCCCATTTTAAGAAAACTACAATTCCCAACACATACAAGTTACTCCGCCCTAAAACCTACGTCACCCGCCCCGTTCCCACGCCCCGCGCCACGTCACAAACTCCACCCCCTCATTATCATATTGGCTTCAATCCAAAATAA
+
EGB@ABBDEDDDBAADDEEEEDGFDEEEDEEFFFFFFFFFFFEEEFFEEFFFFFGEEEEEEFFGFGFGFCDEBDGDDEDFFFGECBFGFFFFEEEEBC:<?BDFDGFFDCDCAC?;=?CCBA8=D@@ACB=9C?EF>9A?BE47966H7?
...
@False_Junction|249|1|0-150=Bacterial_backbone:515-666
agtctagctacgcctagagtgctacgtatcatctattgcgcatatatagccgaatcagatagttatatgcgcgcgatgcagtcgtactgagctagcgtactgacgttactactgagtcgactgatagatgcgcgagctagctagctacga
+
AAAAA?A@A>??CDCDEEFFDECCCDDEEEEDEEEFFFGFDFEFEGGFGFFFFDEEDFHGFFFDEEEEEFHIEHGGFFIDGFHGFEDEEAA<>CFFDDCGIIBHF@HFEFC>;EGACGCCEHECAB<C<@DD9AI>=??DBIGD<5::93Exemple of reads extracted from R1.fastq.gz with condition that allow reads to overlap a junction
...
@True_Junction|194|1|0-97=Bacterial_backbone:1144-1242|98-150=Phage_lambda_genome_J02459.1:41117-41170
CatcgtagctagctGgctcgcgcatctatcagtcgactcagtagtaacgtcagtacgctagctcagtacTactgcatcgcgcgcatataactatctgaCAGCATATTTGCTCTGGCCAATGGAGCAAAAGCGACGGGCAGGTAAAGACGT
+
A@@AABBB@ABCDDCDEEEDEEEEEEECEEEFFFFHGFFFEEFFGHGFFFFFEEEEEEEDDEEEEDEEGFGFFHGDEECCBDFFEFFFBHGB<?ECGFEEIGHHIIFB>HFDGC=<ACBHIAB@A>:C8EG=@83:<IAA?GDBFG><>
...
@False_Junction|107|1|0-4=Bacterial_backbone:1019-1024|5-150=Phage_lambda_genome_J02459.1:35573-35719
TCAGATCTCTCACCTACCATACAATGCCCCCCTGCAAAAAATAAATTCATATAAAAAACATACAGATAACCATCTGCGGTGATAAATTATCTCTGGCGGTGTTGACATAAATACCACTGGCGGTGATACTGAGCACATCAGCAGGttata
+
GECEFDDDBCBEEDDFFFFFFFDEFFGFGFEDEEEEEEEEDEEEEFFGFFFFGEEEFFFFFEEEECDEFHGGFFIIGFFHDEEBCDGGFFFHFEDDCHIGCCE@=:A@AE?DDDDB?=BHEA@?@<?GFAA?<>>CDEB9<@CGD:8FIEMain() was removed to be replaced by top level instructions just bellow if __name__ == '__main__':
It has several avantage over the use of an independant main :
Make a doxygen complient documentation