Создать Dendogram из Геном

Question

Я хотел играть с геномных данных:

Species_A = ctnngtggaccgacaagaacagtttcgaatcggaagcttgcttaacgtag 
Species_B = ctaagtggactgacaggaactgtttcgaatcggaagcttgcttaacgtag 
Species_C = ctacgtggaccgacaagaacagtttcgactcggaagcttgcttaacgtag 
Species_D = ctacgtggaccgacaagaacagtttcgactcggaagcttgcttaacgccg 
Species_E = ctgtgtggancgacaaggacagttccaaatcggaagcttgcttaacacag 

Я хотел создать дендрограммы на основе того, насколько близко эти организмы связаны друг с другом, учитывая последовательность генома выше. То, что я сделал первым было подсчитать число элементов а, C-х, Т х и г-х каждого вида, то я создал массив, а затем нанесены на дендрограммы:

gen_size1 = len(Species_A) 
a1 = float(Species_A.count('a'))/float(gen_size1) 
c1 = float(Species_A.count('c'))/float(gen_size1) 
g1 = float(Species_A.count('g'))/float(gen_size1) 
t1 = float(Species_A.count('t'))/float(gen_size1) 
. 
. 
. 
gen_size5 = len(Species_E) 
a5 = float(Species_E.count('a'))/float(gen_size5) 
c5 = float(Species_E.count('c'))/float(gen_size5) 
g5 = float(Species_E.count('g'))/float(gen_size5) 
t5 = float(Species_E.count('t'))/float(gen_size5) 

my_genes = np.array([[a1,c1,g1,t1],[a2,c2,g2,t2],[a3,c3,g3,t3],[a4,c4,g4,t4],[a5,c5,g5,t5]]) 
plt.subplot(1,2,1) 
plt.title("Mononucleotide") 
linkage_matrix = linkage(my_genes, "single") 
print linkage_matrix 
dendrogram(linkage_matrix,truncate_mode='lastp', color_threshold=1, labels=[Species_A, Species_B, Species_C, Species_D, Species_E], show_leaf_counts=True) 
plt.show() 

Виды А и В являются вариантами одного и того же организма и I я ожидаю, что и то, и другое должно расходиться с общей формой клада корня, то же самое происходит с видами C и D, которые должны расходиться с другой общей кладой от корня, а затем с видами E, расходящимися от основного корня, потому что это не связано с видами A-D К сожалению, результат дендрограммы был смешан с видами A и E, расходящимися от общей клады, затем виды C, D и B в другой кладе (довольно перепутались).

Я читал об иерархической кластеризации для последовательности геномов, но я заметил, что он вмещает только 2-мерную систему, к сожалению, у меня есть 4 измерения, которые являются a, c, t и g. Любая другая стратегия для этого? Спасибо за помощь!

Answer 1

Это довольно распространенная проблема в биоинформатике, поэтому вы должны использовать библиотеку биоинформатики, такую ​​как BioPython, которая имеет эту функциональность.

Сначала вы создаете мульти-файл FASTA с вашими последовательностями:

import os 
from Bio import SeqIO 
from Bio.Seq import Seq 
from Bio.SeqRecord import SeqRecord 
from Bio.Alphabet import generic_dna 


sequences = ['ctnngtggaccgacaagaacagtttcgaatcggaagcttgcttaacgtag', 
      'ctaagtggactgacaggaactgtttcgaatcggaagcttgcttaacgtag', 
      'ctacgtggaccgacaagaacagtttcgactcggaagcttgcttaacgtag', 
      'ctacgtggaccgacaagaacagtttcgactcggaagcttgcttaacgccg', 
      'ctgtgtggancgacaaggacagttccaaatcggaagcttgcttaacacag'] 

my_records = [SeqRecord(Seq(sequence, generic_dna), 
       id='Species_{}'.format(letter), description='Species_{}'.format(letter)) 
       for sequence, letter in zip(sequences, 'ABCDE')] 

root_dir = r"C:\Users\BioGeek\Documents\temp" 
filename = 'my_sequences' 
fasta_path = os.path.join(root_dir, '{}.fasta'.format(filename)) 

SeqIO.write(my_records, fasta_path, "fasta") 

Это создает файл C:\Users\BioGeek\Documents\temp\my_sequences.fasta который выглядит следующим образом:

>Species_A 
ctnngtggaccgacaagaacagtttcgaatcggaagcttgcttaacgtag 
>Species_B 
ctaagtggactgacaggaactgtttcgaatcggaagcttgcttaacgtag 
>Species_C 
ctacgtggaccgacaagaacagtttcgactcggaagcttgcttaacgtag 
>Species_D 
ctacgtggaccgacaagaacagtttcgactcggaagcttgcttaacgccg 
>Species_E 
ctgtgtggancgacaaggacagttccaaatcggaagcttgcttaacacag 

Далее, используйте инструмент командной строки ClustalW, чтобы сделать выравнивание нескольких последовательностей:

from Bio.Align.Applications import ClustalwCommandline 
clustalw_exe = r"C:\path\to\clustalw-2.1\clustalw2.exe" 
assert os.path.isfile(clustalw_exe), "Clustal W executable missing" 
clustalw_cline = ClustalwCommandline(clustalw_exe, infile=fasta_path) 
stdout, stderr = clustalw_cline()  
print stdout 

Это печатает:

CLUSTAL 2.1 Multiple Sequence Alignments 


Sequence format is Pearson 
Sequence 1: Species_A  50 bp 
Sequence 2: Species_B  50 bp 
Sequence 3: Species_C  50 bp 
Sequence 4: Species_D  50 bp 
Sequence 5: Species_E  50 bp 
Start of Pairwise alignments 
Aligning... 

Sequences (1:2) Aligned. Score: 90 
Sequences (1:3) Aligned. Score: 94 
Sequences (1:4) Aligned. Score: 88 
Sequences (1:5) Aligned. Score: 84 
Sequences (2:3) Aligned. Score: 90 
Sequences (2:4) Aligned. Score: 84 
Sequences (2:5) Aligned. Score: 78 
Sequences (3:4) Aligned. Score: 94 
Sequences (3:5) Aligned. Score: 82 
Sequences (4:5) Aligned. Score: 82 
Guide tree file created: [C:\Users\BioGeek\Documents\temp\my_sequences.dnd] 

There are 4 groups 
Start of Multiple Alignment 

Aligning... 
Group 1: Sequences: 2  Score:912 
Group 2: Sequences: 2  Score:921 
Group 3: Sequences: 4  Score:865 
Group 4: Sequences: 5  Score:855 
Alignment Score 2975 

CLUSTAL-Alignment file created [C:\Users\BioGeek\Documents\temp\my_sequences.aln] 

my_sequences.dnd файл ClustalW создает, является стандартной Newick tree file и Bio.Phylo может разобрать эти:

from Bio import Phylo 
newick_path = os.path.join(root_dir, '{}.dnd'.format(filename)) 
tree = Phylo.read(newick_path, "newick") 
Phylo.draw_ascii(tree) 

который печатает:

 ____________ Species_A 
    ____| 
| |_____________________________________ Species_B 
| 
_|   ____ Species_C 
|_________| 
|   |_________________________ Species_D 
| 
|__________________________________________________________________ Species_E 

Или, если у вас есть matplotlib или pylab вы можете создать графику с помощью функции draw:

tree.rooted = True 
Phylo.draw(tree, branch_labels=lambda c: c.branch_length) 

, который производит:

Это дендрограммы наглядно иллюстрирует то, что вы наблюдали: что виды А и В являются вариантами одного и того же организма, и оба расходятся от общего клады от корня. То же самое происходит с видами C и D, и они расходятся из другой общей клады из корня. Наконец, виды E расходятся с основным корнем, потому что он не связан с видами от A до D.

Answer 2

Ну, используя SciPy, вы можете использовать пользовательскую дистанцию ​​(моя ставка на Needleman-Wunsch или Smith-Waterman в начале). Here - пример того, как играть с вашими входными данными. Вы также должны проверить how to define a custom distance in SciPy. После его установки вы можете использовать более продвинутый подход к выравниванию, например, MAFFT. Вы можете извлечь отношения между геномами и использовать их при создании своей дендрограммы.