Bonjour,
J'aurai besoin d'aide pour l'analyse de fichiers fasta (sequence ADN).
J'ai une centaine de fichiers fasta contenant chacun 600 sequences.
Voila a quoi ressemble les fichiers:
>1-1-1
CAACCCACAAAAACCCAACACAACAAAACCAACCCAACCAACCCCCCAACACACCCAAAACACACACAACCCAACAAACCACAAAACCAAACAACCCAACACACAACACCCCACCAACCAACAACACCCAAACCAACCCAAAACAAAACCACCACCACCACCCCACACAAAACCAAACCCACC
>1-1-2
CAACCCACAAAAACCCAACAAAACAAACCCCACCCAACCAACCCCCCACCACACCCAAAACACACACAACCCAACAAACCACAAAACCAAACAAACCAACACACAACACCCAACCAACCAACAACACCCAAACCAACCCAAAACAAAACCACCACCACCACCCCACACAAAACCAAACCCACC
>1-1-3
AAACCCACAAAAACCCAACACAACACAACCACCACAACCAACCCCCCCCCACACCCAAAACACACCCAAAACAACAAACCACAACACCAAACAAACCAAAACACAACACCCCCCCAACCAACACCACCAAAACCAACCCAAAACAACCCCACCACCACCACCCCAAACAAAACCACACCCACC
>1-1-4
CAACCCACAAAAACCCAACACAACAAAACCAACCCAACCAACCCCCCAACACACCCAAAACACACACAACCCAACAAACCACAAAACCACACAACCCAACACACAACACCCCACCAACCAACAACACCCAAACCAACCCAAAACAAAACCACCACCACCACCCCACACAAAACCAAACCCACC
J'aimerai calculer pour chaque fichier de 600 sequences le nombre de differences par pair de sequences.
Par example, le nombre de difference entre 1-1-1 et 1-1-2 puis le nombre de difference entre 1-1-1 et 1-1-3, etc
A la fin, j'aimerai obtenir un fichier output qui me donne le nombre de pairs pour chaque nombre differences ainsi que le pourcentage.
Par example:Nombre de Difference Nombre de pairs Pourcentage
1 25 2%
2 50 4%
... .... ....
Voila le script que j'ai ecrit:
Code :
- import sys
- def readFas(inf):
- seqs = []
- count = 0
- for line in inf:
- if line.find(">" ) >= 0:
- #seq = line
- count += 1
- else:
- #seq += line
- if count == 0:
- return seqs
- def count_Diff(seqs, real):
- numDiff = 0
- for s in range(len(seqs[0])):
- for idx_seq in range(len(seqs)-1):
- if seqs[idx_seq][s] != seqs[idx_seq+1][s]:
- numDiff += 1
- break
- return numDiff
- #################################################
- filehead = sys.argv[1]
- startRepl = int(sys.argv[2])
- endRepl = int(sys.argv[3])
- real = 0
- outf = open("p-seg1_"+filehead+".txt", "w" )
- for sim in range(startRepl, endRepl):
- for repl in range(1, 51):
- name = "seg1_"+ filehead + "_" +str(sim) + "_" + str(repl) + "_pick.fas"
- inf = open(name, "r" )
-
- seqsWindow = readFas(inf)
- inf.close()
-
- Diff = count_Diff(seqs, real)
- outf.write(str(Diff)"\n" )
|
Cependant je n'arrive pas a faire un script qui fonctionne et je ne vois vraiment pas comment obtenir le tableau que je veux.
Si quelqu'un pouvait m'aider ou avait des suggestions a me proposer, j'en serai tres heureuse car je suis vraiment perdue.
Merci d'avance
