Discussion :

[benoit3166]

Bonjour à tous,

Je suis en train de me réorienter en bioinformatique. En plein apprentissage de python, j'ai decide de passer un peu à la pratique. J'ai fait un script qui ouvre un fichier au format fasta. Il s'agit d'un fichier texte dans lequel on stocke des sequences nucléique (ATCG). Un fichier peut contenir une séquence nucléique ou plusieurs sequences.

Voici un exemple de fichier fasta avec plusieurs entrées :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
>gi|373251181|ref|NG_001742.2| Mus musculus olfactory receptor GA_x5J8B7W2GLP-600-794 (LOC257854) pseudogène on chromosome 2
AGCCTGCCAAGCAAACTTCACTGGAGTGTGCGTAGCATGCTAGTAACTGCATCTGAATCTTTCAGCTGCT
TGTTGGGCCTCTCACAAGGCAGAGTGTCTTCATGGGACTTTGATATTTATTTTTGTACAACCTAAGAGGA
ACAAATCCTTTGACACTGACAAATTGGCTTCCATATTTTATACCTTAATCATCTCCATGTTGAATTCATT
GATCAACAGTTTAAGAAAAAAAGATGTAAAAATGCTTTTAGAAAGAGAGGCAAAGTTATGCACAATAACT
TCTCATGAAGTCACAGTTTGTTAAAAGTTGCCTTAGTTCACAATAAATAATTATGTATGCTCTATAATTT
CAGTGA
 
>gi|425889810|gb|JZ141974.1| ES4f_H12.esd SL-enriched library from schistosomula Schistosoma mansoni cDNA, mRNA sequence
ATGAAGTGAAAAATTTATTTATTTATAACAAAATGATGAATCAATGTACATTGTATTTATCTAATGTAGA
TATAGCCTAACTTCCCGCTTTTTTCAATGAGCGTATATTGCGTAATTCACGATATAATACATTTTGTTGA
ACAAGTGGATCAAATGCAACTTTCTCTTTACGATTATTAGCTAGACGAGGTCGAAAACCAACTACACAAT
GTCCTGAACCACTGGTACTTTCTAATAATACCATAACAGTTTTTATAAATCCTTGTGATGTTGACATAAA
GGATAAATGTGTATACTTGTACAAAGTTATTTATATTGGGTAATCTTCGTTGATATTATGTAATGATCAT
GCAACAAATCACAAGAGTAAAACCGTGACGG

Mon script (voir ci-dessous) a pour but de :

1. Compter les différents acide nucléique pour chaque sequence.
2. Enregistrer les résultats dans un fichier csv ou txt avec le nom de la sequence, le nombre de chaque acide nucléique.

Mon script lit bien les fichiers fasta et il créer bien le fichier que je lui demande. Seul problème, lorsque j'ouvre le tableau (csv ou txt) à l'aide de libre office ou excel. Les résultats sont décalés. C'est-à-dire que le nom de la séquence apparait bien dans la premiere ligne et à la premiere colonne. Mais les résultats associé à cette sequences sont écrite à partir de la seconde ligne et à la seconde colonne.

Comment résoudre ce problème ?

Petite question supplémentaire: comment faire un header ?

Voici le script que j'ai écris:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
 
#!/usr/bin/python3                                                                                                                                                                                                             
 
import sys, os
 
pyScriptName = sys.argv[0]
dataFileName = sys.argv[1]
outputFileName = sys.argv[2]
 
# Read fasta file                                                                                                                                                                                                              
fasta_dict = {}
 
with open(dataFileName, "r") as fasta_file, open(outputFileName, "w") as output_cpg:
    sequence_id = ""
    for line in fasta_file:
        if line.startswith(">"):
            sequence_id = line
            fasta_dict[sequence_id] = ""
        else:
            fasta_dict[sequence_id] += line.strip()
 
    for id in fasta_dict:                                                                                                                                                                                                     
        seq = fasta_dict[id]
        G_count = seq.count("G")
        C_count = seq.count("C")
        A_count = seq.count("A")
        T_count = seq.count("T")
        CpG_count = seq.count("CG")
        output_cpg.write(id +  ";" + str(G_count) + ";" + str(C_count) + ";" + str(A_count) + ";" + str(T_count) + ";" + str(CpG_count) + ";" + "\n")                                                
 
fasta_file.close()
output_cpg.close()

Merci pour votre aide.

A bientôt !

[wiztricks]

Salut,

Comment résoudre ce problème ?

En regardant le fichier de sortie pour confirmer que vous écrivez 2 lignes au lieu d'une... Puis en corrigeant le code pour ne pas écrire cette fin de ligne en trop.

Petite question supplémentaire: comment faire un header ?

Le header est juste la première ligne du fichier.

- W

[benoit3166]

Bonjour wiztricks,

Merci pour ta réponse.

Salut,

En regardant le fichier de sortie pour confirmer que vous écrivez 2 lignes au lieu d'une... Puis en corrigeant le code pour ne pas écrire cette fin de ligne en trop.

Avant d'ouvrir la discussion, j'avais déjà fait plusieurs essais comme vous me le conseillez, afin d'identifier d'où vient la ligne supplémentaire. Elle vient du nom des sequences. En effet, j'ai obtenu un tableau sans décalage lorsque j'ai écris ce code :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
output_cpg.write(str(G_count) + ";" + str(C_count) + ";" + str(A_count) + ";" + str(T_count) + ";" + str(CpG_count) + ";" + id)

Or, ce que je cherche à avoir dans le fichier de sortie, c'est que le nom des séquences soit dans la première colonne. C'est plus lisible pour analyser les résultats à l'oeil dans un premier temps. Il y a donc un saut de ligne dans le nom des séquences. Quand j'utilise cette ligne de code (voir ci-dessous), j'obtiens un fichier comme je l'ai décris dans mon premier post: le nom des séquence décale tout.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

output_cpg.write(id +  ";" + str(G_count) + ";" + str(C_count) + ";" + str(A_count) + ";" + str(T_count) + ";" + str(CpG_count) + ";" + "\n")

Comment enlever ce saut de ligne dans le nom des sequences (qui apparait après le > dans le fichier au format fasta) pour avoir un fichier de sortie interpretable ?

Le header est juste la première ligne du fichier.

- W

Oui et ce que je cherche à faire, c'est rajouter une ligne informative qui informe sur le contenu de chaque colonne (ex: nom des sequence, %A, %T, %C, %G, etc.). Cette ligne sera utile pour des analyses statistique avec R (header = TRUE avec la fonction read.table ou read.csv) ou avec d'autres libraries en python (quand je le maitriserai mieux).

J'aurai du être plus explicite dans mon premier post. Je m'en excuse.

[wiztricks]

Salut,

Comment enlever ce saut de ligne dans le nom des sequences (qui apparait après le > dans le fichier au format fasta) pour avoir un fichier de sortie interpretable ?

Le code fait déjà cela pour les lignes donnant la séquences de nucléotides.
Si c'est vous qui l'avez écrit, vous savez déjà comment!

ce que je cherche à faire, c'est rajouter une ligne informative qui informe sur le contenu de chaque colonne (ex: nom des sequence, %A, %T, %C, %G, etc.). Cette ligne sera utile pour des analyses statistique avec R (header = TRUE avec la fonction read.table ou read.csv) ou avec d'autres libraries en python (quand je le maitriserai mieux).

Le code écrit déjà plein de lignes dans le fichier, je ne vois pas quelle difficulté vous avez à écrire une première ligne.

- W

[benoit3166]

Salut,

Merci wiztricks pour ton aide. J'ai résolu le problème et cela fonctionne. J'avais oublié le strip sur sequence_id.

Pour ceux que cela intéresse, voici le code final ci-dessous.

A bientôt !

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
 
#!/usr/bin/python3                                                                                                                                                                                                                                                 
 
import sys, os
 
pyScriptName = sys.argv[0]
dataFileName = sys.argv[1]
outputFileName = sys.argv[2]
 
# Read fasta file                                                                                                                                                                                                                                                  
fasta_dict = {}
 
with open(dataFileName, "r") as fasta_file, open(outputFileName, "w") as output_cpg:
    sequence_id = ""
    for line in fasta_file:
        if line.startswith(">"):
            sequence_id = line.strip()
            fasta_dict[sequence_id] = ""
        else:
            fasta_dict[sequence_id] += line.strip()
 
    output_cpg.write("Name" + ";" + "Size" + ";" + "G" + ";" + "C" + ";" + "A" + ";" + "T" + ";" + "CpG" + ";" + "\n")
    for id in fasta_dict.keys():
        print(id)
        win = len(fasta_dict[id])
	seq = fasta_dict[id]
        G_count = seq.count("G")
	C_count = seq.count("C")
        A_count = seq.count("A")
	T_count = seq.count("T")
        CpG_count = seq.count("CG")
 
        output_cpg.write(id + ";" + str(win) +  ";" + str(G_count) + ";" + str(C_count) + ";" + str(A_count) + ";" + str(T_count) + ";" + str(CpG_count) + ";" + "\n")
fasta_file.close()
output_cpg.close()